Xiaolong Du,Yufang Liu,Xiaoyun He,Lin Tao,Meiying Fang,Mingxing Chu
Uterus proliferative period ceRNA network of Yunshang black goat reveals candidate genes on different kidding number trait
goat,uterus,proliferative period,ceRNA network,cell adhesion
Pregnancy loss that occurs in the uterus is an important and widespread problem in humans and farm animals and is also a key factor affecting the fecundity of livestock. Understanding the differences in the fecundity of goats may be helpful in guiding the breeding of goats with high fecundity. In this study, we performed RNA sequencing and bioinformatics analysis to study the uterus of Yunshang black goats with high and low fecundity in the proliferative period. We identified mRNAs, long non-coding RNAs (lncRNAs), and microRNAs (miRNAs) by analyzing the uterine transcriptomes. The target genes of the identified miRNAs and lncRNAs were predicted, and miRNA–mRNA interaction and competitive endogenous RNA (ceRNA) networks were constructed. By comparisons between low- and high-fecundity groups, we identified 1,674 differentially expressed mRNAs (914 were upregulated, and 760 were downregulated), 288 differentially expressed lncRNAs (149 were upregulated, and 139 were downregulated), and 17 differentially expressed miRNAs (4 were upregulated, and 13 were downregulated). In addition, 49 miRNA–mRNA pairs and 45 miRNA–lncRNA pairs were predicted in the interaction networks. We successfully constructed a ceRNA interaction network with 108 edges that contained 19 miRNAs, 11 mRNAs, and 73 lncRNAs. Five candidate genes (PLEKHA7, FAT2, FN1, SYK, and ITPR2) that were annotated as cell adhesion or calcium membrane channel protein were identified. Our results provide the overall expression profiles of mRNAs, lncRNAs, and miRNAs in the goat uterus during the proliferative period and are a valuable reference for studies into the mechanisms associated with the high fecundity, which may be helpful to guide goat to reduce pregnancy loss.
Uterus proliferative period ceRNA network of Yunshang black goat reveals candidate genes on different kidding number trait

Pregnancy loss that occurs in the uterus is an important and widespread problem in humans and farm animals and is also a key factor affecting the fecundity of livestock. Understanding the differences in the fecundity of goats may be helpful in guiding the breeding of goats with high fecundity. In this study, we performed RNA sequencing and bioinformatics analysis to study the uterus of Yunshang black goats with high and low fecundity in the proliferative period. We identified mRNAs, long non-coding RNAs (lncRNAs), and microRNAs (miRNAs) by analyzing the uterine transcriptomes. The target genes of the identified miRNAs and lncRNAs were predicted, and miRNA–mRNA interaction and competitive endogenous RNA (ceRNA) networks were constructed. By comparisons between low- and high-fecundity groups, we identified 1,674 differentially expressed mRNAs (914 were upregulated, and 760 were downregulated), 288 differentially expressed lncRNAs (149 were upregulated, and 139 were downregulated), and 17 differentially expressed miRNAs (4 were upregulated, and 13 were downregulated). In addition, 49 miRNA–mRNA pairs and 45 miRNA–lncRNA pairs were predicted in the interaction networks. We successfully constructed a ceRNA interaction network with 108 edges that contained 19 miRNAs, 11 mRNAs, and 73 lncRNAs. Five candidate genes (PLEKHA7, FAT2, FN1, SYK, and ITPR2) that were annotated as cell adhesion or calcium membrane channel protein were identified. Our results provide the overall expression profiles of mRNAs, lncRNAs, and miRNAs in the goat uterus during the proliferative period and are a valuable reference for studies into the mechanisms associated with the high fecundity, which may be helpful to guide goat to reduce pregnancy loss.

