Emergent Genome-Wide Control in Wildtype and Genetically Mutated Lipopolysaccarides-Stimulated Macrophages

Large-scale gene expression studies have mainly focused on highly expressed and ‘discriminatory’ genes to decipher key regulatory processes. Biological responses are consequence of the concerted action of gene regulatory network, thus, limiting our attention to genes having the most significant variations is insufficient for a thorough understanding of emergent whole genome response. Here we comprehensively analyzed the temporal oligonucleotide microarray data of lipopolysaccharide (LPS) stimulated macrophages in 4 genotypes; wildtype, Myeloid Differentiation factor 88 (MyD88) knockout (KO), TIR-domain-containing adapter-inducing interferon-β (TRIF) KO and MyD88/TRIF double KO (DKO). Pearson correlations computed on the whole genome expression between different genotypes are extremely high (>0.98), indicating a strong co-regulation of the entire expression network. Further correlation analyses reveal genome-wide response is biphasic, i) acute-stochastic mode consisting of small number of sharply induced immune-related genes and ii) collective mode consisting of majority of weakly induced genes of diverse cellular processes which collectively adjust their expression level. Notably, temporal correlations of a small number of randomly selected genes from collective mode show scalability. Furthermore, in collective mode, the transition from large scatter in expression distributions for single ORFs to smooth linear lines emerges as an organizing principle when grouping of 50 ORFs and above. With this emergent behavior, the role of MyD88, TRIF and novel MyD88, TRIF-independent processes for gene induction can be linearly superposed to decipher quantitative whole genome differential control of transcriptional and mRNA decay machineries. Our work demonstrates genome-wide co-regulated responses subsequent to specific innate immune stimulus which have been largely neglected.


Introduction
The innate immune system utilizes pattern-recognition receptors (PRRs), present in phagocytes such as macrophages and dendritic cells, to recognize pathogen associated molecular patterns (PAMPs), such as lipopolysaccarides (LPS). LPS, which is found on the outer membrane of Gram-negative bacteria, through the Toll-like receptor (TLR) 4, triggers a cascade of signaling events initiated mainly by the MyD88-and TRIFdependent pathways. This activates a number of common transcription factors including activator protein (AP)-1, nuclear factor-kB (NF-kB) and interferon regulatory factors (IRF)-3. As a result, a number of cytokines such as IL-1b and TNF-a, and type I interferons such as IFN-a and IFN-b are produced. These proinflammatory mediators activate helper T-cells for the onset of acquired immune defense where foreign intruders are eliminated and immunological memory is created [1][2][3][4][5]. These processes happen at a multi-cellular level implying the coordinated activities of many different tissues. Immunological responses are self-limiting, highly orchestrated systemic processes that if not precisely controlled can lead to major illnesses such as autoimmune diseases, cardiovascular diseases and cancer [6][7][8].
Recent high throughput experimental technologies have enabled the comprehensive analysis of cellular response to a given stimulus. However, molecular immunology still largely follows the tradition of analyzing the snapshot of only a small number of specific statistically significant molecules' response. Although analyzing statistically significant genes can help in the explanation of 'local' specific response, to grasp the global regulatory processes requires the comprehensive understanding of genome-wide response [9].
A previous high throughput study on LPS-stimulated murine macrophages (in wildtype, MyD88 KO, TRIF KO, and MyD88/ TRIF double (DKO)) focused on 148 highly expressed genes out of 22690 ORFs based on 3-fold expression increase and 100 expression unit cut-off from 0 to 4 h after LPS stimulation [10]. Although the study showed novel local insights of immunerelated genes in different KOs, it did not show the capacity of LPS to induce pleiotropic biological processes not directly linked to immunity [11,12]. To infer system-level emergent complexity, we re-investigated the same data without any biased expressions cut-off.
Building upon the accuracy and reliability of a correlation metrics based upon thousands of statistical units (genes, ORFs), we investigated temporal whole genome response to LPS stimulation in the above-mentioned 4 different genotypes of macrophages. In contrast to individually analyzing each microarray elements (ORFs), where weakly expressed ORFs usually has been considered to incur high noise-to-signal ratio, we analyzed the temporal expression changes of entire ORFs set considered as a whole. We confirmed LPS induces pleiotropic biological response and, additionally, found two characteristic response modes: acutestochastic and collective. Acute-stochastic mode is largely immunerelated response, while collective mode participates in diverse regulatory processes normally unrelated to immunity that transform the initial local response to the antigenic stimulus into a general state change of the cell. Overall, we show genome-wide differential response mainly occur through lowly expressed genes. This global LPS response has been previously neglected by biased gene-expression cut-offs. Our work can be used to understand both specific and global responses of biological systems.

