Consensus module analysis of abdominal fat deposition across multiple broiler lines

Despite several RNA-Seq and microarray studies on differentially expressed genes (DEGs) between high- and low-abdominal fat deposition in different broiler lines, to our knowledge, gene coexpression analysis across multiple broiler lines has rarely been reported. Here, we constructed a consensus gene coexpression network focused on identifying consensus gene coexpression modules associated with abdominal fat deposition across multiple broiler lines using two public RNA-Seq datasets (GSE42980 and GSE49121). In the consensus gene coexpression network, we identified eight consensus modules significantly correlated with abdominal fat deposition across four broiler lines using the consensus module analysis function in the weighted gene coexpression network analysis (WGCNA) package. The eight consensus modules were moderately to strongly preserved in the abdominal fat RNA-Seq dataset of another broiler line (SRP058295). Furthermore, we identified 5462 DEGs between high- and low-abdominal fat lines (FL and LL) (GSE42980) and 6904 DEGs between high- and low-growth (HG and LG) (GSE49121), including 1828 overlapping DEGs with similar expression profiles in both datasets, which were clustered into eight consensus modules. Pyruvate metabolism, fatty acid metabolism, and steroid biosynthesis were significantly enriched in the green, yellow, and medium purple 3 consensus modules. The PPAR signaling pathway and adipocytokine signaling pathway were significantly enriched in the green and purple consensus modules. Autophagy, mitophagy, and lysosome were significantly enriched in the medium purple 3 and yellow consensus modules. Based on lipid metabolism pathways enriched in eight consensus modules and the overexpression of numerous lipogenic genes in both FL vs. LL and HG vs. LG, we hypothesize that more fatty acids, triacylglycerols (TAGs), and cholesterol might be synthesized in broilers with high abdominal fat than in broilers with low abdominal fat. According to autophagy, mitophagy, and lysosome enrichment in eight consensus modules, we inferred that autophagy might participate in broiler abdominal fat deposition. Altogether, these studies suggest eight consensus modules associated with abdominal fat deposition in broilers. Our study also provides an idea for investigating the molecular mechanism of abdominal fat deposition across multiple broiler lines.

low-abdominal fat deposition in different broiler lines using RNA-Seq or microarray technology. There were 1687 DEGs identified between high-and low-abdominal fat line (FL and LL) broilers with a 2.8-fold divergence in abdominal fat percentage (AFP) at 7 weeks of age [6]. A total of 2410 DEGs were identified in abdominal fat tissue between high-and low-growth-rate (HG and LG) broilers, with a 19.2-fold divergence in AFP at 7 weeks of age [7]. A total of 286 DEGs were identified between low-and high-feed efficiency (LFE and HFE) broilers, with a 1.6-fold divergence in AFP at 7 weeks of age [8]. A total of 230 DEGs were identified between the Northeast Agricultural University (NEAU) broiler lines divergently selected for abdominal fat content with a 1.9-fold divergence in AFP at 7 weeks of age [9]. In each study, AFP was significantly different between broilers with high and low abdominal fat deposition, but the DEGs were not the same. The different DEGs among these studies may be due to the differences in lines, sample sizes, analytical methods, experimental technologies, and other factors. The differences among these studies may disturb further research on abdominal fat deposition in broilers. Consensus modules are comprised of genes tightly coexpressed in both datasets, so consensus modules exhibit a degree of preservation between the two data networks [10,11]. Therefore, it is important to construct a consensus gene coexpression network across multiple broiler lines and to detect consensus modules correlated with abdominal fat deposition. However, to our knowledge, this work has rarely been reported in broilers.
In the current study, to explore the molecular mechanism of abdominal fat deposition across multiple broiler lines, we constructed a consensus gene coexpression network across GSE42980 (12 FL and 12 LL samples) and GSE49121 (8 HG and 8 LG samples) and detected consensus modules associated with abdominal fat deposition. Then, we tested the preservation of consensus modules in abdominal fat RNA-Seq data of another broiler line (SRP058295). The ClusterProfiler software package was used to predict Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways and Gene Ontology (GO) terms associated with the genes within consensus modules. Our study provides an idea for investigating abdominal fat deposition across multiple broiler lines.