In female agricultural animals, fecundity traits are the main indicator of their economic value, especially the number of fetuses (1). Pregnancy loss that occurs in the uterus is an important and widespread problem in humans and farm animals, especially in low-ovulatory animals (2). During natural conception, approximately 75% of the final pregnancy failures are due to the inappropriate environment in the uterus (3, 4). Differences between uteruses are the main consideration when exploring pregnancy loss. The Yunshang black goat is a specialist mutton goat that was bred in China. These goats have good economic benefits, but intravarietal kidding numbers vary significantly (5). Hence, it is important to study the effect of uterine variability on kidding traits. The uterus comprises the serous membrane, myometrium, and endometrium. The endometrium consists of luminal and glandular epithelial cells surrounded by supporting stromal cells (6). The main function of the uterus is to maintain pregnancy, and the endometrium promotes implantation and the decidua and supports the growth and development of the embryo (7). Embryo implantation in humans and rodents mainly involves embryonic trophectoderm cells penetrating the endometrial epithelial layer and invading the stroma to form the placenta and establish pregnancy (8). However, the pregnancy process in ruminants is entirely different from that of humans and rodents. In ruminants, trophoblast cell fusion involves the development of finger-like villi or papillae between the caruncles of the uterus. These structures penetrate the mouths of superficial ducts of the uterine glands at days 15–18 after mating. Although similar features were described for the cow conceptus from day 15 of pregnancy, goats did not exhibit the same phenomenon (9, 10). The main difference is that, in goats, superficial implantation occurs; that is, the embryonic trophectoderm cells adhere to the endometrial caruncle and endometrial epithelial cells fuse to form multinucleated syncytial plaques, then connect with each other to form a cotyledon type placenta (11). Complex mucins bind the embryo to the uterus. A study using the uterine gland knockout (UGKO) ewe model found that the reduction or absence of endometrial-derived cohesin was the main cause of recurrent miscarriage in UGKO ewes (12). The endometrium proliferates and secretes with the cyclical changes in the ovary, and whether embryo implantation is successful is directly related to regression of the corpus luteum and recovery of the endometrium (13). However, some scientists consider the uterus to be the cause of the cyclical changes in the ovary because the uterus can secrete prostaglandins to promote luteolysis. Indeed, when the uterus was removed during the luteal phase, the estrus period was delayed by 5 months, which is the same as the pregnancy time of goats (10). MicroRNAs (miRNAs) that were differentially expressed in the proliferative and secretory phases of human endometrial were reported in 2008 (14). Subsequent research found that changes in miRNA expression levels in endometrial tissue may be influenced by estradiol and progesterone, indicating that miRNAs may be involved in the reproductive physiological processes of the endometrial tissue (15). The involvement of miR-182 in the establishment of endometrial receptivity by regulating the pleiotrophin (PTN) and homeobox A10 (HOXA10) genes was found in dairy goat endometrial epithelial cells (16, 17). More recently, the roles of long non-coding RNAs (lncRNA) in diverse biological activities, such as embryogenesis and transcriptional regulation, have been recognized, and many lncRNAs with functions in goat endometrium during pregnancy recognition have been discovered (18, 19). However, that lncRNAs and miRNAs jointly affect mRNA expression in goat endometrium has rarely been reported. In this study, we performed miRNA and lncRNA sequencing and constructed an miRNA–lncRNA–mRNA interaction network and competitive endogenous RNA (ceRNA) network to explore the factors that affect goat kidding numbers from a new perspective.

Ethics approval (no. IAS2019-63) was granted by the Animal Ethics Committee of the Institute of Animal Sciences of Chinese Academy of Agricultural Sciences (IAS-CAAS) (Beijing, China).

A total of 10 non-pregnant adult female Yunshang black goats (aged 3–5 years old; weight, 41–63 kg) were selected randomly from the Yixingheng Animal Husbandry Technology Co., Ltd, Tuanjie Township Base in Kunming City (24° 23′ north latitude), Yunnan Province, China, and were used in this study. The goats were separated into a high-fecundity group (HF group; average kidding number, 3.00 ± 0.38, n=5) and a low-fecundity group (LF group; kidding number, 1.32 ± 0.19, n=5) (p <0.01, t-test) according to their kidding number records. All uterine tissue layers from similar segments with distance from the cervix were collected in the proliferative period after synchronous estrus processing and stored at −80°C until used for RNA extraction. All pretreatments were based on those of a previous study (20). Blood samples were collected for hormone determination, which was measured by competitive radioimmunoassay. Total RNA was isolated from 10 uterine tissue samples (containing endometrium and stroma) with TRIzol (Invitrogen, Carlsbad, CA, USA) according to the manufacturer’s protocol. The quality of the RNA samples was determined using a NanoDrop 2000 spectrophotometer (Thermo Scientific, Wilmington, DE, USA), Agilent Bioanalyzer 2100 system (Agilent Technologies, Palo Alto, CA, USA), and RNA Nano 6000 Assay (Agilent Technologies, CA, USA). The OD 260/280 was 1.8–2.2, and the RNA integrity number (RIN) of all samples was >7.