Results and Discussion
Genome-wide invariance between wildtype, single and double KOs LPS stimulates the MyD88-and TRIF-dependent pathways to activate the innate immune response (Figure 1). We evaluated the correlation structure of Affymetrix mouse expression data for wildtype, MyD88 KO, TRIF KO, DKO at 0, 1 and 4 hours (Methods and [10]). The analyses of only highly up-or downregulated gene expressions ignoring lowly-expressed genes may not be sufficient to unravel genome-wide organizational principles [13][14][15]. To investigate the existence of such principles, we analyzed the whole genome cDNA microarray expression without any threshold cut-off. First, we performed Pearson correlation analysis on the entire genome expression vector. This ensemble property of the population of genes is a robust measure that is not biased by noise at the level of individual gene measurement [16,17]. The whole genome correlations of the same cell-type between different genotypes are extremely high (Pearson r above 0.98, Figure 2A-B), indicating a strong common order parameter influencing the expression level of the entire genome, correspondent to the cell-type characterization [17]. High correlation between genotypes may suggest that technical noise of our microarray dataset is rather low [18,19]. Furthermore, our results are coincident with the results obtained on erythroid cell lineages using the same metrics [20]. The presence of such invariant order spanning more than twenty-thousand elements (genes, ORFs) and around four orders of magnitude of expression levels is a signature of general order parameters organizing the entire cell regulation network. This organization is, in our opinion, a 'fact of nature' that, for its dimensions and invariance, asks for a deep thinking in analogy of what happened for other collective phenomena in physics (magnetism, laser coherence, super-fluid helium, hydrodynamic instabilities). Such strong invariance is imposed by the presence of a common attractor correspondent to the cell kind [16,20], between all genotypes especially when we already know that MyD88 KO and DKO show significantly impaired proinflammatory responses [21]. Hence, to investigate specific proinflammatory and global LPS response, we compared the between genotypes Pearson similarities as computed both on the whole genome and on different extractions of gene subsets (random and immune-related).

