Differences in gene expression and variable splicing events of ovaries between large and small litter size in Chinese Xiang pigs

Abstract

Background

Although lots of quantitative trait loci (QTLs) and genes present roles in litter size of some breeds, the information might not make it clear for the huge diversity of reproductive capability in pig breeds. To elucidate the inherent mechanisms of heterogeneity of reproductive capability in litter size of Xiang pig, we performed transcriptome analysis for the expression profile in ovaries using RNA-seq method.

Results

We identified 1,419 up-regulated and 1,376 down-regulated genes in Xiang pigs with large litter size. Among them, 1,010 differentially expressed genes (DEGs) were differently spliced between two groups with large or small litter sizes. Based on GO and KEGG analysis, numerous members of genes were gathered in ovarian steroidogenesis, steroid biosynthesis, oocyte maturation and reproduction processes.

Conclusions

Combined with gene biological function, twelve genes were found out that might be related with the reproductive capability of Xiang pig, of which, eleven genes were recognized as hub genes. These genes may play a role in promoting litter size by elevating steroid and peptide hormones supply through the ovary and facilitating the processes of ovulation and in vivo fertilization.

Summary

Based on analyzing of the transcriptome and alternative splicing events, twelve candidate genes related with fecundity and litter size trait were found out from the ovary of Xiang pig.

Background

In pig industry, the productive power of sows is one of most concerned economic traits in the world [1]. Reproductive traits are extremely intricate and influenced by multifactors originating from heredity and environment especially in litter size of pigs [2,3,4,5]. Lots of genes have functions on the reproduction capability [5, 6]. And number of litter size provides a direct effect on the economic benefits for pig farmer [7]. Ovaries are the main reproductive organs; they perform ovulation and show a direct influence on the efficiency of fecundity. Consequently, the different expression profiles of some important genes in ovary might devote to comprehend the diversity of litter size among breeds [8, 9].

Previous reports focus on the quantitative trait loci (QTL) together with intrinsic genes related with litter size of pigs, and the relationship between genes and the traits [10,11,12] is explored based on recent technology progress at molecular level [4, 13,14,15]. Series of genes connected with the fecundity of pig have been found, including FSH-β, ESR, OPN, MTNR1A, PRLR, GDF9 and BMPs members [16, 17]. Many genes and QTLs have been ascertained to have a linkage with the litter size trait of pig [18,19,20]. Differently expressed genes related with fecundity and litter size have been detected using transcriptome information from gonads in European pigs [4, 13] and Chinese local breeds [15, 21]. Nevertheless, knowledge of previous works aims at several limited pig breeds and they could not make it clear for the huge diversity of reproductive capability in Eurasian pig breeds.

Xiang pig is one of indigenous breeds in China originated from the southeast mountain environment of Guizhou province. It is featured by short stature, early maturity, excellent environmental adaptability as well as with nice meat quality [22, 23]. Furthermore, the populations among Xiang pig herds present great variation in litter size, ranging from 5 to 21 piglets, while most of sows gave 5–8 piglets per litter from the third to seven parities [24]. It was proposed that the cause might due to the specific regulation in the gene expression related with litter size trait. To screen pivotal genes related with litter size, the ovaries were sampled from two groups of Xiang pigs with large and small litter size. The expression profiles and alternative splicing events of transcripts were analyzed by RNA-seq method. The results will benefit to the interpretation of the molecular regulation manner on the diverse reproduction capability and litter size in pig breeds.

Materials and methods

Samples

