Zhuoran Gao
Ran Ding2
Xiangyun Zhai
Yuhao Wang
Zhi-Qiang Du- 1College of Animal Science, Yangtze University, Jingzhou, China
- 2College of Animal Science and Technology, Northeast Agricultural University, Harbin, China
Excessive fat deposition can cause chicken health problem, and affect production efficiency by causing great economic losses to the industry. However, the molecular underpinnings of the complex adiposity trait remain elusive. In the current study, we constructed and compared the gene co-expression networks on four transcriptome profiling datasets, from two chicken lines under divergent selection for abdominal fat contents, in an attempt to dissect network compositions underlying adipose tissue growth and development. After functional enrichment analysis, nine network modules important to adipogenesis were discovered to be involved in lipid metabolism, PPAR and insulin signaling pathways, and contained hub genes related to adipogenesis, cell cycle, inflammation, and protein synthesis. Moreover, after additional functional annotation and network module comparisons, common sub-modules of similar functionality for chicken fat deposition were identified for different chicken lines, apart from modules specific to each chicken line. We further validated the lysosome pathway, and found TFEB and its downstream target genes showed similar expression patterns along with chicken preadipocyte differentiation. Our findings could provide novel insights into the genetic basis of complex adiposity traits, as well as human obesity and related metabolic diseases.
Introduction
The global obesity pandemic and related metabolic syndromes are currently devastating the human society, by threatening human health and decreasing life expectancy (Blüher, 2019). To find effective therapy for metabolic diseases, current research efforts employ large-scale genomics and systems biology methods, to understand better the biology and physiology of adipose tissues (Dyar et al., 2018; Guan et al., 2018; Pellegrinelli et al., 2018; Zhao et al., 2018; Brial et al., 2019; Justice et al., 2019; van der Klaauw et al., 2019). However, adipose tissues are of complex and heterogeneous origins, such as depot-specific (subcutaneous, abdominal, inguinal, etc.), different tissue-types (white, brown, and beige), and composition of different cell types (immune cells, preadipocytes, stromal cells, neurons, etc.) (Lee et al., 2014; Jeffery et al., 2016; Ghaben and Scherer, 2019; Luong et al., 2019; Sebo and Rodeheffer, 2019). It’s rather difficult to disentangle the molecular circuits driving the growth and development of adipose tissue (Tang and Lane, 2012; Ahrends et al., 2014; Cohen and Spiegelman, 2016).
Gene network approach, as one method of integrative analysis on high-throughput transcriptome profiling and a variety of other omics data, helps discover successfully the structural and functional gene modules and molecular signaling pathways for human diseases (Barabási and Oltvai, 2004). Different statistical and machine learning methods were developed and refined for the construction of gene networks, e.g., the gene coexpression network (Langfelder and Horvath, 2008; Ranola et al., 2013), and the Bayesian network (Chen et al., 2008; Yang et al., 2009). Later, integrative genomics methods are effective in combining different omics data (e.g., transcripomics, proteomics, metabolomics), to pinpoint key biochemical and molecular biomarkers (Lee et al., 2016; Emilsson et al., 2018; Ghaemi et al., 2019). In farm animals, the gene-coexpression network were recently used to study cattle reproduction, feed intake, meat quality, immune response, and functional annotation of gene functions (Fortes et al., 2010; Beiki et al., 2016, 2018; Buchanan et al., 2016; Mateescu et al., 2017; de Oliveira et al., 2018; Gonçalves et al., 2018; Alexander et al., 2019). In chickens, network methods were also employed on the investigation of growth and reproduction traits (Zhang et al., 2017; Wang et al., 2018).
Intensive selection on growth rate and feed efficiency traits in the past several decades makes the broiler industry one of the most efficient animal production systems. However, fast growth could bring along excessive fat deposition, and cause economic loss and processing burden to the broiler industry. To breed broiler lines with less fat, recent efforts are focused on understanding the molecular genetics of adipose tissue growth and development in the chicken (Siegel, 2014). Broiler lines were divergently selected for abdominal fat content (Guo et al., 2011; Baéza and Le Bihan-Duval, 2013), and a systematic approach integrating genetic, genomic, cellular, and molecular studies were performed, to discover important genes and molecular pathways involved in adipogenesis (Wang et al., 2007; Zhang et al., 2017; Cui et al., 2018; Na et al., 2018; Wang et al., 2019). Similarly, other chicken lines continuously under divergent selection for abdominal fat content (French lines) or body weights (Virginia lines) were also constructed, to study the molecular genetics of fat deposition or growth rate (Baéza and Le Bihan-Duval, 2013; Siegel, 2014). However, to our knowledge, analysis on gene network construction and functional comparison on chicken adipogenesis is very limited.
In the present study, we constructed gene co-expression networks for four different transcriptome profiling datasets collected on abdominal fat tissues or isolated preadipocytes from chicken lines under divergent selection for adiposity, and compared the structural characteristics of the obtained network modules, to see if common molecular programs exist for adipogenesis. We identified important gene modules for adipogenesis, and interestingly, discovered common gene modules shared by different chicken lines, too. Our results provide evidences that even though genetic underpinnings of adiposity in chickens are complex, common molecular features could still exist, which renders novel insights on animal breeding practices, and also genetic investigation on human obesity.
Materials and Methods
Ethics Statement
All animal work was conducted in compliance with the recommended guidelines described in the Guide for the Care and Use of Laboratory Animals, and was approved by the Animal Care and Use Committee of Hubei Province, China (YZU-2018-0031).
Animals
Two chicken lines under divergent selection for abdominal fat content were used in the present study. One chicken line was from the Northeast Agricultural University broiler lines divergently selected for abdominal fat contents (NEAUHLF) since 1996, using abdominal fat percentage (AFP) and plasma very low-density lipoprotein (VLDL) concentration as selection criteria (Guo et al., 2011). These lines were developed from a common base grandsire line, Arbor Acres. Birds were bi-directionally selected for abdominal fat content, whereas bodyweights were kept the same for both lines. Abdominal fat percentages of the fat and lean broiler lines were significantly different from each other since the 4th generation, and detailed description on the selection procedure and housing conditions could be found in the previous report (Guo et al., 2011). The other French chicken line was from fat and lean chickens bred and raised at INRA UE1295 Pôle d’Expérimentation Avicolede Tours, F–37380 Nouzilly, France, as described previously (Resnyk et al., 2013, 2017).
Transcriptome Datasets
In order to explore the genes and molecular pathways affecting adipose tissue growth and development, we analyzed four sets of transcriptome data in the present study: one from the cellular perspective (in vitro preadipocyte differentiation), and the other three from adipose tissue perspectives (datasets downloaded from the public database). At the cellular level, microarray data on the differentiation of abdominal preadipocytes of NEAUHLF were collected. We briefly described the whole procedure as follows. Male birds at 10 days of age from the lean and fat lines at the 13th generation were selected, which were offsprings of the families with the highest and lowest AFP according to their slaughtered sib information, respectively. For each line, abdominal adipose tissues of 15 male birds were excised and pooled together. All pooled abdominal adipose tissues were then digested, filtrated, centrifuged, thus allowed the separation of floating adipocytes from the preadipocytes. The stage when the primary preadipocytes were collected was defined as the 0001 time point. And preadipocytes were passaged once, and harvested when 50% confluent (termed as -12 h) and 95% confluent (termed as 0 h). From the time point 0 h onward, preadipocytes were divided into two groups, the control group, and the group treated with 160 μM oleate. Cultured preadipocytes were divided into four groups, including preadipocytes from the lean and fat lines cultured without oleate treatment (LC and FC), or with oleate (LO and FO), respectively. Then preadipocytes in both groups were cultured continuously for 120 h, and samples were collected at four time points (12, 24, 72, and 120 h). Preadipocytes at each time point were collected, and total RNA was prepared from preadipocytes using Trizol reagent (Invitrogen) following the manufacturer’s instructions. Then, cDNA was prepared by oligo(dT)-primed reverse transcription. Labeled cRNA probes were prepared using an IVT Labeling Kit (Affymetrix, Inc.) according to the manufacturer’s protocol. The cRNA were fragmented, heated, loaded onto the Affymetrix probe array cartridge (Affymetrix, Inc.), and then hybridized, washed, and scanned at 560 nm using a confocal scanner. Raw data sets were normalized with Microarray Suite 5.0 (MAS5) and limma package in the R statistical environment (Ritchie et al., 2015; R Core Team, 2019). The raw microarray data of chicken preadipocyte differentiation were deposited into the National Center for Biotechnology Information (NCBI) Gene Expression Omnibus, and could be accessed through GEO Series (accession number: GSE51330).
At the tissue level, three datasets for abdominal fat tissue were accessed from the public domain, one microarray dataset from NEAUHLF, and one RNA-seq and one microarray dataset from the French chicken line, respectively. For the NEAUHLF dataset, 10 birds of the 8th generation were chosen based on the AFP values (five had the highest AFP and the other five had the lowest), and were slaughtered at 7 weeks. Abdominal fat tissues were collected for RNA extraction, and submitted for the GeneChip Chicken Genome Array (Wang et al., 2007). Data were downloaded through GEO (accession number: GSE8010). For the French lines, RNA-Seq and microarray datasets were collected from abdominal adipose tissues isolated from 7-week-old fat and lean chickens, as described previously (Resnyk et al., 2013, 2017), and were downloaded from the public GEO database (accession numbers: GSE42980 and GSE37585), respectively.
Construction of Weighted Gene Co-expression Network and Identification of Significant Modules
To identify gene co-expression network modules associated with preadipocyte differentiation and adipose tissue growth and development in chickens, the weighted gene co-expression network analysis (WGCNA) was conducted by the WGCNA package in the R statistical environment based on the expression profile data for each gene (Langfelder and Horvath, 2008). We did not perform merging the RNA-seq and microarray data, since the four transcriptome profiling datasets were collected from different broiler lines based on different selection criteria, and of different tissue origins (abdominal fat and preadipocytes). Instead, we constructed the network separately for each of the four data sets, and then compared the network modules obtained, to identify the common and line-specific modules associated with the adipogenesis trait, following the methods as described previously (Langfelder and Horvath, 2008). First, the raw gene expression data were pre-processed. We normalized and formatted the gene expression data into a data frame in the R statistical environment and filter out the null values of gene expression data. We performed sample clustering and gene expression analyses, respectively, and outliers (sample-wise or gene-/probe-wise) were detected and discarded. Second, soft threshold powers (β) were set first to fulfill the scale-free network assumption. In addition, we selected 5,000 genes with the highest connectivity for subsequent analysis based on the kRank, and calculated the weighted correlation values between genes to build the adjacency matrix. Third, the adjacency matrix was converted to a topological overlap matrix (TOM) to reduce noise and false correlation. Afterward, 1-TOM was calculated and treated as a biological important measure for network interconnectedness, and used as the distance measure for hierarchically clustered genes. Then, the dynamic tree cut algorithm was used to identify modules of co-expressed genes, and the module eigengenes (ME) were computed. Modules were merged based on the MEs, and labeled with different colors. Lastly, the specific modules were identified based on the module-trait relationship (pairwise Pearson’s correlation between ME and phenotypic trait values, and correlation coefficients >0.80). For abdominal fat tissue samples from the French lines, phenotypic traits were defined according to the bodyweights of the samples. As for the NEAUHLF tissue microarray samples, phenotypic traits included the bodyweights, the weights and percentages of abdominal fat tissues. For the preadipocyte microarray samples, we used directly their group properties, i.e., 1 and 0 for fat and lean chicken lines, respectively.
Functional Analysis of Modules
We next performed the functional analysis of important hub genes and network modules. The hub genes were identified in each module based on the module membership (MM) values and the positional importance of genes in the co-expression network (Langfelder and Horvath, 2008). The detected hub genes of specific modules and their relationship were visualized in Cytoscape (version 3.6.1). Modules for each constructed co-expression network were compared, based on their molecular function and pathway importance. Furthermore, functional enrichment analysis was performed for the specific modules based on Gene ontology (GO) enrichment and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses via DAVID (Database for Annotation Visualization and Integrated Discovery)1 (Huang et al., 2009). We then manually curated the identified significant pathways (P < 0.05). Based on the functional analysis results, we searched first common genes, and then common network sub-modules, to understand their relationship with adipogenesis.
Functional Validation of Network Modules
Genes of important roles in modules related to chicken fat development were selected and studied by Real-time Quantitative PCR Detecting System (qPCR). The immortalized broiler cell lines ICP1 (Wang et al., 2017) were cultured and differentiated by oleic acid treatment in vitro. Total RNAs from cells undergoing differentiation were collected and extracted by the TRIZOL method. Genomic DNA was removed by the PrimeScript RT reagent Kit with gDNA Eraser. Using the SYBR Green method, the qPCR experiment system was set up on ABI7500 as follows: the total volume was 10 μL, including the forward and reverse primers (0.2 μL each), the cDNA template (1 μL), Fast Start Universal SYBR Green Master (ROX, 2×) (5 μL), and ddH2O (3.6 μL). The reaction conditions were as follows: 40 cycles of pre-denaturation at 95°C for 10 min, then denaturation at 95°C for 15 s, and lastly, extension at 60°C for 60 s. The TATA-box binding protein (TBP) was used as the reference gene, and data were analyzed by the 2–ΔCt method (Livak and Schmittgen, 2001).
Results
Gene Co-expression Network Construction
To utilize effectively and accurately the public gene expression datasets, we began with standardized data preprocessing for each datasets. For one of the two datasets from the NEAUHLF chicken lines (the preadipocyte differentiation dataset), no outliers were detected, i.e., all 66 samples and 38,536 probes were kept (Supplementary Figure S1A); for the tissue expression dataset, two samples and probes were discarded (8 samples and 18,914 probes remained), respectively (Supplementary Figure S1B). As for the French datasets, none of the samples, but 2,649 probes, were discarded for the tissue RNA-seq dataset (24 samples and 15,286 transcripts remained) (Supplementary Figure S1C). And one sample, but no probe, was discarded for the tissue microarray (23 samples and 17,435 probes kept) (Supplementary Figure S1D).
After the model comparison (the scale free topology model vs. the fitted model) and by considering R2 >0.8, we found soft thresholds (β = 9, 12, 12, 10) for the preadipocyte and tissue microarrays of NEAUHLFF, and the RNA-seq and tissue microarray datasets of French lines, respectively (Supplementary Figure S2). Then, for the construction of co-expression networks, we selected 5,000 genes with the most expression variation, and transformed the adjacency matrix into topology matrix. Dynamic tree cut algorithm was used to merge all modules (correlation coefficients >0.75, and lowest number of genes in the module set at 30) (Figure 1). For the four datasets (the preadipocyte and tissue microarrays of NEAUHLF lines, and the tissue RNA-seq and microarray datasets of French lines), 13, 6, 12, and 17 modules were detected, respectively (Figures 1A–D). And specifically, modules containing the largest and least number of genes were Turquoise (1,173) and Darkgray (42) (Figure 1A), Blue (3,457) and Purple (163) (Figure 1B), Turquoise (1,173) and Tan (33) (Figure 1C), Lightgreen (838) and Mediumorchid (31) (Figure 1D), respectively.
Figure 1. Clustering of co-expression modules. Upper panel: genes were clustered into different groups, and assigned to network modules after dynamic tree cut and merging analyses (Lower panel). (A) NEAUHLF: preadipocyte microarray. (B) NEAUHLF: adipose tissue microarray. (C) French: adipose tissue RNA-seq. (D) French: adipose tissue microarray.
Identification of Modules of Interests
We obtained the module-trait heatmaps by calculating the correlation coefficients between modules and traits for the four datasets, respectively (Figure 2). First, for the NEAUHLF preadipocyte dataset, the Turquoise module was highly positively correlated with the preadipocyte differentiation group after oleate treatment for both fat and lean lines, but negatively correlated with the control groups (Figure 2A). Moreover, the Turquoise module was significantly positively correlated with the induced differentiation group of the fat line, but the early differentiation group showed stronger correlation than the late stage group. As for the lean line, the Turquoise group was significantly negatively correlated with the control group, but showing stronger correlation for the late stage than the early stage. So, genes in the Turquoise module could promote the differentiation, but inhibit the proliferation of preadipocytes. The Lightgreen module, in contrast, was positively correlated with both the induced differentiation and the control groups for the lean chicken line, but negatively correlated with all treatment and control groups of the fat line, indicating that the Lightgreen module could inhibit preadipocyte proliferation, and enhance lipid catabolism. Moreover, we found that the Darkgreen module was positively correlated with the induced differentiation group of the fat line, but negatively correlated with the control group of the lean line. The degree of correlation appeared in a progressively increasing or decreasing manner for the fat and lean lines, respectively. Thus, genes in the Darkgreen module could enhance the preadipocyte differentiation, and promote adipogenesis.
Figure 2. Module-trait heatmaps. Gene modules and trait relationships established by correlation analyses. (A) NEAUHLF: preadipocyte microarray. (B) NEAUHLF: adipose tissue microarray. (C) French: adipose tissue RNA-seq. (D) French: adipose tissue microarray.
As for the NEAUHLF tissue microarray dataset, we didn’t find any strong correlation with the traits for most of the modules (Figure 2B). Only the Blue module correlated strongly with the bodyweights, and relatively strongly with both the abdominal fat rate and weight, too. Whereas for the French tissue RNA-seq dataset, the Blue module correlated positively with the lean line, but negatively with the fat line (all most significantly), respectively. Genes with negative effects on fat deposition might be contained in the Blue module (Figure 2C).
For the French tissue microarray dataset (Figure 2D), the Turquoise and Orange modules were positively correlated with the bodyweight at 1 week of age, but negatively correlated with those at 11 and 9 weeks of age, respectively. In contrast, the Darkseagreen4 and Darkorange modules showed the opposite trend, correlated negatively with the bodyweight at 1 week of age, but positively with bodyweights at 9 and 11 weeks of ages (all most significantly), respectively. Genes in these modules could have potential effects on fat deposition. Thus, our module-trait correlation analyses pinpoint network modules potentially underlying adipogenesis or the growth and development of chicken adipose tissue.
Hub Genes in Important Sub-Modules
Hub gene occupies central position, and is of vital function in the gene network, which is determined by their connectedness with other genes in the network. Our previous analyses helped find nine functional modules of importance to adipose tissue growth and development, after we built the network based on the threshold of correlation coefficients (>0.85) (Table 1). Then, we selected 3–6 hub genes from each module. Functional exploration revealed that these genes were related to adipogenesis, cell cycle, inflammation and protein synthesis. In addition, we extracted and plotted the gene expression values of hub genes, which showed that different hub genes in the same module were of similar expression patterns (Supplementary Figure S3).
Table 1. Hub genes in modules identified to be related to chicken fat development.
Common Modules Identified by Network Comparison
We continued to build the gene co-expression sub-networks around the discovered hub genes, in an attempt to show the importance and implication of these genes on adipogenesis. For the four datasets, 3, 1, 1, and 4 sub-networks were found, respectively (Figure 3). After functional enrichment analyses and literature retrieval on central hub genes (such as ROCK2, RRM2B, PTN, CAV2, FARSB, SH3GLB1, FABP5, and FADS2), interestingly, we found great similarity for biological processes and molecular pathways (Supplementary Table S1). For many modules constructed, development and regulation of adipose tissue and metabolic pathways were found to be enriched, such as lipid storage, lipid modification, negative regulation of lipid synthesis, fatty acid metabolism, PPAR and insulin signaling pathways. In addition, after GO analysis, modules were found to be enriched in a variety of pathways, such as protein hydrolysis, transport and localization, cytoskeleton, apoptosis and programmed cell death. KEGG pathway analysis also identified two pathways, focal adhesion and proteasome metabolism (Supplementary Table S1).
Figure 3. Subnetworks constructed for hub genes. Direct gene neighbors of the hub genes were extracted for building subnetworks. Four subnetworks were built. (A) NEAUHLF: preadipocyte microarray. (B) NEAUHLF: adipose tissue microarray. (C) French: adipose tissue RNA-seq. (D) French: adipose tissue microarray.
Furthermore, based on the functional consistency between identified important modules, common gene modules were found (e.g., Turquoise, Blue, Blue, and Turquoise modules for the four datasets, the NEAUHLF preadipocyte and tissue microarrays, and the French tissue RNA-seq and microarray datasets, respectively, even though labeled with different color codes) (Figure 4). These common modules were enriched in the establishment of cytoskeleton, and protein metabolic transport and localization. Additional analyses found many modules with similar functional enrichment, such as cell apoptosis and programmed cell death (Supplementary Table S2). In contrast, modules specific to chicken lines and having their specific functionality to chicken adipogenesis were also found. Thus, in chicken lines divergently selected for abdominal fat contents, though of different genetic backgrounds, common gene modules with similar functionality were detected, and could potentially play conserved roles in the growth and development of adipose tissues.
Figure 4. Common gene modules discovered for different chicken lines.
Validation of the Lysosome Module
After functional enrichment analysis, the lysosome pathway was among those pathways discovered to be potentially important for chicken adipose tissue growth and development. Lysosome has a vital role in lipid metabolism, and we and others showed previously that the lysosome pathway is fundamental to chicken adipogenesis (Thelen and Zoncu, 2017; Data not shown). Here we selected TFEB (the master regulator of lysosome biogenesis and also involved in lipid metabolism), and also its downstream genes (LAMP1, CTSA, CTSB), to examine their expression patterns during preadipocyte differentiation in the preadipocyte cell line (ICP1) (Supplementary Figure S4). The expression dynamics of TFEB and its downstream genes exhibited a similar pattern, and highly correlated with each other. Furthermore, obvious changes along with the differentiation of preadipocytes in broilers could be seen, since significant differences existed between these genes at different time points (Figure 5). Thus, the lysosome pathway may play a regulatory role in the differentiation of chicken preadipocytes.
Figure 5. The expression and significance test of TFEB and its downstream genes during induced differentiation of ICP1 preadipocyte line. (A) TFEB. (B) LAMP1. (C) CTSA. (D) CTSB. Mean significance at *P < 0.05 and **P < 0.01 levels.
Discussion
In the current study, gene co-expression networks were built for two different chicken lines both under divergent selection for abdominal fat content, and common network modules with similar functionality were discovered to be of potential importance for adiposity. Our findings could provide novel insights into the genetic basis of complex traits, and help understand the outcomes of intensive breeding practices in modern animal production systems.
The recently developed high-throughput RNA-seq technology propels the study on the detection of differentially expressed genes by statistical modeling and analysis on transcriptome profiling data. However, the co-linearity and confounding factors embedded within the gene expression data (usually small number of samples, but large number of parameters, i.e., the n << p problem) could not be discerned by these normal statistical methods (Barabási and Oltvai, 2004; Ghaemi et al., 2019). To cope with the high-dimensional data analysis question induced by large number of genes, methods for the construction of gene networks were developed, such as the gene co-expression network (Langfelder and Horvath, 2008; Banf and Rhee, 2017; Mangul et al., 2019), and have been effectively used in identifying the topological relatedness and interaction of genes in the biological systems of interests (Lee et al., 2016; Emilsson et al., 2018). Nevertheless, the spatio-temporal and heterogeneous characteristics of the transcriptome data require the deployment of novel statistical methods (e.g., machine learning), to effectively solve the statistical and modeling issues related to network dynamics and heterogeneity (Camacho et al., 2018).
We identified network modules significantly associated with fat deposition in chicken lines with different genetic backgrounds, but all under divergent selection for the same trait of interest, abdominal fat content (Guo et al., 2011; Resnyk et al., 2013, 2017). Previously, after differential expression analyses on the same datasets, genes involved in fatty acid metabolism and PPAR signaling pathways were found, such as PPARG, and its direct targets, SCD, ACSL1, and DGAT2 (Guo et al., 2011; Resnyk et al., 2013, 2017), which potentially explained the underlying differences of fat deposition in these chicken lines. However, with gene co-expression network analysis, we found network modules containing interesting genes in molecular pathways (fatty acid metabolism, PPAR and insulin, cytoskeleton, and protein synthesis, etc.) fundamental to adipose tissue growth and development, such as FABP5, FADS2. The insulin pathway was found to be vital to fatty acid content composition in cattle (de Oliveira et al., 2019), and FABP5 and FADS2 are well-known regulatory genes on fatty acid metabolism (Li et al., 2019; Xing et al., 2019). Cytoskeleton has to reorganize with adipose tissue remodeling and expansion, and recently cytoskeletal transgelin 2 was proved to be associated with preadipocyte proliferation and differentiation (Ortega et al., 2019). Protein synthesis is regulated mainly by the master growth regulator mTOR (Thoreen et al., 2012), which plays a vital role in lipid metabolism (Caron et al., 2015). Further expression pattern analysis showed that these genes in the same module did have similar expression levels and dynamics, which could be transcriptionally regulated by a common set of transcription factors (Gachon et al., 2018; Klemm et al., 2019). These discovered genes and signaling pathways seem to be common and have important conserved functions for adipogenesis, which could be partially due to that divergent selection put more forces on gene modules of common fundamental roles in adipogenesis, and were picked up by our network module comparison analysis. Further detailed investigation on the transcriptional regulation of genes in these identified network modules could render novel insights into the underlying molecular regulatory mechanisms of these chicken lines.
Integrated genomics and network science methods are widely employed in systems biology, to fully utilize the large volume of genomics and biological data (Lee et al., 2016; Emilsson et al., 2018; Ghaemi et al., 2019). We herein detected common molecular network and functional pathways involved in abdominal fat deposition after network construction and comparison between different chicken lines. It’s a common phenomenon that similar traits of animals in different populations could be under convergent selection, i.e., different genes but pathways of similar molecular functions are selected by evolutionary forces or artificial selective pressure, such as convergent selection signatures in sheep and goat (Alberto et al., 2018), and the frizzle phenotype in chickens (Dong et al., 2018; Guo et al., 2018). Here, our network analysis also showed that, though different chicken lines were divergently selected for abdominal fat content, network modules of similar molecular functionality were detected. These results could provide further insights on the genetics of complex trait, including human diseases.
Conclusion
Gene co-expression networks were constructed for different chicken lines under divergent selection for adiposity. Common sub-modules of similar functionality for chicken fat deposition were identified after gene functional enrichment analysis and network comparison. Our findings indicate that even in different chicken lines, common molecular pathways could be underlying the growth and development of adipose tissue.
Data Availability Statement
Publicly available datasets were analyzed in this study. This data can be found here: GSE51330, GSE8010, GSE42980, and GSE37585.
Ethics Statement
The animal study was reviewed and approved by Laboratory Animal Management Committee of Yangtze University.
Author Contributions
ZG analyzed data and wrote the manuscript. ZG and RD performed the experiment. XZ, YW, and YC participated in the experiment. Z-QD and C-XY designed the experiment, analyzed data, and wrote the manuscript.
Funding
The funding support was provided by the Chinese National Natural Science Foundation (Grant No. 31472088).
Conflict of Interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Acknowledgments
We appreciate the help from other lab members, and also the Poultry Research Group at Northeast Agricultural University.
Supplementary Material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.2020.00537/full#supplementary-material
FIGURE S1 | Cluster analyses of samples and gene expression data to detect outliers. No outliers found for (A) (NEAUHLF: preadipocyte microarray), whereas two samples and two genes were discarded for (B) (NEAUHLF: adipose tissue microarray). As for the French datasets, 2,649 probes were discarded for the adipose tissue RNA-seq dataset (C) and one sample outlier was discarded for the adipose tissue microarray (D).
FIGURE S2 | Identification of soft thresholds. The soft thresholds β = 9, 12, 12, 10 for preadipocyte (A) and adipose tissue (B) microarrays of NEAUHLFF, and adipose tissue RNA-seq (C) and microarray (D) of French lines.
FIGURE S3 | Gene expression profiles of hub genes in the important modules. Consistent and highly correlated gene expression patterns were identified for hub genes in the network sub-modules, for the four datasets, respectively.
FIGURE S4 | Culture and induced differentiation of chicken preadipocyte cell line. Oil red oxygen staining (A) and gene expression trend of differentiation markers of adipocytes in broilers (B).
TABLE S1 | Enrichment analysis of modules.
TABLE S2 | Summary of hub genes in modules related to fat development.
Footnotes
References
Ahrends, R., Ota, A., Kovary, K. M., Kudo, T., Park, B. O., and Teruel, M. N. (2014). Controlling low rates of cell differentiation through noise and ultrahigh feedback. Science 344, 1384–1389. doi: 10.1126/science.1252079
Alberto, F. J., Boyer, F., Orozco-Terwengel, P., Streeter, I., Servin, B., de Villemereuil, P., et al. (2018). Convergent genomic signatures of domestication in sheep and goats. Nat. Commun. 9:813. doi: 10.1038/s41467-018-03206-y
Alexander, P. A., Naval-Sanchez, M., Porto-Neto, L. R., Ferraz, J. B. S., Reverter, A., and Fukumasu, H. (2019). Systems biology reveals NR2F6 and TGFB1 as key regulators of feed efficiency in beef cattle. Front. Genet. 10:230. doi: 10.3389/fgene.2019.00230
Baéza, E., and Le Bihan-Duval, E. (2013). Chicken lines divergent for low or high abdominal fat deposition: a relevant model to study the regulation of energy metabolism. Animal 7, 965–973. doi: 10.1017/S1751731113000153
Banf, M., and Rhee, S. Y. (2017). Computational inference of gene regulatory networks: approaches, limitations and opportunities. Biochim. Biophys. Acta 1860, 41–52. doi: 10.1016/j.bbagrm.2016.09.003
Barabási, A. L., and Oltvai, Z. N. (2004). Network biology: understanding the cell’s functional organization. Nat. Rev. Genet. 5, 101–113. doi: 10.1038/nrg1272
Beiki, H., Nejati-Javaremi, A., Pakdel, A., Masoudi-Nejad, A., Hu, Z. L., and Reecy, J. M. (2016). Large-scale gene co-expression network as a source of functional annotation for cattle genes. BMC Genomics 17:846. doi: 10.1186/s12864-016-3176-2
Beiki, H., Pakdel, A., Javaremi, A. N., Masoudi-Nejad, A., and Reecy, J. M. (2018). Cattle infection response network and its functional modules. BMC Immunol. 19:2. doi: 10.1186/s12865-017-0238-4
Blüher, M. (2019). Obesity: global epidemiology and pathogenesis. Nat. Rev. Endocrinol. 15, 288–298. doi: 10.1038/s41574-019-0176-8
Brial, F., Le Lay, A., Hedjazi, L., Tsang, T., Fearnside, J. F., Otto, G. W., et al. (2019). Systems genetics of hepatic metabolome reveals octopamine as a target for non-alcoholic fatty liver disease treatment. Sci. Rep. 9:3656. doi: 10.1038/s41598-019-40153-0
Buchanan, J. W., Reecy, J. M., Garrick, D. J., Duan, Q., Beitz, D. C., Koltes, J. E., et al. (2016). Deriving gene networks from SNP associated with triacylglycerol and phospholipid fatty acid fractions from ribeyes of Angus cattle. Front. Genet. 7:116. doi: 10.3389/fgene.2016.00116
Camacho, D. M., Collins, K. M., Powers, R. K., Costello, J. C., and Collins, J. J. (2018). Next-generation machine learning for biological networks. Cell 173, 1581–1592. doi: 10.1016/j.cell.2018.05.015
Caron, A., Richard, D., and Laplante, M. (2015). The roles of mTOR complexes in lipid metabolism. Annu. Rev. Nutr. 35, 321–348. doi: 10.1146/annurev-nutr-071714-034355
Chen, Y., Zhu, J., Lum, P. Y., Yang, X., Pinto, S., MacNeil, D. J., et al. (2008). Variations in DNA elucidate molecular networks that cause disease. Nature 452, 429–435. doi: 10.1038/nature06757
Cohen, P., and Spiegelman, B. M. (2016). Cell biology of fat storage. Mol. Biol. Cell 27, 2523–2527. doi: 10.1091/mbc.E15-10-0749
Cui, T., Xing, T., Huang, J., Mu, F., Jin, Y., You, X., et al. (2018). Nuclear respiratory factor 1 negatively regulates the P1 promoter of the peroxisome proliferator-activated receptor-γ gene and inhibits chicken adipogenesis. Front. Physiol. 9:1823. doi: 10.3389/fphys.2018.01823
de Oliveira, P. S. N., Coutinho, L. L., Cesar, A. S. M., Diniz, W. J. D. S., de Souza, M. M., Andrade, B. G., et al. (2019). Co-expression networks reveal potential regulatory roles of miRNAs in fatty acid composition of Nelore cattle. Front. Genet. 10:651. doi: 10.3389/fgene.2019.00651
de Oliveira, P. S. N., Coutinho, L. L., Tizioto, P. C., Cesar, A. S. M., de Oliveira, G. B., Diniz, W. J. D. S., et al. (2018). An integrative transcriptome analysis indicates regulatory mRNA-miRNA networks for residual feed intake in Nelore cattle. Sci. Rep. 8:17072. doi: 10.1038/s41598-018-35315-5
Dong, J., He, C., Wang, Z., Li, Y., Li, S., Tao, L., et al. (2018). A novel deletion in KRT75L4 mediates the frizzle trait in a Chinese indigenous chicken. Genet. Select. Evol. 50:68. doi: 10.1186/s12711-018-0441-7
Dyar, K. A., Lutter, D., Artati, A., Ceglia, N. J., Liu, Y., Armenta, D., et al. (2018). Atlas of circadian metabolism reveals system-wide coordination and communication between clocks. Cell 174, 1571–1585. doi: 10.1016/j.cell.2018.08.042
Emilsson, V., Ilkov, M., Lamb, J. R., Finkel, N., Gudmundsson, E. F., Pitts, R., et al. (2018). Co-regulatory networks of human serum proteins link genetics to disease. Science 361, 769–773. doi: 10.1126/science.aaq1327
Fortes, M. R., Reverter, A., Zhang, Y., Collis, E., Nagaraj, S. H., Jonsson, N. N., et al. (2010). Association weight matrix for the genetic dissection of puberty in beef cattle. Proc. Natl. Acad. Sci. U.S.A. 107, 13642–13647. doi: 10.1073/pnas.1002044107
Gachon, F., Yeung, J., and Naef, F. (2018). Cross-regulatory circuits linking inflammation, high-fat diet, and the circadian clock. Genes Dev. 32, 1359–1360. doi: 10.1101/gad.320911.118
Ghaben, A. L., and Scherer, P. E. (2019). Adipogenesis and metabolic health. Nat. Rev. Mol. Cell Biol. 20, 242–258. doi: 10.1038/s41580-018-0093-z
Ghaemi, M. S., DiGiulio, D. B., Contrepois, K., Callahan, B., Ngo, T. T. M., Lee-McMullen, B., et al. (2019). Multiomics modeling of the immunome, transcriptome, microbiome, proteome and metabolome adaptations during human pregnancy. Bioinformatics 35, 95–103. doi: 10.1093/bioinformatics/bty537
Gonçalves, T. M., de Almeida Regitano, L. C., Koltes, J. E., Cesar, A. S. M., da Silva Andrade, S. C., Mourão, G. B., et al. (2018). Gene co-expression analysis indicates potential pathways and regulators of beef tenderness in Nellore cattle. Front. Genet. 9:441. doi: 10.3389/fgene.2018.00441
Guan, D., Xiong, Y., Borck, P. C., Jang, C., Doulias, P. T., Papazyan, R., et al. (2018). Diet-induced circadian enhancer remodeling synchronizes opposing hepatic lipid metabolic processes. Cell 174, 831–842. doi: 10.1016/j.cell.2018.06.031
Guo, L., Sun, B., Shang, Z., Leng, L., Wang, Y., Wang, N., et al. (2011). Comparison of adipose tissue cellularity in chicken lines divergently selected for fatness. Poult. Sci. 90, 2024–2034. doi: 10.3382/ps.2010-00863
Guo, X., Li, Y. Q., Wang, M. S., Wang, Z. B., Zhang, Q., Shao, Y., et al. (2018). A parallel mechanism underlying frizzle in domestic chickens. J. Mol. Cell Biol. 10, 589–591. doi: 10.1093/jmcb/mjy037
Huang, D. W., Sherman, B. T., and Lempicki, R. A. (2009). Systematic and integrative analysis of large gene lists using DAVID Bioinformatics Resources. Nat. Protoc. 4, 44–57. doi: 10.1038/nprot.2008.211
Jeffery, E., Wing, A., Holtrup, B., Sebo, Z., Kaplan, J. L., Saavedra-Peña, R., et al. (2016). The adipose tissue microenvironment regulates depot-specific adipogenesis in obesity. Cell Metab. 24, 142–150. doi: 10.1016/j.cmet.2016.05.012
Justice, A. E., Karaderi, T., Highland, H. M., Young, K. L., Graff, M., Lu, Y., et al. (2019). Protein-coding variants implicate novel genes related to lipid homeostasis contributing to body-fat distribution. Nat. Genet. 51, 452–469. doi: 10.1038/s41588-018-0334-2
Klemm, S. L., Shipony, Z., and Greenleaf, W. J. (2019). Chromatin accessibility and the regulatory epigenome. Nat. Rev. Genet. 20, 207–220. doi: 10.1038/s41576-018-0089-8
Langfelder, P., and Horvath, S. (2008). WGCNA: an R package for weighted correlation network analysis. BMC Bioinformatics 9:559. doi: 10.1186/1471-2105-9-559
Lee, S., Zhang, C., Kilicarslan, M., Piening, B. D., Bjornson, E., Hallström, B. M., et al. (2016). Integrated network analysis reveals an association between plasma mannose levels and insulin resistance. Cell Metab. 24, 172–184. doi: 10.1016/j.cmet.2016.05.026
Lee, Y. H., Mottillo, E. P., and Granneman, J. G. (2014). Adipose tissue plasticity from WAT to BAT and in between. Biochim. Biophys. Acta 1842, 358–369. doi: 10.1016/j.bbadis.2013.05.011
Li, X., Li, P., Wang, L., Zhang, M., and Gao, X. (2019). Lysine enhances the stimulation of fatty acids on milk fat synthesis via the GPRC6A-PI3K-FABP5 signaling in bovine mammary epithelial cells. J. Agric. Food Chem. 67, 7005–7015. doi: 10.1021/acs.jafc.9b02160
Livak, K. J., and Schmittgen, T. D. (2001). Analysis of relative gene expression data using real-time quantitative PCR and the 2(-Delta Delta C(T)) Method. Methods 25, 402–408. doi: 10.1006/meth.2001.1262
Luong, Q., Huang, J., and Lee, K. Y. (2019). Deciphering white adipose tissue heterogeneity. Biology 8:23. doi: 10.3390/biology8020023
Mangul, S., Martin, L. S., Hill, B. L., Lam, A. K., Distler, M. G., Zelikovsky, A., et al. (2019). Systematic benchmarking of omics computational tools. Nat. Commun. 10:1393. doi: 10.1038/s41467-019-09406-4
Mateescu, R. G., Garrick, D. J., and Reecy, J. M. (2017). Network analysis reveals putative genes affecting meat quality in Angus cattle. Front. Genet. 8:171. doi: 10.3389/fgene.2017.00171
Na, W., Wu, Y. Y., Gong, P. F., Wu, C. Y., Cheng, B. H., Wang, Y. X., et al. (2018). Embryonic transcriptome and proteome analyses on hepatic lipid metabolism in chickens divergently selected for abdominal fat content. BMC Genomics 19:384. doi: 10.1186/s12864-018-4776-9
Ortega, F. J., Moreno-Navarrete, J. M., Mercader, J. M., Gomez-Serrano, M., Garcia-Santos, E., Latorre, J., et al. (2019). Cytoskeletal transgelin 2 contributes to gender-dependent adipose tissue expandability and immune function. FASEB J. 33, 9656–9671. doi: 10.1096/fj.201900479R
Pellegrinelli, V., Peirce, V. J., Howard, L., Virtue, S., Türei, D., Senzacqua, M., et al. (2018). Adipocyte-secreted BMP8b mediates adrenergic-induced remodeling of the neuro-vascular network in adipose tissue. Nat. Commun. 9:4974. doi: 10.1038/s41467-018-07453-x
R Core Team (2019). R: A Language and Environment for Statistical Computing. Vienna: R Foundation for Statistical Computing.
Ranola, J. M., Langfelder, P., Lange, K., and Horvath, S. (2013). Cluster and propensity based approximation of a network. BMC Syst. Biol. 7:21. doi: 10.1186/1752-0509-7-21
Resnyk, C. W., Carré, W., Wang, X., Porter, T. E., Simon, J., Le Bihan-Duval, E., et al. (2013). Transcriptional analysis of abdominal fat in genetically fat and lean chickens reveals adipokines, lipogenic genes and a link between hemostasis and leanness. BMC Genomics 14:557. doi: 10.1186/1471-2164-14-557
Resnyk, C. W., Carré, W., Wang, X., Porter, T. E., Simon, J., Le Bihan-Duval, E., et al. (2017). Transcriptional analysis of abdominal fat in chickens divergently selected on bodyweight at two ages reveals novel mechanisms controlling adiposity: validating visceral adipose tissue as a dynamic endocrine and metabolic organ. BMC Genomics 18:626. doi: 10.1186/s12864-017-4035-5
Ritchie, M. E., Phipson, B., Wu, D., Hu, Y., Law, C. W., Shi, W., et al. (2015). limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res. 43:e47. doi: 10.1093/nar/gkv007
Sebo, Z. L., and Rodeheffer, M. S. (2019). Assembling the adipose organ: adipocyte lineage segregation and adipogenesis in vivo. Development 146:dev172098. doi: 10.1242/dev.172098
Siegel, P. B. (2014). Evolution of the modern broiler and feed efficiency. Annu. Rev. Anim. Biosci. 2, 375–385. doi: 10.1146/annurev-animal-022513-114132
Tang, Q. Q., and Lane, M. D. (2012). Adipogenesis: from stem cell to adipocyte. Annu. Rev. Biochem. 81, 715–736. doi: 10.1146/annurev-biochem-052110-115718
Thelen, A. M., and Zoncu, R. (2017). Emerging roles for the lysosome in lipid metabolism. Trends Cell Biol. 27, 833–850. doi: 10.1016/j.tcb.2017.07.006
Thoreen, C. C., Chantranupong, L., Keys, H. R., Wang, T., Gray, N. S., and Sabatini, D. M. (2012). A unifying model for mTORC1-mediated regulation of mRNA translation. Nature 485, 109–113. doi: 10.1038/nature11083
van der Klaauw, A. A., Croizier, S., Mendes de Oliveira, E., Stadler, L. K. J., Park, S., Kong, Y., et al. (2019). Human semaphorin 3 variants link melanocortin circuit development and energy balance. Cell 176, 729–742. doi: 10.1016/j.cell.2018.12.009
Wang, H. B., Li, H., Wang, Q. G., Zhang, X. Y., Wang, S. Z., Wang, Y. X., et al. (2007). Profiling of chicken adipose tissue gene expression by genome array. BMC Genomics 8:193. doi: 10.1186/1471-2164-8-193
Wang, W., Wu, K., Jia, M., Sun, S., Kang, L., and Zhang, Q. (2018). Dynamic changes in the global microRNAome and transcriptome identify key nodes associated with ovarian development in chickens. Front. Genet. 9:491. doi: 10.3389/fgene.2018.00491
Wang, W., Zhang, T., Wu, C., Wang, S., Wang, Y., Li, H., et al. (2017). Immortalization of chicken preadipocytes by retroviral transduction of chicken TERT and TR. PLoS One 12:e0177348. doi: 10.1371/journal.pone.0177348
Wang, Z. B., Du, Z. Q., Na, W., Jing, J. H., Li, Y. M., Leng, L., et al. (2019). Production of transgenic broilers by non-viral vectors via optimizing egg windowing and screening transgenic roosters. Poult. Sci. 98, 430–439. doi: 10.3382/ps/pey321
Xing, K., Wang, K., Ao, H., Chen, S., Tan, Z., Wang, Y., et al. (2019). Comparative adipose transcriptome analysis digs out genes related to fat deposition in two pig breeds. Sci. Rep. 9:12925. doi: 10.1038/s41598-019-49548-5
Yang, X., Deignan, J. L., Qi, H., Zhu, J., Qian, S., Zhong, J., et al. (2009). Validation of candidate causal genes for obesity that affect shared metabolic pathways and networks. Nat. Genet. 41, 415–423. doi: 10.1038/ng.325
Zhang, H., Yu, J. Q., Yang, L. L., Kramer, L. M., Zhang, X. Y., Na, W., et al. (2017). Identification of genome-wide SNP-SNP interactions associated with important traits in chicken. BMC Genomics 18:892. doi: 10.1186/s12864-017-4252-y
Keywords: chicken, WGCNA, gene network, modules, fat deposition
Citation: Gao Z, Ding R, Zhai X, Wang Y, Chen Y, Yang C-X and Du Z-Q (2020) Common Gene Modules Identified for Chicken Adiposity by Network Construction and Comparison. Front. Genet. 11:537. doi: 10.3389/fgene.2020.00537
Received: 21 August 2019; Accepted: 04 May 2020;
Published: 29 May 2020.
Edited by:
Xiquan Zhang, South China Agricultural University, ChinaReviewed by:
Zhe Zhang, South China Agricultural University, ChinaLingyang Xu, Chinese Academy of Agricultural Sciences, China
Copyright © 2020 Gao, Ding, Zhai, Wang, Chen, Yang and Du. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Cai-Xia Yang, caixiayang@yangtzeu.edu.cn; Zhi-Qiang Du, zhqdu@yangtzeu.edu.cn