Temporal Pearson correlation reveals genotype differences
We adopted two measures to compare different genotypes by means of a Pearson correlation metrics: i) auto-correlation: Pearson r between 0 h (t 0 ) and other time points (t 1 , t 4 ) of the same genotype and ii) cross-correlation: Pearson r between wildtype and other genotypes at same time point (Methods). The auto-correlation analysis measures progressive response from t 0 for the same genotype, while the cross-correlation analysis measures the temporal difference from wildtype response.
The whole genome auto-correlation of all genotypes shows progressive response, correspondent to a progressive displacement of correlation from unity, to LPS stimulation using RMA normalized data ( Figure 3A and Methods). We further checked the consistency of this result using another normalization process, MAS5 data ( Figure S1A). Wildtype response follows distinctly different trend from MyD88 KO, TRIF KO and DKO that in turn are remarkably similar to each other. These results show i) 1% of auto-correlations variability is sufficient to discriminate different genotypes, and ii) the similar global gene expression behavior between single KOs and DKO suggests DKO possesses LPS response, too (see also section 'Emergence of regulatory signature from scattered expressions in all genotypes' for further proof of DKO response).
To find the source for genome-wide response similarity between DKO and single KOs, we calculated the cross-correlation for all genotypes. Cross-correlation shows the response of TRIF KO is the most similar to wildtype, followed by MyD88 KO and DKO ( Figure 3B and Figure S1B). Although the auto-correlations are similar, the different cross-correlations of other genotypes, with wildtype as a reference, show that the gene expression responses between genotypes are distinct from one another. We find these differences originate from differential activation machinery of gene expressions (see section 'Deciphering gene regulatory mechanisms from emergent signature').
Since LPS is a well known inducer of immune response, we next, specifically focused on the auto-correlations of 157 immunerelated genes (according to GenMAPP [22], Table S1). By this, we shift the focus from the 'whole-genome' response to the 'local immune-related' response of the system. The result shows that DKO has a flat profile (as expected, almost perfect correlation between different times is observed due to the lack of any classical immune response to LPS), followed by TRIF KO and MyD88 KO displaying a linear displacement in time from unity correlation but to a lesser extent than wildtype, pointing to a diminished immune response of single KO with respect to wildtype. ( Figure 3C and Figure S1C). These results are consistent with current literature which suggests that MyD88 is key for LPS stimulation and consequently, for DKO, no activation of immune response is expected to occur [23,24].
The cross-correlations of immune-related genes show TRIF KO was closest to wildtype response, followed by MyD88 KO, while DKO was the farthest ( Figure 3D and Figure S1D). Based on these results, we concur that the i) 157 immune-related genes are dependent on both MyD88 and TRIF, and that they are mainly activated by the MyD88-dependent pathway, ii) MyD88 and TRIF cannot be considered to be acting along the same pathway, otherwise single and DKO's cross-correlation should be equal and overlapping and iii) MyD88-and TRIF-dependent pathways are not completely independent, otherwise DKO's distance from unity in terms of cross-correlation would be equal to a composition of MyD88 KO and TRIF KO individual values. These results imply the existence of synergized action between the MyD88-and TRIFdependent pathways because the sum of single KO responses is different from DKO response [25,26] (see MT in 'Deciphering gene regulatory mechanisms from the emergent signature').
In summary, even though we observed genome-wide strong invariance between genotypes, temporal correlation analyses reveal the difference between them. The temporal whole genome (global) response and immune (local-specific) response show distinct profiles [27]. Thus, we pondered whether there are organization principles that distinguish such different modes of response. To uncover, we next investigated Pearson r of all genotypes between 0-1 h by selectively removing ORFs from each genotype. This procedure will allow us to individualize the different contributors to the whole genome behavior.   Biphasic acute-stochastic and collective modes of LPS response We analyzed the change (i.e., difference of expression) of response based on Pearson r of all genotypes by removing highest up-and down-regulated ORFs (one by one up to 300 ORFs) between 0 to 1 hr. For removing highest up-regulated ORFs, a biphasic (hyperbolic) phenomenon emerges in wildtype and single KOs but not in DKO ( Figure 4A). The gradient of curves is steep up to the removal of about 80 ORFs in wildtype and 50 ORFs in both single KOs after which the slope gentles. In contrast, in DKO, only the gentle gradient exists ( Figure 4B). As control, we next compared randomly removed ORFs in similar steps and found Pearson r (0-1 h) of all genotypes does not change noticeably ( Figure 4C). This result points to a transition in the response: a small fraction of 'acute' responding ORFs, and the majority of ORFs responding 'weakly'. Notably, for downregulated ORFs, biphasic response is not observed for any genotype ( Figure 4D). To investigate further, we evaluated the standard deviation of Pearson r (0-1 h) for randomly chosen ORFs in steps of 10 up to 300 (with each selection repeated 30 times) from whole genome and measured the standard deviation (SD) of auto-correlation. We notice that the mean value of SD of auto-correlation decreases with number of ORFs selected ( Figure 4E). Notably, at around 80 ORFs for wildtype (and 50 for single KOs) the auto-correlation transits to a flat trend (DKO, however, did not show any transition). The amplitude of SD of auto-correlation is large for less than 80 ORFs for wildtype (and around 50 for single KOs) and small otherwise ( Figure 4E). These results reveal the coincidence of biphasic transition at around 50 ORFs. Therefore, we can consider acute response as acute-stochastic mode (high SD) and the weak response as collective mode (low SD) [27]. We also found biphasic response for auto-correlations between 1-4 h ( Figure 4F).
Investigating further the temporal response of acute-stochastic and collective modes we observed autoand cross-correlation profiles of acute-stochastic mode is similar to immune-related genes whereas the average response of randomly selected 80 ORFs in collective mode is scaled to genome-wide response ( Figure 4G-H, and Figure 3A), thus demonstrating that collective mode is scalable. This implies the collective mode corresponds to a global, a-specific adjustment of the expression network that can be appreciated even in the case of a sufficiently populated random sample of genes. The collective mode can be considered as a sort of 'mean field' encompassing the entire genome expression and needs a sample of at least 80 ORFs to be reliably detected.