An Epicenter Ribo-Zero™ Removal Kit (Epicenter, Madison, WI, USA) was used to remove rRNA, and rRNA-free residues were cleaned by ethanol precipitation. Then, sequences that passed the quality inspection were used for library construction and sequencing. The lncRNA and mRNA library was constructed using 3 μg of the total RNA and a NEBNext® UltraTM Directional RNA Library Prep Kit (NEB, Ipswich, MA, USA) according to the manufacturer’s instructions. The Illumina NovaSeq 6000 platform (Illumina, San Diego, CA, USA) was used for PE150 (paired-end 150 bp) sequencing. The miRNA library was built using a Small RNA Sample Prep Kit (Tiangen, Beijing, China). The total RNA was used as a template, and cDNA was synthesized with a special adapter. After PCR amplification, the products were purified by polyacrylamide gel electrophoresis (PAGE), and 140–160-bp sequences were recovered to construct an miRNA sequencing library. The above platforms were used for SE50 (single-end 50 bp, SE50) sequencing. All the sequencing were outsourced to Wuhan Frasergen Gene Information Co., Ltd. (Wuhan, China).

The original image file from the sequencer was converted into a sequence file that contained mainly sequence information and sequencing quality information using CASAVA software. The raw data (FASTQ) were filtered for use in subsequent analysis as follows: 1) reads that contained the adapters were removed; 2) reads with N content >1% of the length of the read itself and reads that constituted the paired-end pair were filtered out; and 3) reads with low-quality bases (≤20) that exceeded 50% of the length of read itself and reads that constituted the paired-end pair were filtered out. Quality evaluation of the clean reads included calculation of sequencing quality (Q20 and Q30) base distribution (GC content). Clean reads that passed the quality evaluation were mapped to the reference genome (Capra hircus, ARS1) using HISAT2 (v2.1.0). StringTie (v1.3.5) was used to assemble the transcripts (21, 22). HISAT2 was run with parameters “-x -no-unal -un-conc”, “-rna - strandness RF” and “-dta -t -p 4,” String Tie with “-e -B”, “-G ref.gtf -rf -1;” the other parameters were set as the default. The fragments per kilobase per million reads (FPKMs) for each gene were calculated based on the length of the gene and read counts mapped to the gene (23).

We screened the assembled transcripts using the following criteria to identify lncRNAs and mRNAs: 1) transcript lengths ≥200 bp and number of exons >1; 2) transcripts were compared with mRNA and other non-coding RNA sequences (rRNA, tRNA, snRNA, and snoRNA) in known databases using Cuffcompare, and rRNA and tRNA sequences were removed (23). When an annotated lncRNA overlapped the exon region of a sequence in the database, it was included as a known lncRNA in the subsequent analysis (3). For novel lncRNAs and mRNAs, three coding potential prediction software [CNCI (24), CPC2 (25), and PLEK (26)] were used to predict the coding potential of the transcripts, and the intersection of the three prediction results was used. To identify miRNAs, clean high-quality reads of 18–35 nt were compared against miRNA sequences in the miRBase (v22) database (27) using Bowtie2. Clean high-quality reads not identified as known miRNA sequences were compared with other non-coding RNA sequences (rRNA, tRNA, snRNA, and snoRNA) in Rfam (14.2) (28) and classified and annotated. RepeatMasker was used to identify repeat sequences in the clean reads that were not annotated as non-coding RNAs or known miRNAs (29). After evaluating these sRNAs and using the location information of exons and introns of the gene from the matching (100% positional overlap) results, the sRNAs were annotated (30). For the sRNA reads that did not match the above-known annotation type, miREvo and miRDeep2 (31) were used to predict novel miRNAs, and the sequence, expression, and structure information of these miRNAs were analyzed. Sequences with different core sequences, precursor sequences, and positions were considered to be new miRNAs.

The count of reads mapped to each transcript was obtained using StringTie and normalized by computing TPM and FPKM values to indicate miRNA and mRNA/lncRNA expression levels. For the paired-ended sequencing, paired-end reads from the same fragment were counted as one fragment. The DESeq2 R package was used for differential expression analysis, and adjusted p <0.05 and |log2FC| >1 were used as the threshold (32).

The Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) functional analyses of the target genes of the differentially expressed lncRNAs (DELs) and differentially expressed miRNAs (DEMs), and the differentially expressed mRNAs (DEGs) were performed using GOSeq (Release 2.12) and KOBAS (v2.0) (33), respectively. The hypergeometric test method was implemented to assess significantly enriched GO terms and KEGG pathways, and p-values (t-test) <0.05 were considered significant.

To analyze the interaction between miRNA and mRNA, we used MiRanda and qTar to predict the target genes of the miRNAs and used the intersection of the two software as the result (34). Correlated miRNA–mRNA pairs were identified by Pearson correlation analysis. Because lncRNAs can act as miRNA sponges, they can inhibit the regulatory effect of miRNAs on their target mRNAs; therefore, we analyzed the relationship between miRNAs and lncRNAs in a similar way.

For the RT-qPCR analysis, we used a PrimeScript™ RT Reagent Kit (TaKaRa, Dalian, China) and miRcute Plus miRNA First-Strand cDNA Kit (TIANGEN) synthetic-detection templates. The expression of the selected mRNAs and miRNAs was quantitatively analyzed using a SYBR Green qPCR Mix Kit (TaKaRa) and miRcute Plus miRNA qPCR Kit (TIANGEN) according to the manufacturer’s instructions. The RT-qPCRs were performed on a RocheLight Cycler®480 II system (Roche Applied Science, Mannheim, Germany). The PCR primers were synthesized by Shanghai Sangon Biotech. The relative expression was calculated by the 2−ΔΔCP method. The goat PRL19 gene and U6 were used as reference genes. The primer sequences are listed in Supplementary Material .

The SPSS 20.0 software was used for statistical analysis. The statistical significance of the data was tested by performing paired t-tests. The results are presented as the means ± SEMs of three replicates, and statistical significance was represented by p-values (*p < 0.05; **p<0.01).

Test samples were selected based on an evaluation of kidding records over two birthings. At the time of sample collection, we also recorded body weight data, hormone levels, and total number of follicles in the left and right ovaries after dissection. The phenotypic data for the five high-fecundity and five low-fecundity Yunshang black female goats in this study are shown in Table 1 . Except for significant phenotypic date, kidding number, and average kidding number (p <0.01, t-test), some not significant phenotypic differences are listed in Table 1 between the high- and low-fecundity groups (p>0.05, t-test).

Total RNA was extracted from 10 uterine tissues and used to construct two sets of libraries: one to identify lncRNAs and mRNAs, and the other to identify small RNAs (sRNAs). After processing and quality control of the raw data, an average of 118.3 million and 24.7 million clean reads were obtained from the lncRNA/mRNA and sRNA libraries, respectively ( Table 2 ). The mapping rates of the lncRNA/mRNA and sRNA reads against the goat genome (Capra hircus, ARS1) were an average of 96.97% and 97.72%, respectively. The quality of the lncRNA and mRNA sequencing results was assessed by calculating the Q30 and GC content of the two reads that constituted the paired-ends. The results showed that the quality of the clean reads met the requirements for subsequent analysis. The other filtering criteria used to identify candidate lncRNAs, mRNAs, and miRNAs are provided in Supplementary Material S1 .

A total of 67,424 mRNA were identified and successfully mapped to the goat genome. StringTie was used to compare these mRNAs to the reference annotations in the goat genome. The results showed that 65% of the mRNAs were annotated in the reference database; therefore, we defined them as known RNAs. The remaining mRNAs (35%) for which no annotations were found were defined as novel mRNAs. We used the same method to identify sequences longer than 200 bp with more than two exons as lncRNAs. We found that only 7% of the identified lncRNAs were known lncRNAs; the remaining 93% were defined as novel lncRNAs ( Figure 1A ). The predicted genomic position of the novel lncRNAs in the goat genome showed that 8,648 were intergenic lncRNAs, 12,389 were intronic lncRNA, and 2,928 were antisense lncRNA ( Supplementary Material S2 ). We used StringTie to assemble the annotated reads into transcripts so that their detailed chromosome locations could be obtained. We found that the proportions of mRNAs and lncRNAs were similar on most of the chromosomes. Among the autosomal chromosomes, chromosomes 3, 5, 7, 18, and 19 had the highest proportions of mRNAs and lncRNAs, followed by chromosomes 1, 2, 10, and 11. Chromosomes 12, 20, 24, 27, and 28 had the lowest proportions of mRNAs and lncRNAs. We also found that 0.02% of the transcripts were mitochondrial mRNAs (n=13); no mitochondrial lncRNAs were identified ( Figure 1C ). To identify miRNAs, we compared the sequenced reads against a number of sRNA databases that had different annotation information and defined the sRNAs in the following order: known miRNAs > rRNAs > tRNAs > snRNAs > snoRNAs > repeats > novel miRNAs > exons > introns. We found that known miRNA accounted for 74.36% and 77.07% of the average total sRNAs in the high- and low-fecundity groups, respectively. Other sRNAs (tRNA, rRNA, snRNA, and snoRNA), which were identified by comparisons against known databases such as Rfam, accounted for a very small proportion of the average total sRNAs in the two groups. Most of the remaining sRNAs accounted for 23.26% and 19.85% of the average total sRNAs in the high- and low-fecundity groups, not including repeats or degraded fragments of genomic exons and introns, which accounted for very small proportions of the average total sRNAs in the two groups ( Figure 1B ; Supplementary Material S2 ).

