Combined Transcriptome and Metabolome Analysis of Chicken Follicles in China Tengchong Snow Chicken Follicle Selection

1. Introduction

Egg-laying performance is an important economic index in the chicken breeding process and has attracted much attention from researchers. The Tengchong snow chicken is a famous domestic chicken breed in China. It is a very valuable chicken resource. It has the advantages of black flesh, a fresh taste, and beautiful feathers. However, the egg production of Tengchong snow chickens is low, and the egg-laying cycle of some chickens is shortened [1]. Therefore, improving the egg production of Tengchong snow chickens has become a matter of concern for the chicken industry. Improving the genetic potential of egg production is an important strategy for increasing egg production.

Egg-laying performance depends on the development and selection of ovarian follicles in chickens [2]. The development and selection of follicles is an extremely complex process that eventually atresia or develops until ovulation [3]. Follicles in chickens start developing with the recruitment of the original follicle, which continues to grow, undergoes selection, and finally matures and ovulates [4]. After hens become sexually mature, their ovaries get filled with many follicles at different stages of development, which are usually arranged according to the strict hierarchy of follicles and are generally divided into the PHF and hierarchical follicle (HF) stages [5]. According to the volume (from small to large) and yolk color, the former can be divided into the small white follicle (SWF, diameter 1–4 mm) stage, large white follicle (LWF, diameter 4–6 mm) stage, and small yellow follicle (SYF, diameter 6–8 mm) stage. The latter is divided into F5 to F1 stages according to the volume from small to large; the follicles in the F1 stage are about to ovulate, followed by the F2 stage, and so on. Each stage is a tightly controlled process [6]. The follicle reserve at each stage may affect the ovulation cycle, and ultimately, ovulation efficiency, especially in the SYF stage [7]. To maintain sustained ovulation, 1–2 SYFs are selected from the HF stage every day to enter the rank follicle stage, and then they start growing rapidly and eventually differentiate to become the dominant follicle, a process also known as follicle selection [8]. The greater the number of 6–8 mm SYFs selected to enter the HF, the longer the laying cycle and the better the egg-laying performance. Thus, the development of PHF, especially SYF, directly affects the recruitment of dominant follicles; therefore, the development of PHF and their regulatory mechanism need to be studied to improve egg-laying performance.

In this study, we focused on the LWF and SYF stages. These are the critical stages of follicle selection and represent the follicles before and after selection. We aimed to identify regulatory genes and metabolites involved in the development and selection of follicles in chickens. Our findings may provide new insights into the morphological changes and molecular mechanisms underlying follicle selection.

2. Results

2.1. Histological Observation of Follicles

One or two layers of GCs were found in the LWF (Figure 1 A) and two layers of GCs were found in the SYF (Figure 1B). During follicle development, the layers of SYF GCs (9.92–1.20 μm; Figure 1B) were significantly thicker than those of LWF GCs (7.92–1.30 μm; Figure 1A). The area of single GCs in the SYF was significantly greater than that in the LWF, and the GCs in the LWF were more densely arranged (Figure 1 A and B). Figure 1. This is a figure. Schemes follow the same formatting.

2.2. Transcriptome Analysis

2.2.1. Data Summary of the Transcriptome

• Eight cDNA libraries were constructed using the samples obtained from the SYF and LWF groups. A total of 379.34 million raw reads were obtained. After removing adaptors and low-quality reads, 374.79 million clean reads were obtained. The average number of clean reads was 97.92% for Q20 and 94.21% for Q30. The average GC content of the clean reads was 49.05%. Eight cDNA library sequences were mapped to the genome, and the comparison range was 88.25–90.52% (Table 1). FASTQ files were submitted to the NCBI database, and when the article is accepted and the data is released immediately.

2.2.2. Analysis of DEGs

A total of 258 DEGs were identified in SYF vs. LWF, of which 172 were upregulated and 86 were downregulated (Figure 2A, Table S2). In the heat map, the color difference reflects the information of the DEGs. The samples from the same group were clustered together to visualize the gene expression patterns of each group (Figure 2B).

2.2.3. KEGG Analysis of DEGs