Emergence of regulatory signature from scattered expressions in all genotypes
To understand the temporal progress of biphasic response, we investigated the changes of whole genome expression from early  Figure S2). Remarkably, we observed the transition from large scatter in expression distributions around the origin for single ORFs to smooth lines for group of 50 ORFs onwards for all genotypes. This is due to the fact that grouping of expression distribution for 50 ORFs or above forms Gaussian distribution and the average value of each group remarkably follows linear line ( Figure S3A and see 'Linear regressions analyses' in Methods); the fact average values follow linear lines reveals the emergence of regulatory signature working at the level of groups of genes. Furthermore, the fluctuations on Gaussian distribution reduce as we increase the grouping size ( Figure S3A-C). From these, we observe genes upregulated at early time points were downregulated at later time points and viceversa, through switching in gene regulatory circuits. Biologically, upregulated ORFs in early signaling, when transcription rate is faster than mRNA decay, are downregulated by late signaling, when mRNA decay is prevalent. The vice versa occurs for downregulated ORFs. Notably, this switching behavior is also observed for DKO, thus reinforcing that DKO possess genomewide LPS response through MyD88-and TRIF-independent manner.
Deciphering gene regulatory mechanisms from emergent signature We observed earlier from auto and cross-correlations that the gene expression responses between genotypes are distinct from one another. To understand this in depths, we compared genome-wide expression changes between genotypes for 0-  (Figure 6A-F). We plotted the expression change (Dx) of single ORF as well as taking the average value of each group of ORFs sorted from highest to lowest expressions. We observed i) large scatter in expression distributions around the origin for ORFs in the collective mode and ii) linear expression distribution of ORFs in the acute-stochastic mode ( Figure 6A-F, left panels). These results further confirm biphasic response of LPS stimulation.
For the collective mode, we found the distribution of averages from expression change (Dx) in each group follows transition from scatter to smooth lines for group of 50 ORFs onwards for all genotypes ( Figure 6A-F, right panels & Figure S4), revealing the emergence of genome-wide control in wildtype and all KOs. This is due to the fact that when we group ORFs of 50 and above, the distribution of Dx becomes Gaussian and the mean values of distribution follow linearity in the same way as mentioned in previous section. Another observation is that DKO emergent signature is invariant from single KOs ( Figure 6E-F, right panels). This further suggests that the contribution from DKO is independent from MyD88 and TRIF and derives from an unknown source (U). Thus, the emergent signature as smooth line can be used to decipher quantitative genome-wide differential control of transcriptional (Tr) and mRNA decay machineries Using DAVID analysis platform [28], the majority of biological processes of acute-stochastic mode is related to immune system, defense response, inflammatory response, etc., specifically activated by M and T (Table 1, p,0.05).
Next, we compared groups of upregulated ORFs in the collective mode of wildtype (WT c + ) (x-axis) with MyD88 KO and DKO (yaxes) ( Figure 6A and C, right panels). Unlike acute-stochastic mode, collective mode possesses average points in plot which follows flat profiles. This could be due to the equilibrium of transcriptional and mRNA decay machineries or shift in experimental or normalization process. To investigate this statistically, we selected grouping of equal number of random ORFs and assessed their average expression change in wildtype against all genotypes. Flat profiles were observed with M+T+MT+U = 0 (wildtype), T+U = 20.024 (MyD88 KO), M+U = 0.013 (TRIF KO) and U = 0.024 (DKO) (Figure 7A-C). Since random sampling displays genome-wide constant shift in MyD88 KO, TRIF KO and DKO, these should be due to experimental shift or normalization process. Thus, we eliminated these shifts from MyD88 KO, TRIF KO and DKO contributions.
In general, the emergent smooth lines in collective mode can be presented by y = ax+b ( Figure 8A-C) where slope a represents net outcome between transcriptional and mRNA decay machineries on average expression change and y-axis interception, b, implies contribution from the equilibrium state of transcriptional and mRNA decay machineries. Note that we ignore analysis around the origin of axes, where the densely distributed average points results in overlapping of their Gaussian distributions between genotypes.
For WT c + , we observed same flat profile for DKO in both WT + and randomly selected wildtype ORFs, which indicates that WT + in DKO possess no response, i.e., U = 0 ( Figure 8A    and U = 0.15x ( Figure 8A-C, x,0). This result points to i) MyD88 & TRIF-dependent (MT) pathways mainly show synergized repression of genes through mRNA decay, i.e., Tr%Dm; ii) TRIF-dependent (T) pathways can activate the transcription of the same processes as MT, i.e., Tr.Dm; iii) M has no regulatory role for wildtype downregulated ORFs and iv) unknown processes (U) can weakly repress the same genes through decay machinery, i.e., Dm.Tr.
To determine biological processes regulated by the whole genome, we selected genes satisfying upregulation in TRIF KO and downregulation in MyD88 KO for WT c + and genes satisfying downregulation in TRIF KO and DKO, and upregulation in MyD88 KO for WT c 2 , due to the observed emergent linearity ( Figure 8D). We obtained 3798 out of 11793 ORFs for WT c + and 1916 out of 10897 ORFs for WT c 2 . Next, combining WT c + and WT c 2 ORFs, and using DAVID platform analysis, we retained the biological processes of genes with significant enrichment (Tables 2 and 3, p,0.05).
We observed biological processes related to immunity (Cytokinecytokine receptor interaction, Toll-like receptor signaling pathway, etc.) were upregulated in WT c + , which indicates immune-related genes are not restricted to the acute-stochastic mode alone (Table 2). Predominantly, however, in WT c + , a) pre-transcriptional and transcription-related genes, such as genes related to cell surface receptor mediated signal transduction (including Protein phosphorylation required for signal transduction and Metabolism of cyclic nucleotides used for kinase activities), and mRNA transcription regulation and b) genes related to post-translational processes (Proteolysis, Exocytosis, etc.) were upregulated (Table 2). In WT c 2 , c) genes related to posttranscriptional processes, (Pre-mRNA processing, mRNA splicing, Protein biosynthesis, Ribosome, etc.) were downregulated (Table 3). We therefore hypothesize that wildtype cells prepare for secondary immune activation (possibly through cytokine receptors) by upregulating signaling and transcription processes ( Figure 8E). To self-regulate the secondary immune response, post-transcriptional processes (mRNA splicing and translation) are repressed. Interestingly, we observe, from Figure 8E, a complete reverse (switching) behavior in MyD88 KO where only T and U are active, e.g. pre-transcriptional processes and transcription become downregulated while post transcriptional processes are upregulated. This suggests MyD88 KO cells compensate lack of activation (transcription) by enhancing post-transcriptional processes through TRIF-dependent (T) pathways [29]. These findings seem to indicate that MyD88 is a key regulator in collective mode.
In summary, we observe differential activation of group of ORFs between genotypes (Table S2). From these results, we obtained the individual effects of M, T, MT and U on a genomewide scale ( Figure 8E). This shows genome-wide differential activation machinery of biological processes. Unlike acute-stochastic mode which is dedicated to immune system, collective mode may also participate in diverse regulatory processes; to prepare the cell for activation (signal transduction), prevent over expression of genes (regulation of post-transcriptional processes), and compensate lack of activation in KO conditions (enhancement of genes related to post-transcription when transcription process is lacking).