A total of 40 Congjiang Xiang pigs were prepared from the farm of Dachang pig breeding, Congjiang, Guizhou, China, which born from sows ever giving birth of large litter size as XL group with the total number born (TNB) larger than twelve, or XS group with TNB less than eight. Using the same way from previous reports [15, 25, 26], we randomly chose fourteen pigs that the third estrous time was synchronous from XL group (n = 7) and XS group (n = 7) to sequence the transcriptome using Illumina next sequencing technology. All sampled pigs were 6 to 6.5 months old, weighing 37.50 ± 3.77 kg with five pairs of nipples. The estrous was first detected by B-ultrasound (KX5000V, XuZhou KaiXin Electronic Instrument Company, China) started from the onset of female standing reflex according to the method reported by Lopes et al. [27]. When the numbers of matured follicles were counted to be 4–8 with diameter larger than 6 mm on one ovary [28,29,30], both ovaries were picked out by standard surgical operation. The follicles numbers and size on two ovaries were directly measured before putting into liquid nitrogen for total RNA isolation. The sampled pigs were alive and kept feeding under routine method together with other pigs.

Library construction and sequencing

Based on protocol of Beijing Genomies Institute (BGI), Shenzhen, China, ovarian RNA was isolated with TRIzol method (Invitrogen) with the values of RNA integrity number (RIN) in the scope from 7.9 to 8.8 to prepare cDNA library. The libraries were sequenced using HiSeqTM 2000 platform (Illumina, USA) and generated 100 bp paired-end reads. The same RNA sample was determined in the RNA-seq and qRT-PCR tests. The cDNA libraries together with RNA sequencing were carried out as described previously [26].

Dataset analysis

The sequencing data from fourteen libraries were taken for analysis of the expression profile and the alternative splicing events of transcripts at estrus stage by RNA-seq method. The raw reads in fastq format were filtered to remove reads under low quality by program Trimmomatic v0.39 using the cutoffs as previous conditions [26], and the clean reads were aligned with the pig reference genomes (Ssc11.1) using HISAT2 software (v2.1.0). The mapped sam files were conversed and sorted into bam format by samtools (v1.1.0). The subread featurecounts software (v2.0.0) was chosen to count the reads amount, which were included in the regions of genes or exons. The expression level of gene was estimated by CPMs values (counts per million mapped reads).