RNA-Seq data collection
To reduce the influence of nongenetic factors, we collected public abdominal fat gene expression profile datasets of broilers under the same age, sex, and normal feed conditions. Broiler chickens were divergently selected over seven generations for either high (FL) or low (LL) abdominal fat at similar feed intake and body weight, and FL and LL broilers exhibited a 2.8-fold divergence in AFP at 7 weeks of age [6,12]. GSE42980 is an abdominal fat RNA-Seq dataset of male broilers at 7 weeks of age that consists of 12 FL and 12 LL samples [6]. Broiler chickens were divergently selected over 30 generations for either high (HG) or low (LG) growth, and HG and LG broilers were different in AFP (19.2fold) and body weight (3.2-fold) at 7 weeks of age [7]. GSE49121 is an abdominal fat RNA-Seq dataset of male broilers that consists of eight HG and eight LG samples. Because two pairwise broiler lines (FL and LL, HG and LG) are divergently selected over several generations and exhibit significant divergence in AFP at 7 weeks of age, the four broiler lines are good genetic models for abdominal fat deposition research. SRP058295 is an abdominal fat RNA-Seq dataset of male broilers [8]. Although broiler chickens from the SRP058295 dataset were not genetically selected, they exhibited divergence in AFP at 7 weeks of age from normal feed. In addition, all broilers in GSE42980, GSE49121, and SRP058295 were fed a normal diet. Therefore, the SRP058295 dataset was used to validate the preservation of consensus modules across two datasets. In summary, we successfully downloaded the GSE42980, GSE49121, and SRP058295 datasets for subsequent WGCNA analyses.

QC and mapping of RNA-Seq reads
The average number of raw reads across the samples within GSE42980, GSE49121, and SRP058295 was 54.97 M (Additional file 1: Table S1). To ensure correct results in further analyses, low-quality reads were filtered, and the average number of high-quality reads was 48.97. The average percentage of high-quality reads mapped to the chicken genome was 91.65% across all samples within the three datasets (Additional file 1: Table S1). The average percentage of high-quality reads mapping to exonic regions was 80.02% (Additional file 1: Table S1).

Hierarchical cluster analyses
In the current study, a hierarchical cluster of GSE42980 showed that 12 FL and 12 LL samples were fully separated (Additional file 2: Figure S1, a). A hierarchical cluster of GSE49121 showed that eight high-(HG) and eight low-abdominal fat deposition (LG) samples were fully separated (Additional file 2: Figure S1, b). Therefore, 24 samples within GSE42980 and 16 samples within GSE49121 were used to construct a consensus network and identify differentially expressed genes (DEGs).

Differentially expressed gene analyses
To investigate abdominal fat gene expression differences across two pairwise comparisons, 5462 and 6904 significant DEGs were identified in FL vs. LL and HG vs.
There were 3324 upregulated and 3580 downregulated DEGs in the HG vs.
Common DEGs across two pairwise comparisons were identified according to the following criteria: common DEGs were significant DEGs (FDR ≤ 0.05) within both pairwise comparisons and had the same expression trend across two pairwise comparisons, i.e., upregulated within both comparisons, and vice versa. A total of 2020 common DEGs were shared across the two pairwise comparisons, among which many lipid metabolism genes were significantly upregulated in the FL vs. LL and HG vs.

Gene coexpression network construction
After genes with low expression were removed, a total of 13,626 genes were used to construct a consensus gene coexpression network across the two datasets (GSE42980 and GSE49121). A soft threshold (power) of 14 was selected to construct a scale-free network within the two datasets (Fig. 1). Then, we separately constructed gene coexpression networks for GSE42980 and GSE49121 with the same parameters (power = 14, min-ModuleSize = 30, deepSplit = 2, and mergeCutHeight = 0.3). A total of 37 and 17 gene coexpression modules were detected within GSE42980 and GSE49121, respectively (Additional file 4: Figure S2, b-c). Finally, we constructed a consensus gene coexpression network across two datasets with these parameters (power = 14, minMo-duleSize = 30, deepSplit = 2, and mergeCutHeight = 0.3), and a total of 45 consensus modules were identified in this consensus network (Table 1; Additional file 4: Figure S2, a).
Clustering dendrograms of GSE42980 and GSE49121 revealed a degree of preservation of consensus modules across the two datasets ( Fig. 2a-b), and heatmaps showed similar results (Fig. 2 c, and f). The density value D = 0.66 reflected an overall preservation of consensus modules across the two networks (Fig. 2d). The adjacency heatmap indicated a preservation network across the two datasets, and each row and column corresponded to a consensus module (Fig. 2e).