To normalize the counts for the clean paired-end reads, we calculated FPKM and transcripts per kilobase million (TPM) to obtain the expression levels of the transcripts for subsequent quantitative analysis. We found that the overall distribution of transcripts in the high- and low-fecundity groups was almost the same, whether by maximum, minimum, or median analysis ( Figures 2A, B ). The consistency of the expression data also confirmed the accuracy of our sequencing data. The FPKM distribution of the lncRNAs and mRNAs was significantly different but in line with the expected result ( Figure 2C ). The number of exons in the lncRNA and mRNA transcripts was very different; most of the lncRNAs had two exons, whereas most of the mRNAs had from two to seven exons ( Figure 2D ; Supplementary Material S3 ).

We identified differentially expressed lncRNAs (DELs) and mRNAs (DEGs) by comparing the average expression levels between the low- and high-fecundity groups with a cutoff threshold of |log2fold change (FC)| ≥1 and q <0.05. A total of 1,674 DEGs and 288 DELs were identified; 914 DEGs were upregulated and 760 DEGs downregulated, and 149 DELs were upregulated and 139 DELs downregulated. We also identified 17 DEMs (four were upregulated and three were downregulated) in the low-fecundity vs. high-fecundity comparison with a cutoff threshold of |log2FC| ≥1 and p <0.05 ( Figure 3D ). The Euclidean distance was calculated based on log base 2 of the DEG, DEL, and DEM expression levels in each sample, and then, the overall clustering results of the samples were obtained using a systematic hierarchical clustering method ( Figures 3A–C ; Supplementary Material S4 ).

GO and KEGG enrichment analysis were performed using the DEGs and the target genes of the DELs and DEMs. The enrichment results were screened using different thresholds (GO, q ≤0.01; KEGG, p ≤0.05) because the numbers of enriched GO and KEGG terms were very different (DEGs GO, n=124; KEGG, n=17; target genes of DELs GO, n=1,135; KEGG, n=100; target genes of DEMs GO, n=304; KEGG, n=41.). By overlapping the enriched GO terms for the DEG and DEL target genes, we identified 54 terms that were common to both, whereas only four GO terms for the DEG and DEM targets overlapped. The top 10 overlapping terms under the three main GO categories, namely, biological process, cellular component, and molecular function, were analyzed, and the results are shown in the column chart in Figure 4A . They include two biological process terms (GO: cellular macromolecule metabolic process, GO: cellular response to epidermal growth factor stimulus), eight cellular component terms (GO: response to epidermal growth factor, GO: cytoplasmic side of lysosomal membrane, GO: extracellular region part, GO: myelin sheath, GO: extracellular space, GO: cell projection membrane, GO: extracellular matrix, GO: macromolecular complex), and one molecular function MF term (GO: translation factor activity, RNA binding). It is worth noting that the cellular component term extracellular matrix and the biological process term cellular response to endogenous stimulus, which are not listed in Figure 4A , overlapped with terms in the top 10 overlapping terms under the three main GO categories. We constructed a bubble chart based on all the enriched KEGG pathways and found that only 17 pathways were significantly enriched (p ≤0.05) ( Figure 4B ); nine and eight of these pathways were for the DEL and DEM target genes, respectively. Notably, four of these pathways were significantly enriched for the DEGs and DEL and DEM target genes, namely, adherens junction, focal adhesion, platelet activation, and PI3K-Akt signaling pathway ( Supplementary Materials S5 , S6 ).