The differently expressed genes with CPM value were calculated by using DESeq2 and edgeR, in which all of CPM values were added 0.001 for logarithm arithmetic. The minimum normalized CPM was 1.0, in which a gene would be eliminated if its CPM value of any sample was not lager than the threshold. The differently expressed genes with CPM values were computed using model featureCounts in subread program (v1.6.3). The threshold for differentially expressed genes (DEGs) was the gene with P < 0.001 (false discovery rate (FDR) < 0.005) and log2 ratio > 1 or < -1. rMATS (v4.0.2) was used to detect the differential AS events with an FDR ≤ 5 % and a |∆PSI | |IncLevelDifference| of ≥ 0.1. The value of PSI or ψ (percentage spliced in) value was calculated by rMATs according to the ratio of the long form on total form present to characterize inclusion of exon, differential splice-site choice, intron retention, etc. The DEGs and differentially spliced genes (DSGs) datasets were uploaded to the platform of KOBAS v3.0 (http://kobas.cbi.pku.edu.cn/kobas3) taking reference Sscrofa11.1 as background based on Ontology Consortium (http://geneontology.org/) and Kyoto Encyclopedia of Genes and Genomes (KEGG) for enrichment categories. The ovary DEGs and DSGs from Xiang pig were compared with previous reports [17, 31] via Venn program online. The DEGs and DSGs of Xiang pig related with reproduction were further analyzed. The gene list related with reproduction was first collected from the Gene annotation deposited in NCBI using reproduction as keyword. Then it was converted into the gene names and Ensembl IDs of pig via BioMart online (http://www.ensembl.org/biomart). Based on Venn analysis on the shared genes, the lists of DEGs and DSGs of Xiang pig related with reproduction were selected out.

RT-qPCR verification

To verify the DEGs and DSGs expression patterns deduced by transcriptome analysis, six samples in each group were used as the same aliquot of total RNA for RNA-seq detection. Seven DEGs (LDLR, SCARB1, HSD3B1, CYP11A1, AKR1C2, StAR and LRP8), two highly expressed genes (SERPINE2 and RARRES1) together with five types of AS events were chosen randomly to verify the RNA-seq analysis by quantitative real-time RT-PCR methods. Primers were devised by Primer3 online (http://primer3.ut.ee/). The PCR reaction conditions and proportion were the same as our previous work [32] with each primer concentration of 10 pM/µL, taking GAPDH and β-actin genes as internal controls. Based on dissociation curve analysis for PCR products, the amplification efficiency was controlled within range of 100 ± 10 %. The relative expression level of target gene utilized the method of 2(-ΔΔCt) as reported by Livak et al. [33]. The different level of gene expression between two groups was tested by software SPSS (v21.0) taking the P < 0.05 as threshold of significant difference. The results were presented as mean ± standard deviation. The presence of five types of splicing evens were determined by RT-PCR method using pairs of primers outside of the AS region. And the DSGs levels were further qualitatively detected by qRT-PCR method using one primer span both junction ends of AS event. The nucleotide sequences of primers were listed in Table S1.

Results

Illumina sequencing and assembly

After quality control, the cDNA libraries of each ovary sample generated 49 ~ 66 million clean reads with 100 base pairs (bp) in length. The fourteen samples showed similar matching results, with mapping ranges of 96.17 ~ 99.96 % onto the genome (Ssc11.1) and 69 ~ 75 % being unique matches (Table 1). The results indicated that all of fourteen libraries present high-quality, together with high percentage of coverage throughout pig reference genome. It enabled us to analyze the transcriptomes profiles from ovary at estrus stage of Xiang pigs with large or small numbers of litter size.

Table 1 Overview of RNA-Seq data

Differential gene expression analysis between XS and XL groups

Based on mapping data to the reference of pig genome, we obtained 17,089 and 15,928 genes from XS and XL groups, respectively. Normalized CPM data in Table S2 were inputted for principal component analysis (PCA). It appeared that each point in seven samples of the first group could gather and separate from the points in another group according to PC1, which accounted for percentage of 45.1 % of total variation in the dataset (Fig. 1A). It meant that the distance within samples in the same group was much close to each other than that in another group. After removing the noncoding RNA and pseudogene transcripts and those genes with CPM < 1.0 in each sample, the sequencing data of 16,476 genes could be used for following analysis. Of them, 15,389 genes were expressed in both groups; 984 genes were specifically expressed in the libraries of XS group, while 103 genes were determined only from XL group (Fig. 1B). Most of the genes specially expressed in XS or XL group presented in low or very low CPM values. GO analysis indicated that the especially expressed genes in XS libraries enriched mainly in the cellular process, regulation of biological process and metabolic process. However, the genes only expressed in XL libraries had no statistically significant GO terms (Table S3).

Fig. 1

Profiles of gene expression in Xiang pig ovaries between XL and XS groups. A: PCA cluster of the gene expression profile of fourteen libraries. B: Venn diagram of expression genes in two groups. Blue color represented genes only expressed in XS group, orange colors showed genes only expressed in XL group, and the intersection is the common genes in both groups

In total, 2,795 genes were differently expressed between two groups after intersection of results from both softwares of edgeR and DESeq2. Compared with the expression in XS group, 1,376 genes were down-regulated and 1,419 genes were up-regulated in XL group (Table S4). The scope of log2FC values was varied from − 8.75 to 9.31 in DEGs. The numbers of genes, with more than four times of difference between two groups, accounted for 26.11 % of the total DEGs. Approximately 18 and 22 % of the DEGs expressed less than 100 normalized CPM in XS and XL respectively, in which 9.26 % of these genes overlapped between two groups. Moreover, a volcano diagram was plotted based on DEGs data (Fig. 2). The expression levels and numbers of DEGs in XL group were more than that in XS group, as displayed in the heat map of Fig. 3. We compared the top ten genes highly expressed in XS group with that in XL groups (Table 2). The expression level in XL group ranged from 3545.875 ~ 15422.571 CPM (logFC from 0.148 to 6.42), which was decreased to 2594.786 ~ 6588.433 CPM (logFC from − 2.238 to 2.307) in XS group. Of those, the expression levels of five genes (StAR, ATP6, COX3, COX1, SELENOP) increased and two genes (MACF1, HSPG2) decreased in XL group.

Fig. 2

Volcano plots of the differently expressed genes (XL vs. XS). Each point in the figures represented one gene. Red points represented up-regulated genes, green points denoted down-regulated genes. Black points were genes without significant difference

Fig. 3

Heatmap showed the most significantly up-regulated and down-regulated genes between XL and XS groups

Table 2 The top 10 genes at highly expressed level in XS and XL groups

More strict thresholds were used to screen out the significant DEGs, which included FDR less than 0.001, the value of normalized CPM larger than 100 and the absolute value of log2FC larger than four. We identified 37 significant differently expressed genes from the total DEGs between two groups, such as StAR, TIMP1, NR4A1, PTX3, CYP11A1, PTGFR, OVGP1, SERPINE1, CLDN11, and MSMO1 (Table 3).

Table 3 The 37 DEGs with great fold changes between XS and XL groups

Compare of DEGs between Xiang and Yorkshire pig

Compared with previous report from Yorkshire pig [17], 71 DEGs were found from ovaries of both pig breeds. Of them, 64 genes were up-regulated and 6 genes were down-regulated in the groups with high litter size (Table 4). And six genes, COX3, STAR, CYTB, CYP11A1, MSMB and SCARB1, were ranked in the highest expression level between two pig breeds.

Table 4 Compare of DEGs between Xiang and Yorkshire pigs (reported by Zhang et al. [17])

Detection of DSGs and AS events

Five basic types of AS events were classified, including A5SS (alternative 5′splice site), A3SS (alternative 3′splice site), SE (skipped exon), RI (retained intron) and MXE (mutually exclusive exon). The results showed in Table 5. We identified 63,837/64,075 AS events including splice junctions only (JC) and splice junctions and reads on target (JC + ROT) in 11,414/11,468 genes from XS and XL datasets. Thus, approximately 69 % of 16,476 expressed protein-coding genes were subject to alternative splicing. The numbers of AS events were from 1 ~ 28 events (JC or JC + ROT) in a gene. The highest number of alternative splicing event was found out from gene ARHGEF7 (ENSSSCG00000009551) with 28 events. SE was the most prevalent AS event, followed by MXE, and RI. Compared with other AS forms, the high frequency of SE indicated that the manner of skipped exon significantly might impact transcription and resulted in various isoforms during gene transcription. A total of 4,009 / 7,441 (JC / (JC + ROT)) significant differential alternative splicing events were identified from 2,763 / 3,936 genes (Table S5). Of these, 542 differently spliced genes (DSG) also exhibited differently expressing (Table S4). Compared with XL group, the number of up-regulated AS events was significantly more than that of down-regulated events in XS group (Table 5).

Table 5 The types of AS events in ovaries of Xiang pig with small and large litter size

Compare of DSGs between Xiang and Yorkshire pig

Based on Venn results of the shared genes between Xiang pig and Large White sows [31], 1,597 DSGs from ovaries of Xiang pigs (Table S6) were also detected as much as 2,236 events by using single-molecule long-read sequencing (SMRT) in 39 tissues of Large White sows [31].

Changes in expression and AS of the reproduction genes between XS and XL

To understand the possible effects of the variations in expression and AS types of the reproduction genes on the litter size trait, crowds of major reproductive genes were picked out and explored the difference in expression pattern and AS of these genes in two groups. The RNA-seq data from 162 genes involved in reproduction processes were listed in Table S7. We found that 22 genes, such as ESR1, ESR2, GNRH1, FSHR, AR, GDF5, IRS1, CCND2, and so on, were down-regulated, and 33 genes, such as PTGS2, LIF, ECM1, BMPR1B, GPX3, C4BPA, MMP19, MMP25, STAT3, ect, were up-regulated in the XL group. The alternative splicing analysis indicated that 24 / 86 (JC / JC + ROT) significant differential AS events were presented in 42 reproduction related genes. However, of these, only 11 DEGs were also differential splicing, including AR, G6PD, ESR1, ECM1, BMPR1B, HEXB, STAT3, DNMT1, C4BPA, MMP23B, and LIN9 (Table 6).

Table 6 The eleven genes related with reproduction harboring differential expression and splicing

Gene ontology and KEGG analysis

To explain the biological effects of DEGs, we carried out GO and KEGG enrichment analysis (Table S3). For the up-regulated genes between XL and XS groups, 59 significantly enriched KEGG pathways were identified (corrected P-value < 0.05), including metabolism pathways (carbon metabolism, citrate cycle (TCA cycle), amino acid metabolism, glycerophospholipid metabolism, cholesterol metabolism etc.), oxidative phosphorylation, illness pathways (rheumatoid arthritis, Parkinson disease, colorectal cancer, type I diabetes mellitus, insulin resistance, hypertrophic cardiomyopathy etc.), physiological process related paths (renin-angiotensin system, complement and coagulation cascades, aldosterone synthesis and secretion, bile secretion, cardiac muscle contraction, adrenergic signaling in cardiomyocytes), immune paths (cell adhesion molecules, antigen processing and presentation), five signaling pathways (MAPK signaling pathway, NOD-like receptor signaling pathway, PI3K-Akt signaling pathway, TNF signaling pathway, PPAR signaling pathway). Notably, three pathways, ovarian steroidogenesis, steroid biosynthesis, FoxO signaling pathway, were included in the regulation of steroid hormone and ovary function. The 904 GO terms were enriched in kinds of physiology process. Significantly, 12 GO terms were gathered in reproductive process, female pregnancy, mammary gland epithelium development and proliferation, placenta development, blastocyst development and embryo development.

Meanwhile, the down-regulated genes enriched in two KEGG pathways, which were oocyte meiosis and progesterone-mediated oocyte maturation (corrected P-value < 0.05). And the GO terms gathered in meiosis I and meiosis I cell cycle processes. It illustrated that the DEGs between XL and XS groups were important in both function and development of reproductive system and hence were probably to contribute to litter size between two Xiang pig groups.

Very few KEGG paths and GO terms in both the up-regulated and down-regulated DSGs were significant at the level of correct P value less than 0.05 (Table S3). It was found that the up-regulated DSGs were enriched in 19 KEGG pathways and 181 GO terms if the threshold was reduced to the P value less than 0.05, including autophagy, endocytosis, and lysosome, phosphatidylinositol signaling system, metabolisms such as protein processing in endoplasmic reticulum, lysine degradation, fructose and mannose metabolism, fatty acid metabolism, ubiquitin mediated proteolysis etc. And the down-regulated DSGs could be enriched in 35 KEGG pathways and 219 GO terms including metabolism, growth cell process and reproduction etc. Of those, the GO term of utero embryonic development, and two KEGG pathways (oocyte meiosis, and progesterone-mediated oocyte maturation) were related with reproduction. These results indicated that the DSGs also participated the regulation processes on oocyte maturation and ovary function in pig.

Candidate genes related with litter size trait in Xiang pig

We performed the Venn analysis on the differently expressed genes, differently spliced genes, the higher expressed genes (CPM ≥ 100), the top expressed genes and the top differently expressed genes between XL and XS groups. Then, we selected these genes that overlapped in three or more datasets for the next GO and KEGG pathway analysis. According to the reported functions in any similar paths that relate to ovarian steroidogenesis, fecundity, pregnant or embryo development, we identified 12 candidate genes (StAR, DHRS4, RLN2, PTX3, HSD3B, MSMO1, SCARB1, COX1, COX3, SELENOP, CYP11A1, and NR4A1) having a linkage with pig reproduction capability and litter size. Of them, eleven candidate genes were identified to be hub genes connected with 12 to 35 genes based on the network relationship of DEGs by using string online platform (Table 7, Fig. S1).

Table 7 The detection of eleven hub genes and the protein-protein network

Tests and verification

The trends of expression of nine genes via qRT-PCR detection were positively related to that from RNA-seq data (Fig. 4). And all of five types of AS events were detected out from the transcripts of ovaries (Fig. S2). It indicated that the analysis based on RNA-seq data was precise and effective.

Fig. 4

Validation of DEGs by qRT-PCR method. The trend was similar in fold change (XL/XS) from RNA-seq and the ratio of expression levels in groups by qRT-PCR method with the linear correlation coefficient of 0.9532, P = 0.037 based on two-tail T-test

Discussion

In this study, we analyzed the estrus ovarian transcriptome and AS of Xiang pig using Illumina next generation sequencing technology. We detected 16,476 genes that expressed in ovaries from libraries. Of these, there were 15,389 genes expressed in common between two groups (Fig. 1B). The expressed amounts of these genes were much diverse, with CPM values changed from 1 to more than twenty thousand (Table S2). Great amounts of genes expressed specifically in XL or XS group detected to be at low or very low level. Further, the top ten genes highly expressed in XL group were compared with that in XS groups (Table 2), and found that the expression levels of some genes increased in XL group, such as StAR, SELENOP, COX3, COX1, and so on. The genes with expression level ranked top ten occupied 5.04 ~ 6.96 % of the total expression values in XS and XL groups, respectively. It indicated that these high expressed genes were considerable necessary in the function and the development of ovary. For instance, the transport of cholesterol into mitochondria dependents on the effect of steroidogenic acute regulatory protein (StAR), which accelerates the transform of cholesterol into the inner membrane of mitochondrion to trigger steroidogenesis reaction [34]. In mitochondrion, cholesterol is changed into pregnenolone by the cytochrome P450 side-chain cleavage enzyme. And pregnenolone is further transformed into progesterone or dehydroepiandrosterone, two hormones essential for endometrial receptivity, embryo implantation, and the successful establishment of pregnancy [35]. Previous works indicate that selenium (Se) regulates the growth of granulosa cells together with 17-estradiol synthesis in ovary [36]. Another report showed that both of Se and its selenoprotein (SELENOP), as antioxidants, promote the growth and proliferation of granulosa cells [37].

Previous study for goat ovaries suggests that some special differently expressed genes based on RNA-Seq data may improve litter size [9]. In present work, we identified 2,795 DEGs, including 37 most differently expressed genes between XS and XL groups (Table 3). It indicated that these genes might be very important for the litter size of pig. Results from enrichment analysis indicated that the up- and down-regulated DEGs were clustered in many GO terms and pathways, including metabolism, growth, development, and reproduction. And the effects of the top thirty-seven DEGs between XL and XS groups mainly included ovarian steroidogenesis, metabolic pathways, oxidation-reduction process, negative regulation of endopeptidase activity. Compared with Yorkshire sow [17], 71 DEGs (Table 4) and both pathways (steroid biosynthesis and ovarian steroidogenesis) were shared with Xiang pigs (Table S3). Moreover, 15 genes were up-regulated in the group with high litter size of Xiang and Yorkshire ovaries (Table S3). Some of them are reported to have a pivotal role in ovary. For example, Cyp11a1 protein catalyzes the transformation from cholesterol to pregnenolone in mitochondrion in luteal cells. Both of StAR and Cyp11A1 genes are taken as two marker of corpus luteum in mice [38]. And the StAR gene governs the rate-limiting step in steroidogenesis described above [34]. The oxidase HSD3B promotes the oxidation of both delta 5-ene-3-beta-hydroxy steroid and the ketosteroids. The enzyme, 3-beta-HSD, is necessary in the anabolism for all kinds of steroid hormones [39]. As the receptor of HDL, SCARB1 participates in the optional absorption of cholesteryl ether together with the transport outside of HDL-dependent cholesterol, and even accelerates the flow of esterified or free cholesterol on cell surface together with modified lipoproteins [40]. MSMO1 takes part in the first reaction to remove the two C-4 methyl from molecule 4, 4-dimethylzymosterol [41].

Furthermore, the other DEGs listed in Table 3 were reported to have a close connection with ovary function, such as NR4A1, DHRS4, and PTX3. Nuclear receptor subfamily 4 group A member 1 (NR4A1) is an orphan receptor in nucleus, which regulates the transcription of androgen biosynthesis and the expression of paracrine factor insulin-like 3 (INSL3) in thecal cell of ovary. Androgens together with another hormone control the follicle growth in ovary [42]. NR4A1 distributes in many cells of ovary including theca cell, luteal cell, and granulosa cell in human. Furthermore, NR4A1 in Leydig cell is reported to affect the expression of gene StAR in mouse [43]. The dehydrogenase/reductase SDR family member 4 (DHRS4) gene, also known as NADPH-dependent retinol dehydrogenase/reductase (NRDR) gene, is a tetrameric protein that is pivotal to the biosynthesis of steroid hormone. DHRS4 functions as NADPH-dependent 3-ketosteroid reductase to produce the 3β-hydroxysteroids from 3-keto-C19/C21-steroids. Types of 3β-hydroxysteroids are reported to transmit signal and participate various physiology functions, such as binding to estrogen receptor β (ER-β) in nucleus and changing the development of prostate [44]. PTX3 is specifically expressed by cumulus cells around oocyte, and mediates the effect of LH or hCG in preovulatory follicle. PTX3 actively participates in the organization of the hyaluronan-rich provisional matrix required for successful fertilization. And PTX3 is taken as a biomarker of oocyte quality and has a role in oocyte maturity and female fertility based on gene deficiency mice [45].

These genes mentioned above, including Cyp11A1, StAR, HSD3B, SCARB1, MSMO1, NR4A1, PTX3, DHRS4 and so on (Tables 2 and 3), were all up-regulated in the ovaries of Xiang pigs with large litter size. It suggested that these genes might play important roles in promoting litter size by increasing the level of steroid and peptide hormones supply through the ovary and facilitating the oocyte ovulation and in vivo fertilization.

It is interested that many alternative splicing events from DEGs were detected based on comparison between XS and XL groups. About 69 % of all expressed genes contained AS events in both of XL and XS groups, which is much near to the AS rates in human [46]. Total of 1,597 nonredundant genes with differentially splicing in Xiang pig also detected many isoforms from tissues of Large White pigs [31] (Table S6). In DSGs of Xiang pig, skipped exon was the most prevalent AS events. The rates of AS events in XL group were not as high as that in XS group. AS is the main reason leading to change the different transcripts together with proteome varieties [47]. Numbers of reports indicate that alternative splicing interferes the functions of animal genes and alters the receptor structure especially in the processes of development and growth [48]. Lots of hereditary disease appear strong relationship with high frequency of alternative splicing in genes [49]. It was deduced herein that the high percentage of AS in pigs of XS group might cause the decrease of fecundity. Moreover, we found that 542 DEGs were differently spliced at AS levels between two groups (Table S4). The DEGs presented different and special patterns of splicing and events. Many tops differently expressed genes, such as StAR, MSMO1, SCARB1 and PDZK1IP1 showed high percentages of differently alternated splicing events (Table 2, Table S5). However, there were 3,693 genes only undergoing differently AS events between XL and XS groups. Therefore, the expression patterns and AS events of 162 genes related with reproductive processes were explored profoundly from the RNA-seq datasets. And eleven genes were found to be differently expressed and differently spliced in ovarian samples between XS and XL groups (Table S7). In addition, 31 reproductive genes only underwent differential splicing, such as CYP19A1 and FMR1. Estrogens are essential for animal fertility, which are catalyzed by aromatase enzyme coded by gene Cyp19a1 [50]. The gene encoding aromatase of mammals contains two promoters, including gonad specific and brain specific promoters. It exists 10 promoters at tissue-specific manner with the first exon to be chosen differently in diverse tissue cells. In kinds of promoters of CYP19A1 gene, the most vigorous one is promoter II (PII), which drives the transcription of aromatase gene in ovary [50]. The studies from rat find that the expression of aromatase transcription present the diverse and active regulation [51]. Gene FMR1 (FMRP translational regulator 1; FMRP: fragile X mental retardation protein) is composed of seventeen exons occupying about 38 kb in genome [52]. The gene undergoes extensive AS which changes the retain of four exons, 12, 14, 15 and 17, producing various FMR1 transcript isoforms, and some of FMRP isoforms have been reported in several species [52, 53]. In rat follicles, the FMR1 gene was transcribed in granulosa cell, theca cell and germ cell. FMR1 mRNA is much less in pre-ovulatory follicles than that in both preantral and antral follicles. FMRP content raises in the development process of follicles, and could be detected more than four bands by Western blotting method [53]. In present work, the expression of CYP19A1 mRNA isoform with MXE event was significantly down-regulated and FMR1 mRNA isoforms with A3SS event were significantly down-regulated in XL group (Table S5). These results indicated that the changes of gene expression between groups with large or small litter size were moderated at many ways and the splicing variants were highly controlled.

Finally, combined with DEGs, DSGs and the higher expressed genes via Venn analysis, we identified 12 candidate genes related with litter size in Xiang pig, including StAR, DHRS4, RLN2, PTX3, HSD3B, MSMO1, SCARB1, COX1, COX3, SELENOP, CYP11A1, and NR4A1. And eleven of them were identified to be hub genes based on the network relationship of DEGs (Table 7, Fig. S1). It indicated that these candidate genes might play important roles in the regulation of reproduction. The effects of StAR, DHRS4, RLN2, HSD3B, MSMO1, SCARB1, COX1, COX3, CYP11A1 and NR4A1 involve in ovarian hormone biosynthesis and regulation, which could regulate the processes of ovary development, oocytes mature and the quality of embryos. The functions of SELENOP are related to follicular growth and oocyte maturation [37]. Gene PTX3 increases the progress of oocyte ovulation and fertilization in vivo [45]. These genes expressed at a high level in ovaries of XL group and may accelerate ovarian hormone biosynthesis and the quality of oocytes. They might connect to the higher reproduction performance in Xiang pig. There are eleven candidate genes related with the litter size of Yorkshire breed reported from transcriptome analysis [17]. Interestingly, five of them (STAR, COX3, HSD3B, SCARB, CYP11A1) were also found to be related with the litter size trait in Xiang pig in this work. But the other candidates from Yorkshire breed are much different from Xiang pigs. It suggested that trait of litter size in pig breeds shared candidate genes and further controlled by diverse genes because of their various genetic background.

Conclusions

In short, this study showed a transcriptome pattern and AS profiles at estrus stage of ovaries from Xiang pigs. We identified 1,419 genes that showed up-regulated and 1,376 genes that appeared to be down-regulated in the large litter size samples. And it was found that approximately 69 % of expressed genes harbored AS treatment. Of 542 DEGs also exhibited differential alternative splicing. Based on previous finding on those genes, total of 12 candidate genes were found to be corresponding to the reproduction capability and litter size in Xiang pig. These genes play important roles in promoting litter size by increasing steroid and peptide hormones supply by ovary and facilitating the oocyte release and in vivo fertilization.