Consensus modules correlated with abdominal fat deposition
To detect consensus modules associated with abdominal fat deposition, we calculated the correlation between each of the consensus modules and abdominal fat deposition across the two datasets. Eight consensus modules (dark green, green, medium purple 3, purple, saddle brown, sky blue, turquoise, and yellow) were significantly correlated with abdominal fat deposition across the two datasets, and the correlation of the yellow consensus module was the highest (r = 0.93, p = 9 × 10 −11 ) (Fig. 3). A total of 4816 genes were identified within eight consensus modules (Table 1). There were 13 consensus modules significantly correlated with fat deposition within GSE42980, and the correlation of the blue consensus module was the highest (r = 0.93, p = 3 × 10 − 10 ) (Additional file 5: Figure S3, a). A total of 27 consensus modules were significantly correlated with abdominal fat deposition within GSE49121, and the correlation of the steel blue consensus module was the highest (r = 0.99, p = 2 × 10 − 12 ) (Additional file 5: Figure S3, b).

Preservation of interested consensus modules
To verify the preservation of eight consensus modules across two datasets, we detected how many genes among the 2020 common DEGs were in the eight consensus modules. The results showed that a total of 1828 genes (90.5% of common DEGs) were detected in the eight consensus modules (Fig. 4).
To test whether the eight consensus modules were stable in other broiler lines, preservation of eight consensus modules was tested in another broiler abdominal fat RNA-Seq dataset (SRP058295). The Zsummary results showed that the turquoise, green, saddle brown, and yellow modules were strongly preserved (Zsummary > 10) and the dark green, purple, medium purple 3, and skyblue modules were moderately preserved (2 < Zsummary < 10) (Fig. 5).

Functional enrichment of interest consensus modules
Based on the KEGG and GO databases, 72 pathways and 749 GO terms were significantly enriched in eight consensus modules ( Table 2; Additional file 6: Table S3). Glycerolipid metabolism, glycerophospholipid metabolism, and adipocytokine signaling pathways were significantly enriched in the purple consensus module, and relevant genes of the three pathways, GPAT3, GPD2, PNPLA7, ETNPPL, and CD36, were identified. Steroid biosynthesis was significantly enriched in the medium purple 3 module, and relevant genes, DHCR7 and sterol-C5-desaturase (SC5D), were identified. Pyruvate metabolism, glyoxylate and dicarboxylate metabolism, and propanoate metabolism pathways were significantly enriched in the green module, and relevant genes of the three pathways, PDHA1, ACOX1, and SUCLA2, were identified. The PPAR signaling pathway was significantly enriched in the green consensus module, and relevant genes of this pathway, RXRG, ACOX1, and ACSL4, were identified. The fatty acid metabolism pathway was significantly enriched in the medium purple 3 and yellow consensus modules, and relevant genes of this pathway, ACSBG2, ACAT2, FASN, ACADS, PPT1, HACD1, FADS2, HACD2, and MCAT, were identified. Autophagy, lysosome, and mitophagy pathways were significantly enriched in the medium purple 3 and yellow consensus modules, and the relevant genes of the three pathways, ULK2, BECN1, MAP1 LC3A, MAP1 LC3B, ATG3, ATG4B, ATG12, and ATG9A, were identified in the consensus modules (Additional file 6: Table  S3).
To detect gene coexpression modules associated with abdominal fat deposition, we performed the enrichment of KEGG pathways and GO terms of modules within GSE42980 and GSE49121, respectively. There were 116 pathways and 1427 GO terms significantly enriched in the genes of modules significantly correlated with abdominal fat deposition within the GSE42980 dataset (Additional file 6: Table S3). A total of 79 pathways and 849 GO terms were significantly enriched in the genes of modules significantly correlated with the GSE49121 dataset (Additional file 6: Table S3).

