Comparative Transcriptomic Profiling in Ovarian Tissues of Lohmann Hens and Chengkou Mountain Chicken

Background: As a crucial economic characteristic and a major indicator of reproductive performance in layers, egg production is controlled by a series of complex regulatory heredity basis. In particular, the interacting regulatory function between noncoding RNAs (ncRNAs) and coding RNA plays important roles in regulating laying performance. Methods: In this study, the RNA sequencing (RNA-seq) of ovarian tissues from Lohmann hens (n = 3) and Chengkou Mountain chicken (n = 3) under the laying peak period was performed to identify RNA transcriptional differences among different laying-performance populations. Results: Results showed that the expression level of 303 mRNAs, 68 long ncRNAs (lncRNAs), 533 circular RNAs (circRNAs), and 79 microRNAs (miRNAs) was significantly different among the groups. Functional enrichment analysis of these differentially expressed (DE) mRNAs revealed that the laying process was implicated in numerous significantly enriched pathways (p $<$ 0.05), such as the neuroactive ligand-receptor interaction, steroid hormone biosynthesis, and calcium-signaling pathway. Furthermore, the lncRNA/circRNA–miRNA–mRNA regulatory networks related to the regulation of laying performance were constructed. Some randomly selective DE RNAs were verified by Real Time Quantitative (RT-qRCR), indicating that the bioinformatics analysis results of RNA-seq data were credible. Conclusions: This study could increase our understanding of the heredity basis of transcriptome in the laying performance of chicken.

Keywords

egg production

noncoding RNA

RNA-seq

chengkou mountain chicken

1. Introduction

Egg production performance as a key economic characteristic in the poultry industry is reflected in the productive capacity of layers caused by population growth and a meet demand of animal products [1]. Egg production performance is regulated by many factors. In particular, the development of the ovary and follicle is the most influential factor on poultry production performance [2]. As an endocrine gland, the ovary is also involved in the production and secretion of related reproductive hormones [3]. Therefore, ovarian function has a direct impact on laying performance.

To date, transcriptional analysis of chicken ovarian tissues has been extensively researched, and a series of candidate genes, such as zona pellucida glycoproteins 2 (ZP2), Wnt family member 4 (WNT4), and anti-Mullerian hormone (AMH), has been identified to involve in ovarian development and its function maintenance [4]. In addition, functional noncoding RNAs (ncRNAs) can affect follicular development through various molecular regulatory pathways in chicken [5]. In particular, several microRNAs (miRNAs) have been confirmed to regulate hormone synthesis and follicular development, such as miR-26a-5p, miR-205b, and miR-23b-3p [6, 7]. Furthermore, several circular RNAs (circRNAs) are implicated in the development of mammalian follicles. For example, a previous study has shown that the host gene (LPAR3) of the circRNA (novel_circ_0001794) in the oviduct between follicular and luteal phases affects sheep reproduction via the Rap1 signaling and PI3K-Akt signaling pathways [8]. Long non-coding RNAs (lncRNAs) are also essential for various biological functions, including reproduction. However, few studies have been conducted on functional lncRNAs regulating ovarian growth and development in chicken.

Subsequently, some lncRNAs or circRNAs have been reported to play important roles in a series of physiological processes in livestock and poultry by influencing miRNAs. For example, circEML1 can release IGF2BP3 expression by sponging gga-miR-449a to facilitate steroidogenesis in chicken follicular granulosa cells [9]. Thus, a new theory has been proposed, that is, competing endogenous RNAs (ceRNAs) can competitively bind miRNA to regulate gene expression through miRNA response elements [10].

For a long time, chicken, as a widely raised animal, can provide high-quality protein products to human. Its breeding performance directly affects the economic benefits of chicken enterprises [11]. The complex topography and diverse ecological environment create a large number of indigenous chicken breeds worldwide. In particular, numerous excellent local chicken breeds in southwest China have perfect ecology adaptability to habitats, and they play an irreplaceable role in local agricultural development. However, most of them have lower egg and meat production because they lacked long-term artificial selection in comparison with commercial chicken [12]. Chengkou Mountain chicken (CK), a distinctive native breed of chicken in Chongqing, has features common to mountain chicken in the southwestern hilly region of China with rough feeding resistance, high adaptability, and superior meat quality. However, the short peak period and low egg production of CK limit its industrial economic benefits.

