In China, PRRS and ring disease are one of the most important factors affecting litter size. Our experiment is to guide the pig breeding industry, so we chose sows without these two diseases. In addition, the sows we selected are of the same age, body condition, and feed with the same nutritional level. Most importantly, they have no history of reproductive diseases. We collected the litter size records of the Dabai Landrace binary cross sows from the Dahua breeding farm (Guangxi Hanshiwei Food Co., Ltd., from 2016 to 2018, a total of 8,557 litters), and used SPSS25.0 to perform a significance test on these data. After normal distribution conversion and inspection, it is found that the total litter size (12.9 ± 2.17) approximately obeys the normal distribution. The critical value of 15% right tail probability is 14.7 heads per litter, and the critical value of 15% left tail probability is 9.3 heads per litter. Ovarian tissues of pig under follicles of the same diameter were taken when static reflex occurred in pigs of the same parity.

Animal welfare and research plan were approved by Animal Care Committee of the Guangxi University. The ovaries of sows (Large White × Danish Landrace binary hybrid pig) were collected from a commercial pig farm (Nanning, Guangxi, China) 4 days after the fourth delivery. Six ovary tissues were chosen to perform RNA-Seq from pig with small (SLS, 8.48 ± 0.53/litter) and large (LLS, 16.19 ± 0.43/litter) litter size after the fourth delivery, respectively. The sample was were obtained and snap-frozen in liquid nitrogen.

We subjected the extracted RNA from tissue samples of six pig ovaries (three pig ovaries from pigs with small litter size; three pig ovaries from pigs with large litter size) to agarose gel electrophoresis, and then RNA concentrations were quantified using Agilent 2100 Bioanalyzer and NanoDrop spectrophotometer. Ribosomal RNA was removed using Ribo-Zero™ rRNA Removal Kit (Epicentre, Madison, WI, United States). Library preparation and sequencing analysis were described as previous study (Shen et al., 2019). The raw data files were cleaned up for quality control through Trim Galore and then mapped to the pig reference genome (Sscrofa 11.1) from UCSC (Kim et al., 2013). CuffLinks were used to linear transcripts assembly and abundance assessment, then we identified reads from fusion transcripts that did not match the linear RNA sequence. The lncRNAs ≤100 kb in length, with at least two support sequences, and sequences with no more than two mismatches were retained for further analysis.