Discussion
Broiler chickens have been selected to increase production performance, but these gains are also accompanied by excessive abdominal fat deposition [4]. The research results are different among different broiler lines according to published literature [6][7][8][9]12]. Therefore, to identify the genes associated with abdominal fat deposition across multiple broiler lines, we constructed a consensus gene coexpression network and identified eight consensus modules significantly correlated with abdominal fat deposition across the GSE42980 and GSE49121 datasets. Then, the preservation of eight consensus modules was verified in another broiler abdominal fat RNA-Seq dataset (SRP058295). Finally, functional enrichment analysis associated with the genes within eight consensus modules was performed.
Adipose tissue expansion is the result of adipogenesis and lipid droplet accumulation in adipocytes, and lipid droplets are an organelle for storing fat with a hydrophobic core of triacylglycerols (TAGs) and cholesterol esters [5]. Based on the pyruvate metabolism pathway in the KEGG database, PDHA1 and DLAT convert pyruvate into acetyl-CoA, which is the ingredient for de novo synthesis of fatty acids. FASN is a key enzyme catalyzing the synthesis of fatty acids, and the FASN concentration of a tissue determines the synthesis capacity of fatty acids in this tissue by a de novo pathway [13]. The expression level of ACSBG2 mRNA is significantly upregulated during abdominal fat-derived preadipocyte differentiation of chickens [14,15]. The fatty acid metabolism genes FADS2 and SCD are significantly upregulated in chickens with high abdominal fat content and associated with the PPAR signaling pathway [16]. In our current study, the pyruvate metabolism pathway and fatty acid metabolism pathway were significantly enriched in the green, yellow, and medium purple 3 consensus modules, and relevant genes of the two pathways, PDHA1, FASN, ACSBG2, FADS2 and SCD, were identified. The five genes were significantly upregulated in both the FL vs. LL and HG vs. LG, and DLAT was significantly upregulated in the FL compared with the LL.
The PNPLA protein family plays key roles in lipid droplet homeostasis, phospholipid metabolism, and triglyceride hydrolysis [17]. PNPLA7 regulates very low density lipoprotein (VLDL) secretion by modulating apolipoprotein E (ApoE) stability and is positively correlated with plasma TAG levels in mice [18]. The expression of GPAT3 is significantly upregulated during adipocyte differentiation in mice, and overexpression of GPAT3 results in increased TAG in mammalian cells [19]. As the predominant enzyme for TAG storage, DGAT2 catalyzes and accounts for nearly TAG synthesis [20]. DHCR7 and SC5D are pivotal enzymes for catalyzing cholesterol synthesis, and two impaired enzymes lead to a deficiency of cholesterol synthesis [21]. In the current study, glycerolipid metabolism, glycerophospholipid metabolism, and steroid biosynthesis pathways were significantly enriched in purple and medium purple 3 consensus modules, and PNPLA7, GPAT3, DHCR7, DGAT2, and SC5D were identified in the three lipid metabolism pathways. PNPLA7, GPAT3, and DHCR7 were significantly upregulated in both the FL vs. LL and HG vs.
LG, and DGAT2 and SC5D were significantly upregulated in the HG compared with the LG lines. Therefore, we inferred that the synthesis of fatty acids, TAG, and cholesterol in broilers with high abdominal fat deposition was greater than that in broilers with low abdominal fat deposition. Abdominal fat deposition is a complex biological process that is controlled by a series of critical genes in signaling pathways. In chickens, PPARG is a well-known critical gene that is reported to positively regulate adipogenesis, lipid metabolism, adipocyte differentiation, and abdominal fat deposition [22][23][24]. RXRG positively regulates the PPAR signaling pathway and promotes the expression of FABP3 and SCD [25]. Active immunization with CD36 results in decreased visceral fat of male broilers [26]. In the current study, the PPAR and adipocytokine signaling pathways were significantly enriched in the green and purple consensus modules, respectively, and relevant genes of the two pathways, PPARG, RXRG, FABP3, SCD, and CD36, were identified and significantly upregulated in both the FL vs. LL and HG vs. LG.
Autophagy is a process in which cytoplasmic macromolecules and organelles are degraded by lysosomes [27] and plays a crucial role in adipose tissue homoeostasis [28]. The biological functions of autophagy in adipose tissue have been studied in humans and model animals, but autophagy studies have rarely been reported in chicken abdominal fat. The autophagy process consists of six main steps: phagophore initiation, nucleation, and expansion; autophagosome formation, autophagosomelysosome fusion, and lysosomal degradation [29]. Autophagy is initiated with activation of the serine/ threonine kinase ULK1/2, which forms a complex with ATGs, and the ULK1/2 complex triggers phagophore nucleation by activating the PI3KC3 complex [30,31]. The LC3 precursor is processed by ATG4, generating LC3-I (soluble LC3 form), which is conjugated to membrane phosphatidylethanolamine (PE) by ATG3, ATG7, and the ATG12-ATG5-ATG16L complex [32]. Conjugated production is membrane-associated lipidated LC3-II that attracts components for phagophore elongation [32]. Phosphoethanolamine cytidylyltransferase (encoded by PCYT2) and ethanolamine phosphotransferase 1 (encoded by SELENOI) are the two critical enzymes for PE synthesis [33,34]. In the current study, autophagy, mitophagy, and lysosome pathways were significantly enriched in the medium purple 3 and yellow consensus modules. Relevant genes of the three pathways included LG. Therefore, we inferred that autophagy initiation might participate in broiler abdominal fat deposition.