The present study identified the RNA expression profiles (mRNA, lncRNA, circRNA, and miRNA) in ovarian tissues under the laying peak period of CK and Lohmann hens (LM). This study aimed to understand the regulatory heredity basis of different ovarian oviposition abilities and provide a new breeding reference for Chinese local chicken.

2. Materials and Methods

2.1 Ovarian Tissue Collection

The ovarian tissues of six healthy layer individuals with ecologically similar habitats were collected during the peak laying period, including three CK chickens from Chengkou Mountain Chicken Genetic Resources Research Institute (108°54 ${{}^{\prime}}$ 48 ${{}^{\prime\prime}}$ , 31°50 ${{}^{\prime}}$ 57 ${{}^{\prime\prime}}$ , Chongqing, China) and three LM hens from Chongqing Fudi Poultry farm (106°28 ${{}^{\prime}}$ 23 ${{}^{\prime\prime}}$ , 29°35 ${{}^{\prime}}$ 37 ${{}^{\prime\prime}}$ , Chongqing, China).

All animals were euthanized with xylazine hydrochloride (8 mg/kg). Then, the ovaries (including the entire ovary of the large and small yellow follicles) of each individual were immediately collected and rapidly frozen in liquid nitrogen. Upon approval from Southwest University’s Ethics Committee for Animal Experiments, the experiments were conducted (IACUC-20211020-08).

2.2 Extraction of Total RNA, Quality Check, and Library Preparation

A Trizol reagent kit (Invitrogen, Carlsbad, CA, USA) was used for the extraction of total RNA in ovarian tissues, and DNase I was used (TaKaRa, Dalian, China) for potential genomic DNA contamination. RNA quality estimation was performed using an Agilent 2100 bioanalyzer (Agilent Technologies, Santa Clara, CA, USA). The NEBNext Ultra Directional RNA Library Prep Kit (NEB, Ipswich, MA, USA) was used for Illumina sequencing of coding RNA and lncRNA libraries. Using the NEBNext Multiplex small RNA library preparation set (NEB, Ipswich, MA, USA), the miRNA library was constructed on the basis of identical samples for Illumina sequencing. The Illumina HiSeq ${}^{\text{TM}}$ 4000 (Gene Denovo Biotechnology Co., Guangzhou, China) was used to sequence all libraries at high throughput.

2.3 Filtering Data, Quality-Control Transcripts, and Annotating Genes