Conclusion
In this report, we found quantitative genome-wide differential control of transcriptional and mRNA decay machineries through signaling processes superimposing over the general strong coregulation of expression levels that is largely invariant between genotypes and related to the global expression attractor correspondent to the specific cell type [15,20]. Moreover, each genotype, except DKO, possess two modes of responses; acutestochastic (small number of immune-related genes) and collective modes (rest of ORFs). The collective mode, which consists of myriad cellular processes, is often ignored in most analyses as they are made of ORFs displaying small expression changes in time and hence cannot be captured if high cut-off thresholds (e.g. 3-fold) are used. Also in collective mode, for all genotypes, we notice scalable response and emergent linear behavior arise when ORFs are  grouped. Other manifestations of such collective behavior, arising from the functional relations between gene expressions, were observed in terms of coordinated genome-wide expression waves [30,31]. The transition from scatter to linearity was observed in the distribution of whole genome expression when grouping of 50 ORFs and above. Notably similar transition occurs for the distribution of single gene expression of cells in culture when the intrinsic (uncorrelated) noise becomes low [32]. Thus, our work also reveals the regulation by correlations in gene expression fluctuations [33]. However, it is important to stress that our data refers to large ensembles of cells, unlike single cell measurements, and thus exact discrimination between intrinsic and extrinsic noise sources cannot be performed in a similar manner [32,33]. Nevertheless, the strong invariance between different conditions of the same cell-type, can be considered as a sort of dynamical attractor encompassing the entire transcriptome, reflecting hidden genome-wide differential regulations [27]. Understanding the link between the ordered behaviors observed for i) single gene expression when intrinsic noise is low [33] and ii) genome-wide conditions, promises to be a very fruitful future direction.
The discovery of two modes of response has also been shown recently for protein dynamics to a drug perturbation where a rapid translocation of specific proteins and a slower, wide-ranging temporal wave of protein degradation and accumulation occurred [34]. Our work points to the presence of a highly-ordered, coordinated, genome-wide expression dynamics of LPS stimulation, thereby requesting the need to consider global phenomena when interpreting immune response. In general, the consideration of the general rearrangement of the entire expression network after a specific stimulus, with the consequent activation of functions not directly linked to the original perturbation, could be the basis for rationalizing the onset of unexpected side-effects after drug treatments.