Conclusions
In our current study, we predicted eight consensus modules across four broiler lines (FL and LL, HG and LG), the preservations of which were verified in another broiler line dataset (SRP058295). Pyruvate metabolism, fatty acid metabolism, glycerolipid metabolism, steroid biosynthesis, the PPAR signaling pathway, and the adipocytokine signaling pathway were significantly enriched in the eight consensus modules. Relevant genes of these lipid metabolism pathways included PDHA1, FASN, ACSBG2, FADS2, SCD, GPAT3, PNPLA7, DHCR7, PPARG, RXRG, FABP3, and CD36, all of which were identified in eight consensus modules and significantly upregulated in both the FL vs. LL and HG vs.
LG. Based on lipid metabolism pathways enriched in eight consensus modules and overexpression of numerous lipogenic genes, we hypothesize that more fatty acids, TAGs, and cholesterol might be synthesized in broilers with high abdominal fat than in broilers with low abdominal fat. Autophagy, mitophagy, and lysosome were significantly enriched in the eight consensus modules, and the relevant genes of the three pathways, SELE NOI, PCYT2, ULK2, BECN1, MAP1 LC3A, MAP1 LC3B, ATG3, ATG4B, ATG12, and ATG9A, were significantly upregulated in both FL vs. LL and HG vs.
LG. Therefore, we inferred that autophagy initiation might participate in broiler abdominal fat deposition. These genes identified in eight consensus modules can provide targets for further research into abdominal fat deposition in multiple broiler lines. To the best of our knowledge, this is the first report on abdominal fat deposition across multiple different broiler lines via a consensus module analysis method.

RNA-Seq data collection
To reduce the influence of nongenetic factors, such as analytical methods, versions of the genome and annotated files, experimental technologies (RNA-Seq or microarray), sex, diet condition, and weeks of age, we only collected abdominal fat RNA-Seq data of male broilers at 7 weeks of age with normal feed. In the current study, we downloaded GSE42980 [6], GSE49121 [7], and SRP058295 [8] from the GEO (Gene Expression Omnibus) and SRA (Sequence Read Archive) databases of NCBI (National Center for Biotechnology Information).

Read QC and mapping
To ensure high-quality reads for downstream analyses, quality control (QC) of the RNA-Seq reads was performed using fastp (0.19.7) with default parameters and workflow [35]. After filtering out low-quality reads, the high-quality reads were mapped to the chicken genome using HISAT2 (2.2.0) with default parameters and workflow [36]. The chicken reference genome sequence (GCF_000002315.6_GRCg6a) and chicken genome annotation files (GCF_000002315.6_GRCg6a_genomic.gtf.gz) were downloaded from the NCBI genome assembly website (ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/ 000/002/315/GCF_000002315.6_GRCg6a/). Then, the read count of each gene was calculated using feature-Counts (v2.0.1) with default parameters [37]. A summary of the RNA-Seq datasets, QC, and mapping is provided in Additional file 1: Table S1.