Raw data were processed with FASTQ (v0.11.8, http://www.bioinformatics.babraham.ac.uk/projects/fastqc/) software to remove low-quality adapter reads [13]. Bowtie2 (v2.3.2) [14] and TopHat2 (v2.1.1) [15] were used to align reads with GCF_000002315.5, a chicken reference genome. Cufflinks (v2.1.1) [16] was used to assemble transcripts and estimate abundances. BLAST (v2.2.25, https://www.ncbi.nlm.nih.gov/Class/BLAST/) was used for mRNA NR and KEGG annotation (E-value $<$ 10 ${}^{-5}$ ).

Using default parameters, three software were used to assess the protein-coding potential of novel transcripts: CPC [17], CNCI [18], and Pfam [19]. LncRNAs were selected from the junction of non-protein-coding potential results. Using GenBank (Release 209.0) and Rfam (Release 11.0) databases, rRNA, tRNA, snRNA, scRNA, and snoRNA were identified and removed from the clean tags. Clean tags were filtered from the miRbase database (Release 22.0) by Bowtie to identify miRNAs. Using Mireap (v2.0.5, http://sourceforge.net/projects/mireap/) and comparing the unannotated tags to the reference genome, novel miRNA candidates were predicted. The intersected predicted target gene of miRNAs was determined by MiRanda (v3.3) [20] and TargetScan (Release 7.2, http://www.targetscan.org/). MRNA and miRNA expression correlation was evaluated using the Spearman rank correlation coefficient. The miRNA–mRNA network was visualized using Cytoscape (v3.9.1, http://www.cytoscape.org/). CircRNAs and their host genes were predicted using Findcirc (v2.0) software [21].

2.4 Identification of Significantly and Differentially Expressed (DE) RNAs and Functional Enrichment of Candidate Genes

By calculating fragments per kilobase of transcript per million mapped reads, mRNA and lncRNA transcript abundance was quantified. Transcripts per million were used to calculate miRNA expression. Counted reads per million indicated circRNA expression.

Differentially expressed transcripts among groups were identified using edgeR package in accordance with the following criteria: DE mRNAs (DEMs) and DE lncRNAs (DELs) with $|{}$ log ${}_{2}$ (Fold Change) $|{}$ $>$ 1 and threshold false discovery rate $<$ 0.05; DE circRNAs (DECs) and DE miRNAs with $|{}$ log ${}_{2}$ (Fold Change) $|{}$ $>$ 1 and p $<$ 0.05. GO (Gene Ontology) and KEGG (Kyoto Encyclopedia of Genes and Genomes) enrichment analyses were performed using KOBAS 3.0 (KEGG Orthology Based Annotation System, http://bioinfo.org/kobas), and p $<$ 0.05 was determined to be significant.

2.5 CeRNA Network Construction and RT-qPCR Validation

lncRNA/circRNA–miRNA–mRNA ceRNA networks were constructed on the basis of all DEMs, DELs, DECs, and DE miRNAs. Their genetic pattern was visualized using Cytoscape.

The validation of RT-PCR was performed using the same RNA samples as RNA sequencing (RNA-seq). Takara’s PrimeScript ${}^{\text{TM}}$ RT reagent kit (TaKaRa, Dalian, China) was used to remove genomic DNA from RT-PCR reactions for further verification of mRNA, lncRNA, and circRNA. In validating miRNAs, the Mir-X miRNA First-Strand Synthesis Kit (TaKaRa, Dalian, China) was used. All reagents used for RT-PCR in 96-well plates were detected using the Bio-Rad iQ5 Real-time PCR Detection System (Applied Biosystems, Foster City, CA, USA) and SYBR Green qPCR Mix (Life Technologies, Carlsbad, CA, USA).

The mixture consisted of 20.0 µL, which included 10 µL of TB Green Premix Ex Taq II (2 $\times{}$ ), 2.0 µL of cDNA, 7.0 µL of H ${}_{2}$ O, and 0.5 µL of forward and reverse primers (10 µM). In PCR amplification, 10 min of predenaturation at 95 °C is followed by 40 runs at 95 °C for 15 s and at 60 °C for 45 s. The melting curves were analyzed after amplifications, and the experiments were performed in triplicate. Negative controls without a cDNA template were included in this experiment. The primer information of all candidate genes and housekeeping genes (GAPDH and U6) is shown in Supplementary Table 1. Relative quantitation among three biological replicates of samples was calculated using the 2 ${}^{-\Delta{}\Delta{}\text{Ct}}$ method. Data were analyzed using two-tailed Student’s t test and presented as means $\pm{}$ standard deviation.

3. Results

3.1 RNA-Seq Results

All mRNA and lncRNA libraries obtained 70.84 GB of raw data, and 70.01 GB of high quality (HQ) clean reads were retained after being filtered (Table 1). The proportion of Q30 ranged from 92.34% to 93.27%, and approximately 93.68% of the reads could be mapped to the chicken reference genome. Moreover, small RNA-seq data resulted in 80,795,346 raw reads, and 80,034,358 clean reads were achieved after being filtered (Table 2). Over 72% of mapped reads were consistent with the reference genome. An upload of this study’s data was made to the SRA database (PRJNA838653 and PRJNA838418).

Table 1.The sequencing data of lncRNA and mRNA libraries.

Sample	Raw reads (bp)	High quality (HQ) clean reads (bp)	Q20	Q30	Mapping ratio (%)
LM_1	11,273,313,000	11,137,064,634	97.12%	92.34%	93.48
LM_2	11,008,686,600	10,882,438,725	97.54%	93.27%	93.88
LM_3	12,079,826,100	11,923,677,723	97.44%	93.18%	93.55
CK_1	11,790,882,000	11,674,353,315	97.22%	92.54%	93.80
CK_2	11,848,987,500	11,713,521,638	97.44%	93.09%	93.83
CK_3	12,844,688,700	12,686,296,548	97.11%	92.37%	93.59

Table 2.The sequencing data of miRNA libraries.

Sample	Raw reads (bp)	HQ clean reads (bp)	3’adapter_null	5’adapter	Match_abundance (%)
LM_1	15,992,789	15,820,964	0.1580%	0.0515%	74.00
LM_2	11,878,104	11,780,578	0.1061%	0.0328%	72.08
LM_3	10,655,308	10,546,608	0.1318%	0.0236%	73.31
CK_1	10,421,875	10,336,230	0.1184%	0.0225%	72.81
CK_2	15,005,429	14,879,262	0.2307%	0.0512%	75.14
CK_3	16,841,841	16,670,716	0.2478%	0.0271%	74.92

3.2 DE Patterns of mRNA, lncRNAs, circRNAs, and miRNAs between the LM and CK Groups

With high-throughput sequencing, 303 significant DE mRNA transcripts (55 upregulated and 248 downregulated) were identified between the two groups, including 287 known genes (Supplementary Table 2). Second, 68 (33 upregulated and 35 downregulated) lncRNA transcripts (Supplementary Table 3) and 533 (248 upregulated and 285 downregulated) circRNA transcripts (Supplementary Table 4) were differentially expressed within CK compared with LM. Third, the results of small RNA-seq showed that 79 miRNAs were significantly and differentially expressed between the two groups (Supplementary Table 5). The expression level of DE RNAs between the two chicken groups and the statistical significance of the differences are shown by volcano plots (Fig. 1A–D).

Fig. 1.

Volcano plots of transcripts in ovarian tissue. Volcano plots of differentially expressed (DE) mRNAs (A), lncRNA (B), circRNA (C), and miRNAs (D) in ovaries between CK and LM groups. Red, orange, and blue dots represent significantly downregulated, upregulated, and expression with no significance.

3.3 RT-qPCR Validation

A total of 16 DE RNA transcripts were randomly selected to verify the accuracy of the sequencing results using RT-qPCR. Fig. 2 shows the RT-qPCR results of 16 DE RNAs. Although the relative expression levels of two mRNAs and four lncRNAs were significant, their expression trends were consistent with those of RNA-seq, indicating the accuracy of the RNA-seq data.

Fig. 2.

Validation of the RNA-seq results was performed by RT-qPCR. DE mRNAs (A), lncRNA (B), lncRNA (C), and miRNAs (D) were confirmed by RT-qPCR. Expression levels of mRNAs (E), lncRNAs (F), circRNAs (G) and miRNAs (H) were displayed by RNA-seq. “*” means p $<$ 0.05; “**” means p $<$ 0.01; “***” means p $<$ 0.001.

3.4 Functional Enrichment Analysis of DE mRNAs

In this study, using high-throughput sequencing, 303 genes were differentially expressed between the two groups, of which 284 genes were annotated. In addition, the GO enrichment result of DE mRNAs showed that 303 DE mRNAs were enriched with 1019 GO terms, and significant enrichment was observed for 210 GO terms (p $<$ 0.05). Most of these terms were related to neurotransmitter and ovarian development, for example, neuropeptide signaling pathway, regulation of dopamine secretion, and calcium ion binding (Fig. 3A and Supplementary Table 6). Based on the findings of this study, some DE mRNAs were found to be enriched in GO terms related to reproduction, including those describing neurotransmitter secretion, the dopamine receptor signaling pathway, and the dopamine biosynthetic process for downregulated genes as well as the regulation of hormone activity and dopamine secretion for upregulated genes.

Fig. 3.

Functional analysis by DE mRNA and DE lncRNA. (A) DE mRNAs significantly enriched in the top 20 GO terms. (B) DE mRNAs significantly enriched in the top 20 KEGG pathways. (C) DE lncRNAs significantly enriched in the top 20 GO terms. (D) DE lncRNAs significantly enriched in the top 20 KEGG pathways.

KEGG pathway annotation results revealed that 73 of 303 DE mRNAs were enriched in 57 KEGG pathways, some of which were related to the formation of ovarian follicles, such as neuroactive ligand–receptor interactions, steroid hormone biosynthesis, calcium signaling, and extracellular matrix (ECM) receptor interaction (Fig. 3B and Supplementary Table 7).

3.5 Functional Enrichment Analysis of DE lncRNAs

In the present work, 328 potential target protein-coding genes of differential lncRNA transcripts were detected. Based on the results of enrichment analysis, 141 target genes were significantly enriched in 173 GO terms, including 20 terms for cellular component, 35 terms for molecular functions, and 118 terms for biological process (p $<$ 0.05; Fig. 3C and Supplementary Table 8). Some of these terms were related to binding, ion transport regulation, and signal transduction gene expression. In particular, 35 DE lncRNA target genes were annotated as components of 8 significant pathways (Fig. 3D and Supplementary Table 9), including cytokine–cytokine receptor interaction, steroid hormone biosynthesis, and neuroactive ligand–receptor interaction.

3.6 Functional Enrichment Analysis of DE circRNAs

The source genes of circRNA were used for GO enrichment analysis in this experiment. A total of 255 circRNA source genes were enriched to 43 significant GO terms (p $<$ 0.05; Supplementary Fig. 1A and Supplementary Table 10). Analysis indicated that the majority of circRNA source genes were involved in cellular fraction and molecular function, such as the cytoplasm, insulin receptor binding, and cyclin-dependent protein serine/threonine kinase activity. Likewise, a total of 74 source genes of DE circRNA were enriched in 19 KEGG pathways. In particular, the KEGG results indicated that most DE circRNAs source genes were enriched in ECM receptor interaction, focal adhesion, MAPK, and Wnt signaling pathways (Supplementary Fig. 1B and Supplementary Table 11).

3.7 Target Prediction of DE miRNAs and miRNA–mRNA Network Construction

A total of 16,541 target genes were predicted for 1028 miRNAs, 1259 of which were functionally annotated. Based on the abovementioned results, functional analysis showed that 576 target genes were significantly enriched in 360 GO terms (Supplementary Fig. 1C and Supplementary Table 12). In addition, 152 target genes were significantly enriched in 22 KEGG pathways (Supplementary Fig. 1D and Supplementary Table 13). Target genes were primarily enriched in the nucleus, cytoplasm, and membrane in GO terms. In addition, three known reproduction-related pathways were found to be abundant, such as the p53 signaling pathway, focal adhesion, and Wnt signaling pathway. The results showed that the Wnt signaling pathway and focal adhesion functioned in relation to reproduction in chicken.

Moreover, the miRNA–mRNA regulatory network (MMRN) was constructed to illustrate the interaction regulation between DE miRNAs and DE mRNAs, which plays a role in the physiological function of the ovary in chicken (Fig. 4A and Supplementary Table 14). Based on the MMRN, a single miRNA could negatively regulate multitarget genes. For example, upregulated gga-miR-6643-5p could downregulate more than 10 of its target genes (SYBU, CLEC2B, NECTIN4, STK32C, SOX10, EXOC3L1, EVA1C, MGLL, and RGS4).

Fig. 4.

Critical networks of DE coding and non-coding RNA interactions in ovarian tissue. (A) Key miRNA–mRNA regulatory diagram between comparisons. Negatively co-expressed mRNA-miRNA pairs with SCC $<$ –0.7 were selected. MiRNAs are represented by rectangular nodes and mRNAs by circular nodes. The illustration of circRNA–miRNA–mRNA (B) and lncRNA–miRNA–mRNA (C) network diagram. The lncRNAs/circRNAs, miRNAs, and mRNAs were represented by square, triangular and circular respectively. Transcripts in red show upregulation, and transcripts in green show downregulation.

3.8 lncRNA/circRNA–miRNA–mRNA Regulatory Network Construction

DE ncRNAs and mRNAs tracked from RNA-seq were used to construct the circRNA–miRNA–mRNA (Fig. 4B) and lncRNA–miRNA–mRNA (Fig. 4C) interaction networks (Supplementary Table 15). First, the circRNA–miRNA–mRNA network shows 76 nodes and 54 ceRNA pairs among 23 mRNAs, 17 miRNAs, and 36 circRNAs. Consequently, the lncRNA–miRNA–mRNA network contains 43 nodes and 18 ceRNA pairs, interacting with 17 mRNAs, 13 miRNAs, and 13 lncRNAs. Within both networks, a series of genes and ncRNAs (ENSGALG00000020615 [STK32C], ENSGALG00000050091 [CLEC2B], ENSGALG00000051188 [MR1], gga-miR-1736-3p, and gga-miR-6643-5p) may affect the laying production.

The gene functional enrichment result of 37 DE mRNAs from the ceRNA network further revealed that 21 DE mRNAs were present among 120 significantly enriched GO terms (p $<$ 0.05) with 67 biological process terms, 19 cellular component terms, and 34 molecular functional terms (Supplementary Fig. 2A and Supplementary Table 16). Only two significantly enriched pathways (steroid hormone biosynthesis and neuroactive ligand–receptor interaction) were annotated by CYP2D6, GRIA1, and TSPO2 genes (Supplementary Fig. 2B and Supplementary Table 17). Several genes, such as CYP2D6, FST and BF2, have also been primarily involved in steroid hormone biosynthesis, TGF- $\beta{}$ signaling pathway, and phagosome.

4. Discussion

Reproductive regulation in layers is a complicated biological process, and it is regulated by tight genes, including coding and non-coding RNAs [22], as well as various signaling pathways. In this study, various genes that were related to the regulation of ovarian development were identified. For example, VIP promotes the secretion of Prolactin (PRL) by the anterior pituitary gland through relative receptors on the membrane of the adenopituitary gland to participate in ovulation [23]. In poultry, PRL is a key hormone in the onset and maintenance of nesting behavior, and inoculation with PRL can inhibit the development of ovarian follicles and reduce egg production [24].

Furthermore, many DEmRNAs were enriched in neuroactive ligand–receptor interactions, steroid hormone biosynthesis, and calcium signaling pathways, all of which have been linked to egg production. These pathways play an important role not only in ovarian follicular development and normal ovarian function, but also in the regulation of egg yolk deposition, particularly the neuroactive ligand–receptor interaction pathway [25]. Ye et al. [26] found that the neuroactive ligand–receptor interaction was important for hormonal regulation in reproductive cycle transitions of duck. Another study has shown that this pathway affects egg production in geese through ovarian metabolic function. In the present study, AGTR1 was enriched in the neuroactive ligand–receptor interaction pathway. In addition, AGTR1 can constrict blood vessels in atretic follicles, thereby reducing vascular density and causing poor nutrient transport and follicular development [27]. Meanwhile, another DE mRNA, MC2R, has been proven to regulate the physiological response of animal ovaries by participating in the synthesis of steroid hormones [28].

LncRNAs have been implicated in numerous studies pertaining to reproduction, such as gonadal development and hormone regulation [29]. Notably, the target genes (HSD3B1, NR5A2, and INHBA) of three lncRNAs (ENSGALT00000105270, ENSGALT00000091851, and ENSGALT00000107167, respectively) were found from this study. A previous study has suggested that HSD3B1 is a key member of the steroid hormone family, which is associated with reproductive traits in livestock and poultry [30]. Moreover, NR5A2 has been proven to regulate not only the expression level of the steroid factor, but also the physiological process of steroid hormone secretion, follicle growth, and ovulation [31, 32]. The litter size of mice was significantly decreased with NR5A2 knockout [33, 34], and the variants of this gene were related to the lamb size in Songliao black pigs [35], indicating that the NR5A2 gene played a major role in regulating female reproductive performance. Finally, INHBA, as a member of the TGF- $\beta{}$ family, could affect follicle growth and development by regulating the synthesis and secretion of follicle-stimulating hormones [36]. A previous study has found that knocking out INHBA and inhibiting subunit $\beta{}$ B (INHBB) in granulosa cells affect activin and hinder production, resulting in sterile mice with a complex ovarian phenotype [37].

In this study, a large number of the host genes of DEcircRNAs are widely recognized as associated with animal reproductive regulation. In detail, MAPK3, which is the host gene of novel_circ_001023, was confirmed to be involved in primordial follicle activation [38]. MAPK3 is a major signal transduction molecule in the ERK/MAPK pathway, which can regulate a variety of nuclear transcription factors, cell proliferation, and apoptosis by phosphorylating cytoplasmic proteins [39]. Reports indicated that the expression or functional activity of MAPK3 is associated with the metastasis, occurrence, and drug resistance of various tumors [40]. The host gene EGFR of novel_circ_005081 was widely involved in mammalian reproduction, promoting cell differentiation and proliferation, and oocyte maturation [41, 42]. In addition, EGFR serves as a messenger molecule to activate the PKC signaling pathway that promotes oocyte division and growth in the preovulatory follicle, thereby affecting follicular development [43]. The overexpression of EGFR can hinder embryonic implantation and development by promoting progesterone and estrogen secretion [44].

A series of studies demonstrated that miRNAs were involved in many biological regulatory processes, including ovarian function, follicular development and atresia, cell proliferation and apoptosis, steroid biosynthesis, and ovarian cancer [45, 46, 47]. For example, miR-26b increased the number of DNA breaks and facilitated apoptosis of porcine granulosa cell by targeting the ATM gene in vitro [48]. MiR-31 and miR-143 act on synergistically steroid hormone synthesis to inhibit apoptosis in bovine granulosa cells via FSHR [49]. Further analysis of DE miRNA target genes screened in this assay revealed that miR-3525, which was significantly upregulated in CK groups, was negatively correlated with the expression of its predicted target gene, GRIA1. Previous studies have demonstrated that miR-3525 targets PDLIM3 to inhibit the proliferation and differentiation of skeletal muscle satellite cells through the p38/MAPK signaling pathway [50]. Interesting research has revealed that GRIA1 polymorphism was related to the number of antral follicle and fertility in cows [51].

In addition, miRNAs play a significant role in the regulation of ovarian function and create a sophisticated network of regulatory control systems in living things [52]. As an important factor of the P53 pathway, miR-34 family members can regulate various cellular processes, including apoptosis, aging, and differentiation. For example, miR-34b and miR-34c are targets of two well-known miRNAs that inhibit adherence-independent growth and proliferation in mouse ovarian surface epithelial cells [53]. Studies have shown that miR-34c-5p overexpression results in the increased level of apoptosis and decreased proliferation of granulosa cells, whereas adding a miR-34c-5p inhibitor to granulosa cells led to the opposite trend [54]. This finding indicated that gga-miR-34b and gga-miR-34c played important roles in the reproductive management of layers.

Rich studies confirmed that the ceRNA regulatory network participates in the maintenance of normal physiological state as well as the occurrence and development of diseases [55, 56]. For example, LncRNA-MALAT1 regulated apoptosis of granulosa cells and 17 $\beta{}$ -estradiol synthesis in mice by regulating the miR-205/CREB1 axis [57]. LncRNA-m433s1 regulates FSH secretion with a miRNA sponge in anterior pituitary cells of male rats [58]. Notably, novel_circ_000985-miR-449b-5p-LCK may be an important regulatory network in the circRNA–miRNA–mRNA network. LCK is a crucial protein-coding gene screened in the ceRNA network and an important reproduction-related gene, which is associated with follicular development and differentiation [59]. In pigs, LCK was considered as a candidate gene to fecundity [60]. MiR-449b regulates not only DNA synthesis [61], but also steroid synthesis in granulosa cells by targeting IGF2BP3 [62].

Moreover, novel_circ_002475-miR-460b-3p-FST may be another reproduction-related ceRNA network. In particular, novel_circ_002475 in this network was expressed only in CK. Previous studies have shown that the expression level of miR-460 was higher in broody groups than in egg-laying groups [63]. Meanwhile, Qiu et al. [64] confirmed that miR-460 is a key factor associated with egg production in ducks. In addition, Wasti et al. [65] found that FST was expressed in the lumen and glandular epithelium of the uterine tissue in layers, which affected the mineralization of eggshells. Furthermore, FST has been shown to suppress the proliferation of granulosa cells and inhibit the growth of primary cells [66]. The results indicate that these genes may influence egg production in layers through the ceRNA regulatory network.

5. Conclusions

This study described the expression patterns of various types of RNA in detail. Many identified lncRNAs, circRNAs, miRNAs, and mRNAs provide a reference for further research on the regulatory mechanism of ceRNAs. These findings also provide new insights into the molecular mechanism of ovarian egg production, although further investigations are necessary to validate the ceRNA mechanism of ncRNAs and mRNAs.

Availability of Data and Materials

All data generated or analyzed during this study are included in this published article.

Author Contributions

These should be presented as follows: ML and GE designed the research study. ZW performed the research. CheL and ChaL helped analyze and validate the data. DS provided funding. ZW analyzed the data. ZW wrote the manuscript. All authors contributed to editorial changes in the manuscript. All authors read and approved the final manuscript.

Ethics Approval and Consent to Participate

The animal study protocol was approved by the Committee on the Ethics of Animal Experiments of Southwest University (IACUC-20211020-08) and the Animal Protection Law of China.

Acknowledgment

The authors would like to thank all staff at the chicken experimental for animal care and sample collection.

Funding

This research was funded by Chongqing Technology Innovation and Application Development Special Key Projects, grant number cstc2019jscx-gksbX0107.

Conflict of Interest

The authors declare no conflict of interest.

Supplementary Material

Supplementary material.rar

References

[1]

Du Y, Liu L, He Y, Dou T, Jia J, Ge CJBPS. Endocrine and genetic factors affecting egg laying performance in chickens: a review. British Poultry Science. 2020; 61: 538–549.