Gene Ontology analysis (GO; http://www.geneontology.org) was used to explore potential functions of differentially expressed mRNA and lncRNA. GO terminology analyzes information about biological processes involving genes, and cellular components or metabolic pathways, and highlights molecular functions. In addition, we also performed the Kyoto Encyclopedia of Genes and Genomes pathway analysis (KEGG; http://www.kegg.jp), using DAVID to gain insight into the interaction and reaction network of differentially expressed lncRNA molecules. The -log₁₀ P value represents a significant enrichment of a given GO term or KEGG pathway between upregulated and downregulated entities.

Bnd lncRNA sequences, an interaction regulation network was constructed between mRNA-miRNA-lncRNA. The RNAhybrid (https://bibiserv.cebitec.uni-bielefeld.de/rnahybrid/) and TargetScan (http://www.targetscan.org/) were used to predict the interactions of miRNA-mRNA and miRNA-circRNA. The correlation between the expression levels of lncRNA and miRNA was calculated with SPSS Pearson correlation assay.

Invitrogen Trizol reagent (Thermo, MA, United States) was used to extract total RNA from porcine ovarian tissue, and PrimeScript™ RT kit (TaKaRa, Dalian, China) with gDNA eraser was used to reverse transcription of total RNA while removing genomic DNA. Real-time quantitative PCR (qPCR; 95°C 30 s, 95°C 5 s, 60°C 30 s for 40 cycles, following 70°C 10 min for elongation) of lnc RNA was performed using SYBR Green kit (TaKaRa, Dalian, China) in triplicate on a Bio-Rad CF96 system (Bio-Rad, United States), and the data was normalized using β-actin. The primers are listed in Supplementary Table S1.

The data showed mean ± SEM. P-values were calculated by a Student’s t-test. P < 0.05 was considered as significant difference.

In this study, we selected six ovarian samples to perform Ribo-Zero RNA-Seq from pig with small (SLS, 8.48 ± 0.53/litter) and large (LLS, 16.19 ± 0.43/litter) litter size after the fourth delivery, respectively (Supplementary Table S2). The significance of large litter size was higher than that of small litter size (P < 0.001). We used Pfam and Cpat databases to filter, and analyzed the protein coding ability predictions of new lncRNAs through CPC and CNCI tools, and identified 3,556 potential lncRNA transcripts (Figure 1A; Supplementary Table S3). As shown in Figure 1B, Pearson correlation analysis showed that there is a good correlation between SLS and LLS ovarian samples. We obtained 75–83 and 78–85 million unique mapped clean reads from the SLS and LLS sample libraries (Table 1), and found that 69.7% of the reads are located in the exon and sequence coding for aminoacids in protein (CDS) regions, while the distribution in the non-coding (intergenic, intronic, UTR, transcription initiation site (TSS), and transcription end site (TES) regions markedly reduced (30.3%; Figure 1C).

As shown in Figure 1D, the average expression level of the identified RNA including lncRNA and mRNA is 12.97 (FPKM), showing a lower expression level. In addition, we found that every chromosome could produce lncRNA, and the longer the chromosome, the more likely it is to produce lncRNA (Figure 1E).

In this study, 16,428 genes were detected (Supplementary Table S4), and 2,392 unigenes were significantly differently expressed at the LLS and SLS ovarian samples (P < 0.05; Figure 2A); all differentially expressed genes (DEGs) were shown in Supplementary Table S5. In order to verify the accuracy of RNA-seq, we randomly selected 8 genes (4 upregulated genes and 4 downregulated genes) for qPCR detection, and compared the qPCR data with the sequencing data (Table 2). It can be found that the qPCR results of the selected DEGs are basically similar to the RNA-seq results, indicating that the sequencing data has high accuracy.

As shown in Figure 2B, the cluster heat map of DEGs and the Pearson correlation in each sample were displayed. To further explore the potential functions of DEGs, we used volcano plots to compare their transcript expression abundance in SLS and LLS ovarian samples (Figure 2C).

To explore the possible biological functions of DEGs identified in SLS and LLS porcine ovarian tissues, GO analysis was performed (Figure 2D). The results showed that these DEGs were significantly enriched by 20 biological processes, including steroid biosynthesis process, lipid metabolism process, ovarian follicle development, granular cell differentiation, negative regulation of apoptosis process, and cell response to follicle stimulating hormone stimulation.

In order to better understand the molecular regulation of DEGs, we performed KEGG pathway analysis (Figure 2E). It showed that a total of 50 pathways are enriched, among which the most significant is the hippo signaling pathway (path: 04390), with 32 annotation genes. From the three enrichment pathways of steroid biosynthesis (path: 00100), ovarian steroid production (path: 04913) and steroid hormone biosynthesis (path: 00140), it could be seen that the synthesis and metabolism of steroids play an important role in the ovaries of pigs. In order to accurately obtain the key genes that may affect the pig litter size, we combined GO analysis and KEGG pathway analysis. The results showed that 15 genes may be candidate genes affecting piglet size traits. Among them, INHBA, BMP4 are related to ovarian follicle development, CAPN5 and FOXL2 are related to the regulation of granulosa cell differentiation; CYP21A2, STAR, HSD17B2, HSD17B8 and HSD17B12 are involved in steroid biosynthesis; In addition, the key genes TGFBR1 and TGFB3 in TGF signaling pathway are also significantly different (Figure 2F). It is worth noting that TGF signaling pathway is highly related to the regulation of reproductive performance.

From the fold change and the adjusted level of significance, a total of 302 differentially expressed lncRNAs were screened (P < 0.05; Figure 3A; Supplementary Table S6). We found that 206 lncRNAs were downregulated and 96 lncRNAs were upregulated (Figure 3B), and LLS group (FPKM = 346.38) had an obvious trend of high expression of lncRNA when comparing SLS group (FPKM = 256.71). It is worth mentioning that ALAS1 (FPKM = 32,089.67) and LOC100513133 (FPKM = 10,768) have the highest expression levels among all upregulated and downregulated lncRNAs, respectively. In order to better understand the expression changes and potential functions of lncRNA, a volcano plot is shown in Figure 3C.

LncRNA regulates the transcription and translation of coding genes in cis and trans regulation: when the role of lncRNA is limited to the same chromosome (for example, acting as a molecular scaffold), it performs cis regulation; when it affects gene expression on other chromosomes, trans regulation plays a role (for example, acting as a miRNA sponge). As shown in Figures 3D, E, the first 20 enriched GO terms and KEGG pathways are shown. Interestingly, metabolic processes in both GO and KEGG analyses had the most cis and trans’ target gene clusters, indicating that metabolic processes may play an important role in ovarian function.

In order to verify the accuracy of sequencing data, we randomly selected 12 genes (6 upregulated lncRNAs and 6 downregulated lncRNAs) for qPCR detection, and compared the qPCR data with the sequencing data (Figure 4; Supplementary Table S7). It can be found that the qPCR results of the selected differentially expressed lncRNAs are basically similar to the RNA-seq results, indicating that the lncRNA-Seq data are highly accuracy.

lncRNA can act as a molecular sponge of miRNA to sequester miRNA from its target gene, that is, miRNA is the common target of lncRNA and mRNA. We selected 9 miRNAs, 9 lncRNAs and 342 mRNAs related to ovarian development to construct an mRNA-miRNA-lncRNA (ceRNA) network (Supplementary Figure S1). For example, LOC100513133 and LOC102168075 have multiple binding sites for microRNAs related to ovarian development, such as miR-26b (Liu et al., 2016; Lin et al., 2012), let-7g (Zhou et al., 2016; Cao et al., 2015a) and miR-125b (Zhang D. et al., 2020; Yao et al., 2018; Yao et al., 2019). This ceRNA network may provide valuable information for ovarian function.