To predict the targeting relationship of the known and novel miRNAs with the mRNAs and lncRNAs, we used MiRanda and qTar software. We calculated Pearson correlation coefficients and used negative correlations to construct miRNA–lncRNA and miRNA–mRNA regulatory networks. LncRNAs can bind to miRNAs and can act as sponges to regulate the effect of miRNAs on their target miRNA. Therefore, we also constructed a ternary ceRNA interaction network of lncRNA–miRNA–mRNA ( Figure 5 ; Supplementary Material S7 ). After screening, 49 pairs miRNA–mRNA pairs were predicted and the participating miRNAs were significantly differentially expressed between the low- and high-fecundity groups ( Figure 5A ). The novel_179–ITPR2, novel_179–PER3, novel_179–Novel_011082, and novel_179–RAVER2 pairs were of particular interest because the target mRNAs were also significantly differentially expressed between the low- and high-fecundity groups (ITPR2 q ≤0.01, PER3 q ≤0.05, Novel_011082 q ≤0.05, and RAVER2 q ≤0.05). Similarly, 45 miRNA–lncRNA were predicted and three of the participating lncRNAs were significantly differentially expressed between the low- and high-fecundity groups, namely, LNC_019255 (q ≤0 01), LNC_006767 (q ≤0 01), and LNC_019300 (q ≤0 01) ( Figure 5B ). On the basis of the relationships between the miRNAs, lncRNAs, and mRNAs, we successfully constructed a ceRNA interaction network with 108 edges that contained 19 miRNAs, 11 mRNAs, and 73 lncRNAs ( Figure 5C ). The transcripts of 10 genes in the network were significantly differentially expressed between the high- and low-fecundity groups, namely, PLEKHA7-XM_018059600.1 (log2FC=−9.332, q ≤0.01), FAT2-XM_018050131.1 (log2FC=−8.679, q ≤0.01), EDEM3-XM_018060712.1 (log2FC=−2.254, q ≤0.01), KANK2-XM_018051102.1 (log2FC=−1.727, q ≤0.05), MGAT5B-XM_018063721.1 (log2FC=1.543, q ≤0.01), FN1-XM_018059847.1 (log2FC=3.909, q ≤0.01), DIXDC1-XM_018059788.1 (log2FC=−3.580, q ≤0.01), SYK-XM_005684185.3 (log2FC=4.380, q ≤0.05), NRCAM-XM_018047218.1 (log2FC=8.804, q ≤0.01), and CGNL1-XM_018054269.1 (log2FC=8.938, q ≤0.05). The Pearson correlation coefficient results indicated that it was highly possible that the binding of novel_281 and LNC_007083 affected the expression of INVS mRNA and that the binding of chi-miR-326-3p and lncRNA affected the expression of SYK and NRCAM mRNAs. However, the differential expression analysis of the transcripts of INVS showed that only one of the INVS variants (GenBank: XM_013966062.2) was significantly differentially expressed; the transcript of the other INVS variant (GenBank: XM_013966063.2) was not significantly differentially expressed ( Supplementary Material S7 ).

To verify the accuracy and stability of sequencing results, RT-qPCR was used to compare the RNA-seq data. From the results, the quantitative results of mRNA and miRNA are consistent with the sequencing results ( Figure 6 ), which certified the reliability of the data produced by sequencing. Relative gene expression was calculated using the 2−△△CP method.