Hierarchical cluster analyses of samples
Low-expression genes within the GSE42980 and GSE49121 datasets were filtered using the filterByExpr function in the EdgeR package (3.24.3) [38,39]. A total Preservation of eight consensus modules tested in the SRP058295 dataset using Zsummary. Zsummary > 10 indicates strong evidence of preservation; 2 < Zsummary < 10 indicates weak to moderate evidence of preservation; Zsummary < 2 indicates no evidence of preservation (b)  [40]. To detect sample clusters within two datasets, we performed hierarchical cluster analyses of samples within GSE42980 and GSE49121, with read counts of 13,626 genes using the hclust function of R software (average method). The hierarchical cluster results of GSE42980 and GSE49121 are shown in Additional file 2: Figure S1.

Identification of differentially expressed genes
Within DESeq2, read counts from each library were normalized for differences in sequencing depth, used to estimate dispersions from the mean and to fit a linear model according to a negative binomial distribution for each pairwise comparison. Differentially expressed genes (DEGs) within the FL vs. LL (GSE42980) and HG vs.
LG (GSE49121) were identified by testing the two factor variables of high-and low-abdominal fat deposition and batch using DEseq2 (1.22.2) with default parameters [41]. The significance threshold for DEGs was set at a false discovery rate adjusted p-value (FDR) ≤ 0.05. The significant DEGs in each pairwise comparison are listed in Additional file 3: Table S2.

Gene coexpression network construction with WGCNA
Gene coexpression networks were constructed using the WGCNA (version 1.68) package in R [42]. The genes were selected for further analyses after passing the quality test by a "goodGene" function in the WGCNA package. The soft threshold (power) was screened to form a scale-free network in both dataset matrices using the "pickSoftThreshold" function in the WGCNA software package. Then, gene coexpression networks for GSE42980 and GSE49121 were separately constructed using the blockwiseModules function in the WGCNA package [10] with these parameters: power = 14, minMo-duleSize = 30, deepSplit = 2, and mergeCutHeight = 0.3. A consensus gene coexpression network across two datasets was constructed using the blockwiseConsen-susModules function in the WGCNA package [10] with these parameters (power = 14, minModuleSize = 30, deepSplit = 2, and mergeCutHeight = 0.3). Consensus modules shared by two gene coexpression networks were detected using the dynamic tree cut method [10,43].

Correlation analysis between consensus modules and abdominal fat deposition
The module eigengene is defined as the first principal component of a gene module [10]. To identify modules of interest, the correlation between each module eigengene of consensus modules and abdominal fat deposition traits was calculated using the correlation analysis function in the WGCNA package [10], and significant consensus modules were identified with the threshold (p-value ≤0.05).

Preservation test of consensus modules
To test the preservation of consensus modules significantly correlated with abdominal fat deposition in other broiler lines, abdominal fat RNA-Seq data of male broilers (SRP058295) were downloaded to construct a test network for preservation analysis. The preservation analysis of consensus modules was performed using the modulePreservation function in the WGCNA package with the parameters nPermutations = 200, randomSeed = 1, quickCor = 0, verbose = 3 [11]. The medianRank was used to compare relative preservation among consensus modules, and a module with higher median rank exhibits weaker observed preservation statistics than a module with a lower median rank [11]. The thresholds of Zsummary were set as follows: Zsummary > 10 means strong evidence of preservation; 2 < Zsummary < 10 means weak to moderate evidence of preservation; Zsummary < 2 means no evidence of preservation [11].

Function enrichment of consensus modules
Based on KEGG (Kyoto Encyclopedia of Genes and Genomes) and GO (Gene Ontology) databases, the pathways and GO terms associated with the genes within significant consensus modules were predicted using the clusterprofiler software package (3.10.1) [44]. The significance threshold for the pathways and GO terms associated with the consensus modules was set at a p-value ≤0.05.
Ethics approval and consent to participate Not applicable.

Consent for publication
Not applicable.