Biological datasets
We re-analyzed microarray data obtained from time-series experiments (0, 1, and 4 hours) performed on peritoneal macrophages from wildtype, MyD88 KO, TRIF KO, and MyD88/TRIF DKO mice treated with 100 ng/ml of LPS (Salmonella Minnesota Re595, Sigma) [10]. Affymetrix mouse expression array A430 microarray chips were used for gene expression detection. The microarray dataset obtained from these experiments contains expression levels for 22690 Affymetrix probe set IDs. We reprocessed our Affymetrix microarray chip data using Robust Multichip Average (RMA) for further background adjustment and to reduce false positives of our Affymetrix microarray chip [35][36][37]. The complete experimental details can be found in Hirotani et al [10].

Statistical Analysis
Auto-and cross-correlation analysis for interpreting LPS response. To investigate the correlation between any two expression vectors, x and y with n ( = 22690) dimension and mean values of expressions x and y, we calculate their mutual Pearson, r~x,y ð Þ by where X~x 1 {x,x 2 {x,::,x n {x ð Þ , Y~y 1 {y,y 2 {y,::,y n {y ð Þ and h is angle between two expression vectors. Geometrically, Eq. 1 shows the correlation coefficient, r(x,y), can be viewed as the cosine of the angle (h) on n-dimensional space between the two vectors of data representing a measure of response. However, when h = 0 (i.e., r = 1), generally X = aY (a.0). In the case when a = 1 i.e., X = Y, this implies X and Y has the same response, otherwise, X and Y have different but globally proportional response.
We extend the Pearson correlation analysis to measure global temporal gene expression response to a given stimulation,  comparing Pearson correlation between 0 h (X t 0 ð Þ) and all time points (Y~X t 0 ð Þ,X t 1 ð Þ,X t 4 ð Þ) at 0 h, 1 h, 4 h of the same sample vectors, auto-correlation. Therefore, the auto-correlation profiles measures progressive divergence of expression from t 0 for each genotype in terms of decreasing correlation in time, if response of LPS stimulation occurs. Since Pearson correlations of whole genome for all condition are close to one (i.e., h>0), we need to distinguish whether a = 1 or not. We add 0 h vector elements of X into both X and Y, resulting in 2n-dimension; X = (X(t 0 ), X(t 0 )) and Y = (X(t), X(t 0 )). If X(t) = X(t 0 ), that is, no response, then X : , thus, auto-correlation = 1 (h = 0). On the other hand, if X(t) = a X(t 0 ) (a?1 & a.0), auto-correlation, r?1 (i.e. h?0). Biologically, autocorrelation with control profiles will show progressive divergence from 0 h expression for each genotype if dynamic response to LPS exists,.
Similarly to auto-correlation, cross-correlation is a temporal Pearson correlation measure. However, instead of measuring between the same genotypes with different time points, cross-correlation measures between different genotypes from wildtype at the same time points.

Linear regressions analyses
To obtain reliable linear regressions in the analysis of expression change distributions ( Figure 8A-C) for the collective mode, we determined the minimal number of ORFs (N) per group, to observe the formation of a Gaussian distribution centering on the average expression change of the group of ORFs in all genotypes. Shapiro-Wilk test was performed to assess the normality of the distribution of expression changes in each group of N (0,N,2000 by step of 10) ORFs, which provides a statistic W coupled with pvalue of the data (Higher W and lower p-values for normally distributed data). The average of p-values computed were lower than 0.1 (2-tailed) for groups of about N.300 for MyD88 KO distribution, N.800 for TRIF KO and N.1000 for DKO ( Figure  S5). We therefore used N = 1000 to determine the linear regressions of expression changes distributions in collective mode.

Functional enrichment analyses
DAVID functional annotation platform [28] was used to identify the functional categories (Gene Ontology (GO) [38], Panther gene classification [39] or KEGG pathways [40]) enriched in groups of ORFs. Among the 22690 ORFs, 10264 genes with annotations in Gene Ontology, 13824 genes with annotations in Panther gene classification, and 3778 genes with annotations in KEGG pathways were identified. To evaluate functional category enrichment under control of False Discovery Rate, Benjamini-Hochberg adjusted p-values were obtained for each term (GO, Panther, KEGG), and terms scoring p-values,0.05 were retained.

Supporting Information
Table S1 List of immune-related genes. List of 157 immunerelated genes selected from GenMAPP used for analysis.