The Yunshang black goat is a new breed of mutton sheep that was bred in China using the black Nubi goat as the male parent and the Yunling black goat as the female parent after decades of crossbreeding (5, 20). Yunshang black goats have good meat production performance and delicious meat and are very popular in China, but high kidding number is the most wanted reproductive trait. Although kidding traits of Yunshang black goat have been subjected to strong positive selection during decades of cultivation (5), intravarietal significant differences in kidding number are still found. A statistical analysis of kidding numbers showed that the high-fecundity group had an average kidding number of 3.00 ± 0.38 and the low-fecundity group had an average kidding number of 1.32 ± 0.19 (p <0.05) (35). There was also a significant difference in the kidding number in this study but no significant difference in the number of follicles. Thus, pregnancy loss may be the key to affect fecundity. Many candidate genes have been shown to affect the litter size of sheep, including the main fecundity gene FecB, but when CRISPR-Cas technology was used to knock out one of the candidate genes BMPR-IB, it had no significant effect on the kidding number of goats (36). Until now, the major gene that affects the kidding number of goat has not been found. Many candidate genes related to kidding number have been discovered in goat, but the related molecular mechanism is not very clear; for example, some targeted genes (e.g., AMHR2, FGFR1, and SMAD2) under selection were reported in a previous study (37). Genomic regions linked to fecundity were indirectly identified by comparing two typical breeds and combining F-statistics (Fst) and heterozygosity (Hp) (38). Genomic selection sweep analysis and selection signatures analysis have been performed in single goat breeds with differing kidding numbers (39–42). In this study, we used Yunshang black goats with low and high kidding numbers to explore these differences and find candidate genes. Successful pregnancy depends on a well-developed embryo and uterus with the right environment (43). Before the embryo and uterus can establish a connection, they undergo a preparation process. For example, the cytokine IFNT (IFNT a type I IFN), which is secreted by the trophoblast cells of the conceptus and has strong anti-viral, anti-proliferation, and immune-regulating effects, is required by the ruminant mother to recognize the pregnancy signal (10, 44). The mother also needs to secrete various hormones to regulate the uterine environment to cope with implantation of the conceptus; estrogen and progesterone are the most critical of these hormones (45). Successful attachment of the embryo to the uterine lining is needed before subsequent development can occur. Different from other ruminants, in goats, the attachment method is superficial implantation, and therefore, the embryo does not invade the endometrium (11). Opposition of the conceptus involves the trophectoderm becoming closely associated with the endometrial luminal epithelium followed by unstable adhesion. After day 14 of mating, the filamentous conceptus appears to be immobilized in the uterine lumen, and the trophectoderm maintains close contact with the endometrial luminal epithelium (46). Close association of the apical membranes of both these cell types takes place; however, the conceptus can still be recovered intact from the uterus by lavage. Opposition of the blastocyst is ensured by the interdigitation of cytoplasmic projections of the trophectoderm cells and uterine epithelial microvilli (9). In ruminants, the openings of uterine glands are also sites of apposition (47). The trophoblast develops fingerlike villi or papillae between the caruncles, which penetrate the mouths of the superficial ducts of the uterine glands at days 15–18 after mating (9, 48). During their short life (the villi/papillae disappear by day 20), these trophoblastic differentiations were hypothesized to anchor the per-attachment conceptus and absorb phototrophic secretions of the glands (9). Similar features have been described for the cow conceptus from day 15 of pregnancy, but, curiously, the goat conceptus lacks trophoblast papillae. In this study, we performed enrichment analysis on DEGs and differentially expressed target genes and found that adhesion junction, focal adhesion, and related pathways such as epidermal growth factor were significantly enriched. A recent study that used cervical mucus of different breeds of ewes with high and low fecundity for physicochemical analysis and RNA sequencing found significant differences in the terminal modification of mucin between the different fecundities and that the differentially expressed genes were related to adhesion (49). Our sequencing results show that there were significant differences in the expression of hub mRNA between the miRNA–mRNA interaction network and the miRNA–lncRNA–mRNA ceRNA network. Among them, PLEKHA7, which encodes PLEKHA7 (pleckstrin homology domain-containing A7), is a zonula adhesion component closely connected with Nezha, which anchors the microtubules of epithelial cells to the zonula adhesion junction (50). Pulimeno et al. (51) showed that PLEKHA7 was located only in the apical part of epithelial cells in several tissues and that PLEKHA7 was involved in the recruitment of E-cadherin to enhance the potential of adhesion (52). The endometrium is made up of epithelial cells and, in goat, which mainly show surface attachment during embryo implantation, a cotyledon placenta is finally formed by adhesion to the endometrium. Therefore, adhesion ability has a key role for successful pregnancy. Whether our results can confirm that the difference between the high- and low-fecundity goats in our study is due to a difference in endometrial adhesion ability needs further analysis. Besides PLEKHA7, which is involved in adhesion, it has been reported that CircPlekha7, which is derived from PLEKHA7, also plays an important role in resisting pathological intrauterine adhesion (53). Fibronectin 1 (FN1), a large extracellular matrix adhesion glycoprotein composed of two similar subunits, has a variety of cell surface and extracellular ligand-binding domains and is involved in a many biological processes, including intrauterine adhesion (54). For example, FN1 and its receptors have different functions in mammalian reproduction and embryogenesis (sperm–oocyte interaction, implantation, and placenta formation), and most FN1 receptors are integrins that transmit cellular signals in both directions (55–57). The mRNA that encodes the spleen tyrosine kinase (SYK) protein, which is part of the integrin cluster (58), was differentially expressed in our study. One of its functions is to stimulate the phosphorylation of Toll-like receptor 4 (TLR4) (59), which recognizes bacterial lipopolysaccharide and induces inflammation and immune response. SYK also induces pathogen recognition, cell adhesion, and platelet activation (60). Embryo implantation is known to be an immune response to the mother. The surface attachment mode of goat embryos is also the attachment mode of other ruminant embryos and is very different from that of other mammals. The main difference lies in the fetal–maternal cell fusion process. There are three types of trophoblasts in the bovine placenta: mononucleate trophoblast cells (MTCs), binucleate cells (BNCs), and trinucleate cells (TNCs) (61). MTCs and BNCs have consistently been identified in the placental trophectoderm of many ruminants, whereas TNCs have been found only in bovine species, including cattle (62). Instead of TNCs, sheep and goats, which belonging to subfamily Caprinae (family Bovidae), develop multinucleated syncytial plaques (SyPs) (63). BNCs are formed by differentiation of MTCs by endoreduplication, whereas TNCs and multinucleated syncytial plaques are thought to be generated by cell-to-cell fusion between BNCs and maternal endometrial cells (64). Nakaya et al. (62) found that a retrovirus located in intron 18 of the bovine FAT2 gene played a significant role in bovine TNC formation. Coincidentally, we detected a difference in the expression of FAT2 between the high- and low-fecundity groups. In addition, the differentially expressed ITPR2 gene in our study encodes inositol 1,4,5-trisphosphate receptor type 2 (IP3R2), which is a member of the family of intracellular Ca2+ release channels and is located on the endoplasmic reticulum membrane (65). Studies using a gene knockout mouse model have shown that IP3R plays important roles in regulating many physiological processes, including embryonic survival (66), extraembryonic vascular development (67), T-cell development (68), and B-cell function (69). Allantoic-placental defects were detected in all IP3R2 and IP3R3 knockout embryos, and the fetal cardiovascular development was abnormal under the influence of placental defects (65). It has also been suggested that the Ca2+ signaling pathway may be involved in regulating the establishment of fetal–maternal connections (70).