The functional enrichment of the DEGs was performed using KEGG. We obtained 17 KEGG signaling pathways in SYF vs. LWF, which mainly included gap junctions, the calcium signaling pathway, and neuroactive ligand-receptor interactions (Figures 2C, Table S3).

2.2.4. Verification of DEGs by RT-qPCR

We randomly selected the DEGs (GRXCR2, AVPR1B, BMP15, DMRT2, STRA6, and GP9) for RT-qPCR validation. The results showed that the relative levels of expression of three genes (DMRT2, STRA6, and GP) were greater in SYF than in LWF. The relative expression levels of three genes (GRXCR2, AVPR1B, and BMP15) were lower in SYF than in LWF (Figure 2D). The relative mRNA levels of the selected genes were consistent with the results of the transcriptome analysis.

2.3. Metabolomics Analysis

2.3.1. Analysis of PCA and OPLS-DA

• Principal component analysis (PCA) is an unsupervised multivariate statistical analysis method used to determine overall metabolic differences between groups as well as changes in samples within each group. The PCA score plot results showed that the overall trend of metabolite distribution in the SYF and LWF samples was separated in the negative model (Figure 3A) and the positive model (Figure 3D). Additionally, by plotting the OPLS-DA score, a clear separation of metabolites between the SYF and LWF groups was found in the negative model (Figure 3B) and the positive model (Figure 3E). Displacement tests were conducted to determine the accuracy of the OPLS-DA model (Figure 3C and F).

2.3.2. Analysis of DMs

In the negative mode, 77 significant DMs were identified, including 13 upregulated DMs and 64 downregulated DMs (Figure 4A and Table S4). In the positive mode, 52 significant DMs were identified, including 23 upregulated DMs and 29 downregulated DMs (Figure 4D and Table S5). In the heat map, the color difference reflects the information of the DMs. The samples from the same group clustered together, and the heat map visually reflected the metabolomic expression patterns of the three groups (Figure 4B and E).

2.3.3. KEGG Analysis of DMs

The functional enrichment of DMs in the SYF and LWF groups was performed using KEGG. In total, seven and two significantly enriched pathways were identified in the negative model (Figure 4C and Table S6) and the positive model (Figure 4F, Table S7), respectively. The pathways mainly included glutathione metabolism, alpha-linolenic acid metabolism, tryptophan metabolism, linoleic acid metabolism, etc.

2.4. Integrated Analysis of Transcriptomics and Metabolomics

The KEGG analysis of the DEGs and DMs showed that the metabolome and transcriptome shared seven KEGG pathways (Figure 5A and B). Four pathways associated with metabolism, including glutathione metabolism, alpha-linolenic acid metabolism, linoleic acid metabolism, and pyrimidine metabolism, were enriched (Figure 5B). We examined the correlations between DMs and DEGs related to glutathione metabolism, alpha-linolenic acid metabolism, linoleic acid metabolism, and pyrimidine metabolism. A Pearson correlation analysis of 10 metabolites and 14 genes showed that the 13 genes were significantly positively or negatively correlated with one or more of the 10 metabolites (Figure 5C and D).

Based on the above results, we identified a network that may regulate follicular development in Tengchong snow chickens (Figure 6), which involves mainly glutathione metabolism pathways associated with SDMs and DEGs. However, the regulatory network of follicular development in chickens is very complex, and only a small part of this network is understood. Therefore, further studies are needed to investigate and confirm these findings.

3. Discussion

