Deregulation in lipid metabolism leads to the onset of hepatic steatosis while at subsequent stages of disease development, the induction of inflammation, marks the transition of steatosis to non-alcoholic steatohepatitis. While differential gene expression unveils individual genes that are deregulated at different stages of disease development, how the whole transcriptome is deregulated in steatosis remains unclear.
Using outbred deer mice fed with high fat as a model, we assessed the correlation of each transcript with every other transcript in the transcriptome. The onset of steatosis in the liver was also evaluated histologically.
Our results indicate that transcriptional reprogramming directing immune cell engagement proceeds robustly, even in the absence of histologically detectable steatosis, following administration of high fat diet. In the liver transcriptomes of animals with steatosis, a preference for the engagement of regulators of T cell activation and myeloid leukocyte differentiation was also recorded as opposed to the steatosis-free livers at which non-specific lymphocytic activation was seen. As compared to controls, in the animals with steatosis, transcriptome was subjected to more widespread reorganization while in the animals without steatosis, reorganization was less extensive. Comparison of the steatosis and non-steatosis livers showed high retention of coordination suggesting that diet supersedes pathology in shaping the transcriptome’s profile.
This highly versatile strategy suggests that the molecular changes inducing inflammation proceed robustly even before any evidence of steatohepatitis is recorded, either histologically or by differential expression analysis.
Non-alcoholic steatohepatitis (NASH) develops in livers that have accumulated histopathological changes associated with hepatic steatosis and are reflected to the differential expression of genes linked to the induction of inflammation [1,2,3]. During disease progression extensive transcriptional reprogramming occurs that underscores its different stages. This multistage process can be recapitulated with relatively high accuracy in animal models receiving special diets, alone or combined with other stimuli triggering liver injury . Among them, outbred models may be of special value since they can mimic the different courses of disease progression in human patients at which steatosis develops stochastically .
Essential for the molecular characterization of different subtypes of liver disease is differential expression which points to specific transcripts that are enriched or depleted at different disease stages [6,7,8]. Such quantitative changes in expression, usually illuminate full-fledged pathology while subtle alterations, despite their potential significance may remain unnoticeable. To overcome these limitations, we have proposed a novel strategy that instead of the expression levels of individual transcripts takes into consideration the degree of correlation of expression of each transcript with every other transcript in the transcriptome. Evaluation of this coordination profile of the whole liver transcriptome at different disease stages, may provide hints regarding the underlying molecular changes that conventional, differential expression analysis cannot. Such changes in coordination profile appeared relevant in characterizing different liver pathology stages by focusing on the unfolded protein response [9,10,11]. Furthermore, such analysis applied to the most highly expressed genes in the transcriptome was shown capable of illustrating changes in patients with frailty syndrome . More recently, application of this strategy to the transcriptome of deer mice (Peromyscus) indicated changes in gene expression profiles and associated biological processes that occur in the brain during aging in different species .
In the present study we used outbred deer mice (Peromyscus maniculatus) as a model to evaluate how the liver transcriptome is collectively reorganized in specimens with or without steatosis. Our studies were based on our earlier findings indicating that upon high fat diet administration (HFD) P. maniculatus do not develop steatohepatitis, but a subset, about 50%, develops hepatic steatosis . This provides an animal model at which disease development can recapitulate human genetically diverse populations at which disease susceptibility varies considerably among individuals .
The premise of the present analysis is that genes belonging to the same transcriptional networks are co-expressed, but when pathology emerges the profile of co-expression is collectively changed [15,16,17,18,19,20,21]. To address the coordination profile, we calculated the composite correlation for each gene in the transcriptome with every other gene and compared it in the controls, the animals that received HFD but did not develop steatosis and the animals that received HFD but developed steatosis. The results were coupled to gene ontology (GO) analyses  to reveal transcripts that more prominently abolished their coordination with the whole transcriptome. This analysis was recently applied to the brain transcriptome of deer mice and showed that during aging a loss of the perception of smell occurs in different species despite that de-coordination of the transcriptome exhibit interspecific variations . Our present results, besides describing the overall coordination profile of the transcriptome at different conditions, also show that HFD triggers a robust induction of an inflammatory response, irrespectively of the onset of steatosis. This change was not apparent by conventional analyses of the differentially expressed transcripts. Furthermore, it showed that what differentiates the liver transcriptomes with and without steatosis is the preference of the former for T cell activation, myeloid leukocyte differentiation and engagement of genes involved in cell cycle regulation.
Variable response to HFD in outbred deer mice
Earlier studies showed that P. maniculatus exhibits variable response to HFD. In the present study, a panel of 3–4 months old, outbred deer mice (P. maniculatus) received HFD for about 6 months (n = 10). Six animals received regular diet. Body weight was increased in the animals that received the HFD but body weight gain remained highly variable, consistently with the genetically diverse nature of the experimental animals (Fig. 1a). Histology revealed the presence of steatosis in 5 out of 10 animals that received HFD (Fig. 1b). No evidence of fibrosis, ballooning degeneration or lobular inflammation was recorded [23,24,25], suggesting that under these conditions the disease has not progressed to more advanced stages of non-alcoholic liver steatohepatitis (NASH) (Fig. 1b).
Distinct profile of expression coordination in livers with or without steatosis
RNAseq was performed in the liver of P. maniculatus that received regular diet or HFD and results have been deposited to NCBI (GSE146846). To test how the transcriptome in each group is coordinated at these conditions we calculated the composite correlation (Pearson’s composite, Pc) index as follows: Initially we calculated the correlation coefficient for each gene with every other expressed gene in the transcriptome, in all 3 pairwise comparisons being control vs steatosis (C vs S), control vs non-steatosis (C vs NS) and steatosis vs non steatosis (S vs NS) (Supplementary Tables 1, 2, 3). The analysis was performed by using a code in R language we have developed and involved 13,340 transcripts in the NS vs C comparisons, 13,434 transcripts in the S vs C comparisons and 12,041 transcripts in the S vs NS comparisons. In all cases, the expression of each of these transcripts was evaluated versus all other transcripts in the transcriptome of the experimental specimens. Pc of each transcript reflects the composite correlation coefficient of all correlation coefficients that were calculated as described above, for each given transcript, in the 3 pairwise comparisons (C vs S, C vs NS, and S vs NS). Therefore, high Pc values indicate that coordination is retained for the given comparison, while lower Pc values indicate loss of coordination. Conversely, negative Pc values show that the profile of coordination is inversed in a manner according to which positive correlation with the transcriptome in one experimental condition is inversed to negative in the other, and vice versa.
As shown in Fig. 2a, in all three groups, most of the transcripts showed positive Pc values, which indicates that coordination is retained between the animals with or without steatosis, for most of the genes. Higher Pc values (average Pc = 0.17) were seen in the S vs NS groups suggesting that most genes retained their coordination upon HFD administration and irrespectively of the development of steatosis (Fig. 1b,c). Conversely, lowest Pc (=0.086) was seen in the comparison between S and C suggesting extensive transcriptional reprogramming. In the comparison between NS and C, average Pc had intermediate magnitude (Pc = 0.13). All differences were statistically significant (ANOVA, P < 0.0001). Similar were the findings when instead of the whole transcriptome only genes common in the 3 groups were evaluated suggesting that the findings do not reflect a bias towards transcripts that are present only in some experimental groups (Fig. 2).
These results suggest that during HFD administration, extensive reprogramming of the transcriptome occurs, which is more pronounced in the livers that develop steatosis as compared to those that did not. The differences in the transcriptomic profile between the livers that did and those that did not develop steatosis at HFD, were more modest (Pc = 0.17). Thus, special diet such as HFD induces more changes in the transcriptome than the pathology (steatosis) per se.
Gene ontology analyses reveal engagement of inflammation by HFD
To obtain insights regarding the biological processes that are enriched for the transcripts exhibiting the most pronounced changes in Pc values in the 3 comparison groups we utilized the publicly available Gene Ontology Platform (Gene Ontology http://geneontology.org/). For this analysis, the Pc values were sorted for each group in descending order and the genes exhibiting Pc < − 0.2 were analyzed (Suppl Table 4). The results for the top 10 processes are shown in Table 1 and were derived by using 752 genes, 700 genes and 854 genes for the S vs C, the NS vs C and the S vs NS genes respectively. Both comparisons involving administration of HFD (S and NS) vs C exhibited an enrichment for processes associated with a proinflammatory response. Thus, robust transcriptional reprogramming, consistent with the induction of inflammation, occurs irrespectively of steatosis and despite that no histopathological evidence of inflammation was seen. Comparison between the S vs NS specimens revealed that the most prominent processes were associated with cell cycle regulation.
In order to identify genes that collectively exhibit changes in their transcriptomic profile across the different groups, we calculated a cumulative Pc index by adding the 3 independent Pc indices for the genes that were common between the 3 individual pairwise comparisons. Then we sorted the genes in descending order according to their cumulative Pc and selected the top 5% which corresponded to about 600 genes for GO analysis (Table 1 and Suppl Table 5). Not surprisingly, GO analysis indicated that processes associated with metabolism were more prominently enriched.
Differential expression only reveals a fraction of processes linked to steatosis
To appreciate the discovery power of the proposed transcriptome coordination strategy, we also performed conventional differential expression and GO analysis (Suppl Table 6). For this analysis the iDEP online platform was used by using FDR cut-off 0.1 and minimum fold change 2. As shown in Table 2 the results were not highly consistent and informative in the different groups and suggested enrichment of processes associated with lipid metabolism in the NS vs C, while in the S vs C comparison, processes associated with cytokinesis, mitotic spindle formation and epithelial to mesenchymal transition. In the S vs NS comparison, the most prominent processes enriched were related to apoptosis and regeneration. The number of differentially expressed genes and the top 3 upregulated and downregulated transcripts in each comparison are shown in Fig. 3.
Immune cell markers do not exhibit differential expression in the livers of deer mice receiving HFD
Despite that immune cell activation was not detectable by the conventional differential expression analysis it remains plausible that individual mediators of the proinflammatory response are differentially expressed. We tested this by evaluating the expression of a panel of immune cell markers including IL1b, IL15, IL18, CCL2, CCL6, CCL20, CXCL1, CXCL9, CXCL10, CXCL12, CXCL16, and CX3CR1 but none showed differential expression levels between the experimental groups (Fig. 4, not significant, Ordinary one-way ANOVA, Tukey’s multiple comparisons test).
The assessment of differential expression is an important indicator of genes associated with pathology however its value can be limited when genetically diverse specimens are analyzed. Genes highly relevant to disease may remain masked if the variation in gene expression between individuals reduces the statistical power of differential expression studies. For example, in hepatic steatosis, despite the established role of endoplasmic reticulum stress in disease development, genes associated with the unfolded protein response are not usually detected by differential expression analysis . Such role though can be revealed when their coordination with the whole transcriptome is examined . When the role of inflammation is studied in liver disease, it marks only its more advanced stages and is frequently dissociated from steatosis, especially in some animal models [26,27,28,29]. While the deregulation of pro-inflammatory cytokines is detected in benign steatosis in the absence of typical liver inflammation they are generally considered as the direct outcome of aberrant lipid metabolism, occasionally originating from visceral fat, are linked to insulin resistance and are not representative of an orchestrated inflammatory response occurring in the liver [30,31,32].
By using a novel, unbiased whole transcriptome analysis, that relies on the extent of expression of all transcripts in livers from outbred rodents that developed or not steatosis after HFD administration, we were able to show that both in the specimens that did not show pathology and those that exhibited steatosis, a robust engagement of proinflammatory processes occurred. What however differentiated the two entities was the engagement of T cell activation and myeloid leukocyte differentiation processes that was detected only in the fatty livers and could be due to the presence of genetic modifiers that differentiates the responses recorded between the groups. Conventional differential expression analysis that focuses on transcripts exhibiting quantitative differences in the experimental groups, failed to reveal evidence of immune cell activation. This was further supported by the lack of significant changes in the expression of a panel of established mediators of inflammation in liver disease. Probably this limitation is related to the genetically diverse nature of the specimens in combination with the fact that such changes may be below the thresholds of significance of such analysis. Furthermore, it may indicate that in this animal model, under the present conditions and at this stage of disease development, only some transcriptional reprogramming had occurred that has primed cells for the mobilization of the immune cells. The latter, however, has not occurred substantially, as yet, because differential expression was not implemented. Such priming for transcriptional reprogramming may have committed the liver tissue into a permissive state at which pro-inflammatory changes occurred, even prior to the engagement of specific markers that are diagnostic for immune cell engagement. Coordination analysis, especially at the whole transcriptome level, leverages such diversity in gene expression among individual specimens and can extract meaningful information even when subtle changes occurred. To that end, it is conceivable that prolonged high fat diet administration and/or diet richer in fat than the one presently used, may eventually cause histologically detectable changes consistent with hepatic inflammation and steatosis in the animals that received HFD but did not develop steatosis so far. Consistently with this notion, it is plausible that application of this analysis in human specimens may be capable of revealing changes prior to the manifestation of specific histological or molecular alterations.
This coordination analysis also indicated in pairwise comparisons that a major difference of the livers with and without steatosis, as compared to the controls is that in those with steatosis, the transcriptome underwent more extensive reorganization compared to those without. Comparison though of the two, exhibited the higher retention in the profile of coordination. This suggests that diet supersedes pathology in shaping the profile of the transcriptome. From a different perspective this observation implies that for the emergence of pathology, such as steatosis, a more restricted roster of changes in gene expression is sufficient, while more global changes, likely imply the operation of homeostatic mechanisms that maintain normal cellular function.
Conventional differential expression analysis of gene expression was only able to reveal inconsistent changes between the experimental groups, among which the most prominent were related to lipid metabolism in NS vs C, cytokinesis, mitotic spindle formation and epithelial to mesenchymal transition in S vs C, and processes related to apoptosis and regeneration in S vs NS comparison.
Collectively, these results suggest that inflammatory engagement is robustly triggered by HFD even before inflammation is detectable in the histopathological analysis or by the differential expression studies, and illustrate the power of the proposed gene coordination approach to reveal changes that conventional strategies cannot.
Materials and methods
Deer mice, P. maniculatus were obtained from the Peromyscus Genetic Stock Center (PGSC) of the University of South Carolina (UofSC), Columbia, SC (RRID:SCR_002769). Deer mice, P. maniculatus bairdii (BW Stock), were bred as a closed colony in captivity, since 1948 and descended from 40 ancestors wild-caught near Ann Arbor, Michigan . Deer mice were fed either a regular chow diet (n = 6) or a high fat diet (HFD, 58 kcal% fat and sucrose, Research Diets D12331) (n = 10) for 28 weeks, starting at 3–4 months of age. A larger number of animals received HFD to obtain cohorts with and without steatosis . Body weight was measured every 2 weeks. Animals were then sacrificed using isoflurane as an anesthetic followed by cervical dislocation, and the livers were collected. All animal procedures were approved by the Institutional Animal Care and Use Committee (IACUC) and the Department of Health and Human Services, Office of Laboratory Animal Welfare, University of South Carolina (Approval No. 2349–101,211-041917). All methods were carried out in accordance with relevant guidelines and regulations. The study was carried out in compliance with the ARRIVE guidelines.
RNA and library preparation, sequencing, and postprocessing of the raw data and data analysis were performed by the USC CTT COBRE Functional Genomics Core. RNAs were extracted with a Qiagen RNeasy Plus Mini kit as per manufacturer’s recommendations (Qiagen, Valencia, CA). RNA integrity was assessed using the Agilent Bioanalyzer and samples had a quality score ≥ 8.6. RNA libraries were prepared using established protocol with NEBNext Ultra II Directional RNA Library Prep Kit for Illumina (NEB, Lynn, MA). Each library was made with one of the TruSeq barcode index sequences and pooled together into one sample to be sequenced on the HiSeq 2x150bp pair-ended sequencing platform (Genewiz). Sequences were aligned to the P. maniculatus genome (HU_Pman_2.1 (GCA_003704035.1)) in ensembl.org using STAR v2.7.2 . Reads were counted using the featureCounts function of the Subreads package  using Ensembl GTF and summarized at exon, transcript, or gene level. Only reads that were mapped uniquely to the genome were used. The differentially expressed gene analysis was conducted with iDEP.91 (iDEP Platform http://bioinformatics.sdstate.edu/idep/) . The criteria used for this analysis were FDR with cut-off 0.1 and minimum fold change in gene expression of 2.
Upon termination of the study the liver of the animals was removed, fixed in 10% neutral buffered formalin and paraffin embedded. The livers were stained with H&E and were histologically evaluated. Histological examination of the liver specimens was performed blindly for the presence of hepatic steatosis according to the scoring system designed by the Pathology Committee of the NASH Clinical Research Network, which addresses the full spectrum of lesions of NAFLD . Images shown were obtained by a Leica ICC50 HD.
The correlation coefficients (R, Person’s) values of each gene against all other genes in the transcriptome were calculated in specimens of steatosis vs nonsteatosis, steatosis vs control and nonsteatosis vs control, respectively. The composite correlation (Pc) index was evaluated by calculating the Person’s R of all R coefficients, of each gene in each group combination. This transformation resulted in the development of a unique Pc value corresponding to each of these genes. This reflects the extent by which coordination with the whole transcriptome changes in the corresponding gene. All calculations were conducted with R 3.6.3.
For differential expression results were analyzed by using the iDEP platform default values (FDR cut-off = 0.1 and minimum fold change = 2). For GO analyses Fisher test was used by applying FDR correction (details in corresponding supplementary Tables). For correlation studies, R value from Pearson’s correlation was calculated. Comparisons between groups for both Pc values and differential expression of qPCR results were performed by ordinary one-way ANOVA and by applying Tukey’s multiple comparisons test. In all case results were considered significant when P < 0.05.
Availability of data and materials
Deer mice are available from the Peromyscus Genetic Stock Center. Data were deposited to NCBI. Accession No. GSE146846.
Farrell G, Schattenberg JM, Leclercq I, Yeh MM, Goldin R, Teoh N, et al. Mouse models of nonalcoholic steatohepatitis: toward optimization of their relevance to human nonalcoholic steatohepatitis. Hepatology. 2019;69(5):2241–57. https://doi.org/10.1002/hep.30333. PMID: 30372785.
Havighorst A, Zhang Y, Farmaki E, Kaza V, Chatzistamou I, Kiaris H. Differential regulation of the unfolded protein response in outbred deer mice and susceptibility to metabolic disease. Dis Model Mech. 2019;12(2). https://doi.org/10.1242/dmm.037242.
Hoang SA, Oseini A, Feaver RE, Cole BK, Asgharpour A, Vincent R, Siddiqui M, Lawson MJ, Day NC, Taylor JM, Wamhoff BR, Mirshahi F, Contos MJ, Idowu M, Sanyal AJ. Gene expression predicts histological severity and reveals distinct molecular profiles of nonalcoholic fatty liver disease. Sci Rep. 2019;9:12541. https://doi.org/10.1038/s41598-019-48746-5, 1.
Bertola A, Bonnafous S, Anty R, et al. Hepatic expression patterns of inflammatory and immune response genes associated with obesity and NASH in morbidly obese patients. PLoS One. 2010;5(10):e13577. Published 2010 Oct 22. https://doi.org/10.1371/journal.pone.0013577.
Morrison MC, Kleemann R, van Koppen A, Hanemaaijer R, Verschuren L. Key inflammatory processes in human NASH are reflected in ldlr−/−.leiden mice: a translational gene profiling study. Front Physiol. 2018;9:132. Published 2018 Feb 23. https://doi.org/10.3389/fphys.2018.00132.
Zhang Y, Chatzistamou I, Kiaris H. Coordination of the unfolded protein response during hepatic steatosis identifies CHOP as a specific regulator of hepatocyte ballooning. Cell Stress Chaperones. 2020;25(6):969–78. https://doi.org/10.1007/s12192-020-01132-x.
Soltanmohammadi E, Farmaki E, Zhang Y, Naderi A, Kaza V, Chatzistamou I, et al. Coordination in the unfolded protein response during aging in outbred deer mice. Exp Gerontol. 2020;144:111191. https://doi.org/10.1016/j.exger.2020.111191.
Zhang Y, Lucius MD, Altomare D, Havighorst A, Farmaki E, Chatzistamou I, et al. Coordination analysis of gene expression points to the relative impact of different regulators during endoplasmic reticulum stress. DNA Cell Biol. 2019;38(9):969–81. https://doi.org/10.1089/dna.2019.491.
Zhang Y, Chatzistamou I, Kiaris H. Identification of frailty-associated genes by coordination analysis of gene expression. Aging (Albany NY). 2020. https://doi.org/10.18632/aging.102875
Soltanmohammadi E, Zhang Y, Chatzistamou I, Kiaris H. Resilience, plasticity and robustness in gene expression during aging in the brain of outbred deer mice. BMC Genomics. 2021;22(1):291. https://doi.org/10.1186/s12864-021-07613-2.
Roy S, Bhattacharyya DK, Kalita JK. Reconstruction of gene co-expression network from microarray data using local expression patterns. BMC Bioinformatics. 2014;15(Suppl 7):S10. https://doi.org/10.1186/1471-2105-15-S7-S10.
Luo J, Xu P, Cao P, Wan H, Lv X, Xu S, et al. Integrating genetic and gene co-expression analysis identifies gene networks involved in alcohol and stress responses. Front Mol Neurosci. 2018;11:102. https://doi.org/10.3389/fnmol.2018.00102. PMID: 29674951; PMCID: PMC5895640.
van Dam S, Võsa U, van der Graaf A, Franke L, de Magalhães JP. Gene co-expression analysis for functional classification and gene-disease predictions. Brief Bioinform. 2018;19(4):575–92. https://doi.org/10.1093/bib/bbw139. PMID: 28077403; PMCID: PMC6054162.
Amar D, Safer H, Shamir R. Dissection of regulatory networks that are altered in disease via differential co-expression. PLoS Comput Biol. 2013;9(3):e1002955. https://doi.org/10.1371/journal.pcbi.1002955. Epub 2013 Mar 7. PMID: 23505361; PMCID: PMC3591264.
Ashburner M, Ball CA, Blake JA, Botstein D, Butler H, Cherry JM, Davis AP, Dolinski K, Dwight SS, Eppig JT, Harris MA, Hill DP, Issel-Tarver L, Kasarskis A, Lewis S, Matese JC, Richardson JE, Ringwald M, Rubin GM, Sherlock G, And the gene ontology consortium. Gene ontology: tool for the unification of biology. Nat Genet. 2000;25:25–29. https://doi.org/10.1038/75556.
Lackner C, Gogg-Kamerer M, Zatloukal K, Stumptner C, Brunt EM, Denk H. Ballooned hepatocytes in steatohepatitis: the value of keratin immunohistochemistry for diagnosis. J Hepatol. 2008;48(5):821–8. https://doi.org/10.1016/j.jhep.2008.01.026. Epub 2008 Feb 22. PMID: 18329127.
Rodriguez-Suarez E, Mato JM, Elortza F. Proteomics analysis of human nonalcoholic fatty liver. In: Josic D, Hixson D, editors. Liver proteomics. Methods in molecular biology (methods and protocols), vol. 909. Totowa: Humana press; 2012.
Wang W, Xu M-J, Cai Y, Zhou Z, Cao H, Mukhopadhyay P, Pacher P, Zheng S, Gonzalez FJ, Gao B. Inflammation is independent of steatosis in a murine model of steatohepatitis. Hepatology. 2017;66:108–123. https://doi.org/10.1002/hep.29129.
Bradbury MW. Lipid metabolism and liver inflammation. I Hepatic fatty acid uptake: possible role in steatosis. Am J Physiol Gastrointest Liver Physiol. 2006;290:G194–8.
Targher G, Bertolini L, Scala L, Zoppini G, Zenari L, Falezza G. (2005), Non-alcoholic hepatic steatosis and its relation to increased plasma biomarkers of inflammation and endothelial dysfunction in non-diabetic men. Role of visceral adipose tissue. Diabet Med. 2005;22:1354–1358. https://doi.org/10.1111/j.1464-5491.2005.01646.x.
Kleiner DE, Brunt EM, Van Natta M, Behling C, Contos MJ, Cummings OW, et al. Design and validation of a histological scoring system for nonalcoholic fatty liver disease. Hepatology. 2005;41(6):1313–21. https://doi.org/10.1002/hep.20701.
YZ performed experiments, designed study, developed bioinformatic tools, analyzed data and edited manuscript; IC performed histology, analyzed data and edited manuscript. HK designed study, analyzed data and wrote manuscript draft. The author(s) read and approved the final manuscript.
The study was approved by the Institutional Animal Care and Use Committee (IACUC) and the Department of Health and Human Services, Office of Laboratory Animal Welfare, University of South Carolina (Approval No. 2349–101211-041917). All methods were carried out in accordance with relevant guidelines and regulations. The study was carried out in compliance with the ARRIVE guidelines.
Consent for publication
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
GO analysis for differentially expressed genes in all experimental groups comparisons.
Rights and permissions
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.