We identified 1,674 differentially expressed mRNAs (914 upregulated and 760 downregulated), 288 differentially expressed lncRNAs (149 upregulated and 139 downregulated), and 17 differentially expressed miRNAs (4 upregulated and 13 downregulated) in low-fecundity vs. high-fecundity comparisons. By constructing interaction networks, we discovered 49 pairs in the miRNA–mRNA network and 45 pairs in the miRNA–lncRNA network. We successfully constructed a ceRNA interaction network with 108 edges that was composed of 19 miRNAs, 11 mRNAs, and 73 lncRNAs. Five candidate genes (PLEKHA7, FAT2, FN1, SYK, and ITPR2) related to cell adhesion or calcium membrane channel protein were identified. Our study provides the overall expression profiles of mRNAs, lncRNA, and miRNA in goat uterus during the proliferative period and offers new insights into the mechanisms associated with the high fecundity of goats, which may be helpful to guide goat to reduce pregnancy loss.

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/ Supplementary Material .

This manuscript follows the ARRIVE animal experiment report guidelines. The animal study was reviewed and approved by the Animal Ethics Committee of the Institute of Animal Sciences of the Chinese Academy of Agricultural Sciences (No. IAS2019-63). Written informed consent was obtained from the owners for the participation of their animals in this study.

Conceptualization: MC and MF. Methodology: XD, YL, MC, and MF. Validation: XD, YL, and XH. Formal analysis: XD and MC. Resources: XD, YL, LT, and MC. Investigation: XD, YL, LT, XH, MF, and MC. Writing—original draft preparation: XD and MC. Data curation: XD, YL, LT, MF, and MC. Supervision: MC. Project administration, MC. Funding acquisition: MC. All authors have read and agreed to the published version of the manuscript.