As the development of SYF directly affects the recruitment of dominant follicles, it is important to study the development of SYF and its regulatory mechanism for improving egg-laying production performance. As a key component of ovarian follicles, GCs proliferate rapidly during follicular development to provide nutrients to these follicles [[10[M1] ,11]. Second, GCs proliferate rapidly and differentiate into grade GCs that can secrete various factors to regulate the maturation of oocytes and the production of steroid hormones [12]. We made paraffin sections of chicken follicles and observed their morphological characteristics. We found that the GC layer of SYF was significantly thicker than that of LWF during the selection of follicles. The thickness of the GC layer increased due to the proliferation of GCs, which promotes the transition from pre-layered follicles to layered follicles [13,14]. With the development of the GC layer, the area of a single GC increased, and the layer showed a loose arrangement. Therefore, the development of SYF and its regulatory mechanism are highly important for improving egg-laying performance. However, the specific molecular regulatory mechanism involved needs to be elucidated.

The KEGG pathway analysis revealed that the DEGs identified in the SYF and LWF groups were associated with 17 significant signaling pathways, such as the calcium signaling pathway, gap junction, neuroactive ligand-receptor interaction, etc. The calcium signaling pathway that initiates development in many mammalian and non-mammalian follicles has been extensively studied [15,16,17]. The calcium signaling pathway stimulates the maturation of the oocyte [18]. This pathway involves the GRIN2A, EGF, and FGF8 genes. Ovarian steroid hormones can increase GRIN2A expression, thus affecting follicle development [19]. Additionally, GRIN2A may be used as a biomarker to predict oocyte maturation in CCs [20]. EGF directly affects follicle development [21]. FGF8 directly affects the number of bovine ovarian sinus follicles [22]. Gap junctions affect the development of ovarian folliculation [23]. This pathway involves the TUBA1C, ADCY7, and EGF genes. No study has reported the connection between the TUBA1C and ADCY7 genes and follicle development. As EGF is closely related to the development of follicles, the EGF gene may play an important role in ovarian development. The neuroactive ligand-receptor interaction pathway also plays an important role in the development of follicles. Studies have indicated that neuroactive ligand-receptor interactions are related to ovarian follicle development and egg production in poultry [24,25,26]. This pathway involves the P2RY8, P2RX2, and GRIN2A genes. No study has reported the relationship between the P2RY8, P2RX2, and GRIN2A genes and follicle development. These unreported genes need to be further investigated.

Metabolomics is a relatively new field that has emerged in response to genetics and proteomics. It can provide information on the physiological state of organisms by monitoring changes in endogenous metabolites. Metabolomics is extensively used in various fields to study a wide range of animals. In this study, the metabolome analysis revealed 129 significant DMs, including 36 upregulated and 93 downregulated DMs. These DMs may also play important roles in the development of chicken follicles.

The DMs were mapped to nine KEGG pathways in the SYF and LWF groups, which included glutathione metabolism, alpha-linolenic acid metabolism, pyrimidine metabolism, etc. Glutathione is an antioxidant and electrophilic toxicant that prevents the destruction of ovarian follicles and oocytes by toxicants [27,28]. Glutathione helps protect mature ovarian follicles from apoptotic stimuli, including PAHs [29,30]. Additionally, glutathione deficiency affects germ cell apoptosis in cultured embryonic mouse ovaries [31]. These findings suggest that glutathione metabolism may significantly affect follicle development in chickens. In cattle, alpha-linolenic acid significantly affects oocyte development [32]. High alpha-linolenic acid intake can promote ovulation in rats [33]. Therefore, alpha-linolenic acid metabolism also plays an important role in follicular development in chickens. Pyrimidine metabolism is closely related to follicle development in Chinese hamsters [34]. Our results suggested that these pathways were enriched, indicating that they regulate follicle development in chickens.

The Venn diagram plotted in this study indicated that the DMs and DEGs that shared four pathways were associated with glutathione metabolism. Our findings also showed that glutathione metabolism plays a role in regulating follicle development.

4. Materials and Methods

4.1. Animals and Sample Collection

Tengchong snow chickens obtained from Baoshan, Yunnan Province, China, were raised at the experimental chicken farm of Kunming University. All chickens were kept isolated in the same environment under a 10-h/14-h light/dark cycle. Egg-laying in chickens reaches peak performance when they are 26 weeks old. The hens were euthanized, and follicles were collected at 26 weeks of age. The follicles used for RNA sequencing and RT-qPCR were immediately frozen in liquid nitrogen and then stored at -80 °C.

4.2. Histological Observation of Follicles

At the age of 26 weeks, 24 hens (n = 12 from each group) were slaughtered, and their ovaries were collected and dissected. The follicles were fixed in a 4% (volume/volume) buffered formaldehyde solution for 24 h, routinely dried, and embedded in paraffin wax. The sections were stained with hematoxylin and eosin. All sections were observed under a microscope.

4.3. Transcriptome Analysis

4.3.1. RNA Extraction, cDNA Library Construction, and Sequencing

Total RNA was extracted from LWFs and SYFs of chickens (n = 4, each group) using TRIzol reagent (Invitrogen, CA, USA) following the manufacturer's instructions. Quality control of RNA was performed using an Agilent 2100 bioanalyzer (Agilent Technologies, CA, USA) to detect the integrity and total volume of RNA. A cDNA library was constructed, high-throughput sequencing was performed on the Illumina HiSeq 5000 platform of Beijing Novo Technology Co., Ltd., and the paired-end 150 bp sequencing platform was used for sequencing. By capturing the fluorescence signal and converting the light signal into a sequencing peak using computer software, the sequencer can obtain the sequence information of the fragment to be tested. Eight biological samples were used for RNA-seq.

4.3.2. Bioinformatics Analyses

The sequence fragments were converted into sequence data (reads) by CASAVA (Version 1.8) base recognition. The RNA-seq data were pretreated with fastq (version 0.19.7) to remove adapters and low-quality sequences. The clean reads were mapped to the chicken reference genome using the HISAT2 software (version 2.0.5) (ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCF_000002315.4_Gallus_gallus-5.0_genomic.gff.gz). The result obtained by the Hisat2 software was a sam format file, and we converted it into a bam format file using samtools (Version 1.5). We visually browsed the bam file with the Integrative Genomics Viewer (IGV) browser. The featureCounts tool in subreads was used to quantify gene expression. The Cufflinks (version 2.2.1) software was used to normalize gene expression levels based on the number of transcript fragments per kilobase per million mapped reads (FPKM). The number of reads per gene map was counted using the HTSeq (Version 0.6.1) software. DEGs were identified using the DESeq2 (version 1.32.0) software. The screening criteria for DEGs were a |log2-fold fold change| ≥1 and padj < 0.05. The DEGs of all comparison groups were collected and combined into differential gene sets. GO functional enrichment analysis and KEGG pathway enrichment analysis were performed using clusterProfiler.

4.3.3. Quantitative Reverse Transcription PCR (qRT-PCR)

To determine the accuracy and repeatability of the RNA-seq results, six genes were selected for RT-qPCR to measure their expression levels in different follicles. Primers for the amplification of candidate genes, including NPC intracellular cholesterol transporter 2 (NPC2), glutaredoxin and cysteine rich domain containing 2 (GRXCR2), arginine vasopressin receptor 1B (AVPR1B), bone morphogenetic protein 15 (BMP15), doublesex and mab-3 related transcription factor 2 (DMRT2), stimulated by retinoic acid gene 6 (STRA6), and glycoprotein IX platelet (GP9), were designed using Prime-Blast tools from the NCBI database. These primers are listed in Supplemental Table S1.

Total RNA was extracted from the follicle tissue of chickens in each group (n = 6) using a method identical to that used for transcriptome analysis, and cDNA was synthesized using a reverse transcription kit (TaKaRa, Dalian, China). RT-qPCR was performed on a CFX96TM Real-Time System (Bio-Rad, CA, USA) using a SYBR premix Ex TaqTM II kit (Tiangen, Beijing, China). The 2-ΔΔCT method was used to determine gene expression levels, and beta-actin was used as the endogenous reference gene for normalization.

4.4. Metabolomics Analysis

4.4.1. Metabolite Extraction

Metabolites were extracted from chicken SYF and LWF tissues. Four replicates were used for each group, and eight libraries were constructed. The tissue samples (100 mg) were mixed with liquid nitrogen and placed in an EP tube, and 500 μL of 80% methanol aqueous solution was added. The mixture was incubated in an ice bath for 5 min and then centrifuged at 15,000 ×g at 4 °C for 20 min. Then, a certain amount of the supernatant was collected and diluted with ultrapure water until the methanol content decreased to 53%. The mixture was centrifuged at 15,000 ×g at 4 °C for 20 min, the supernatant was collected, and LC-MS/MS analysis was performed.

4.4.2. LC-MS/MS Conditions

The metabolites were separated using a Vanquish UHPLC (Thermo Fisher, MA, USA) Hypesil Gold column (Thermo Fisher, MA, USA). The column temperature was 40 °C, and the flow rate was 0.2 mL/min. The mobile phase consisted of buffer A (0.1% formic acid) and buffer B (methanol). Throughout the analysis, the samples (2 µL) were placed in an automatic injector at 4 °C. QC samples were inserted into the sample queue to monitor and evaluate the stability of the system and the reliability of the test data. The samples were analyzed by a high-resolution mass spectrometer (Q Exactive™ HF-X; Thermo Fisher, MA, USA) in both positive and negative ion modes of electrospray ionization (ESI).

4.4.3. Metabolomics Data Analysis

The KEGG database (https://www.genome.jp/kegg/pathway.html), HMDB (https://hmdb.ca/metabolites), and LIPIDMaps database (http://www.lipidmaps.org/) were used to identify the metabolites.

For multivariate statistical analysis, the metabolomics data processing software metaX [9] was used to transform the data, and then principal component analysis (PCA) and partial least square discriminant analysis (PLS-DA) were performed to obtain the VIP value of each metabolite. In the univariate analysis, the statistical significance (P value) of each metabolite between the two groups was calculated by conducting t-tests, and the fold change (FC value) of the metabolite between the two groups was calculated. The default criteria for differential metabolite screening were VIP > 1, P < 0.05, and FC ≥ 2 or FC ≤ 0.5. The volcano map was drawn using the R package (version 3.4.3) ggplot 2, and three parameters of the VIP value, log2 (fold change) and -log10 (P value), of the metabolites were integrated to screen the metabolites of interest.

A bubble map was drawn using ggplot2, and the KEGG database was used to study the functions and metabolic pathways of the metabolites. When x/n>y/n, the metabolic pathway was considered to be enriched. The metabolic pathway was considered to be significantly enriched when the P-value was <0.05.

4.5. Conjoint Analysis of the Metabolome and Transcriptome

The metabolome and transcriptome were analyzed together. Common significant signaling pathways between the transcriptome and metabolome were identified and displayed using the online Wien tool (https://www.omicstudio.cn/tool). In order to investigate the correlation between DEGs and DMs, a combined analysis of the selected DEGs and DMs from both the metabolome and transcriptome was performed, followed by the construction of a transcriptional-metabolite network with a Pearson correlation coefficient exceeding 0.7 (https://www.omicstudio.cn/tool).

4.6. Statistical Analyses

All data were statistically analyzed using SPSS 22.0. The data were presented as the mean ±standard error. The differences in parameters between the two groups of samples were determined by conducting Student’s t-tests. All differences among and between groups were considered to be statistically significant at P < 0.05; one asterisk (*) indicates P < 0.05, and two asterisks (**) indicate P < 0.01.

5. Conclusions

To summarize, in this study, we analyzed the genetic landscape of ovarian follicle development and selection in chickens. In total, 258 DEGs were identified in SYF vs. LWF. Functional annotation analysis revealed that gap junctions, the calcium signaling pathway, neuroactive ligand-receptor interactions, etc., were involved in chicken follicle selection and development. In total, 129 significant DMs were identified in SYF vs. LWF. The results of a functional annotation analysis showed that glutathione metabolism, alpha-linolenic acid metabolism, and pyrimidine metabolism were involved in chicken follicle selection and development. By performing an association analysis between transcriptomics and metabolomics, four metabolically related pathways, including glutathione metabolism, alpha-linolenic acid metabolism, linoleic acid metabolism, and pyrimidine metabolism, were identified. Glutathione metabolism is a critical signaling pathway that may affect ovarian follicle development to regulate the SYF reserve in chickens. However, the genetic mechanisms underlying follicle development and selection in chickens are very complex, and the functions and regulatory networks associated with these important metabolites and genes need to be further studied. The results of this study help us understand the molecular processes of follicle development and selection in chickens and provide a theoretical basis for studying the molecular regulatory mechanisms underlying the development of chicken follicles.