Identification and characterization of mRNAs and lncRNAs in the uterus of polytocous and monotocous Small Tail Han sheep (Ovis aries)

Samples

All experiments were performed following the relevant guidelines and regulations set by the Ministry of Agriculture of the People’s Republic of China (No. IASCAAS-AE-03).

Based on the TaqMan assay using the FecB mutation probe, a total of 12 pluriparous ewes with FecB ++ genotypes were selected from nucleus herds of Small Tail Han sheep in the southwest region of Shandong Province, China. These ewes were all approximately three years old and weighed 63 kg. Ewes were divided into two groups: a polytocous group (PG, n = 6, litter size and ovulation number ≥ 2) and a monotocous group (MG, n = 6, litter size and ovulation number = 1) based on ovulation number and three lambing records. All animals had free access to water and food under natural lighting.

All ewes were treated with synchronous estrus. A vaginal sponge was first implanted for 12 days. Estrus was then tested by a ram each day. These ewes were then divided into follicular and luteal phase groups. Uteri from six ewes (three polytocous ewes and three monotocous ewes) were collected between 45 and 48 h (follicular phase; PF and MF, respectively). Additionally, uteri from the other six ewes (three polytocous ewes and three monotocous ewes) were collected on a ninth day (luteal phase; PL and ML, respectively). All samples were immediately stored at −80 °C for total RNA extraction.

RNA extraction, library construction, and RNA-seq

Total RNA was extracted from the uterus of 12 ewes using TRIzol (Invitrogen, Carlsbad, CA, USA) according to the manufacturer’s instruction. RNA concentration was measured using the Kaiao K5500 spectrophotometer (Beijing Kaiao Technology Development Co., Ltd, Beijing, China). RNA integrity was assessed using the RNA Nano 6000 Assay Kit of the Agilent Bioanalyzer 2100 System (Agilent Technologies, Santa Clara, CA, USA).

The rRNA was depleted from 3 µg of total RNA using Ribo-Zero™ Gold Kits (Epicentre, Madison, WI, USA). Sequencing libraries of the 12 samples (PF, n = 3; MF, n = 3; PL, n = 3; ML, n = 3) were generated using NEB Next Ultra Directional RNA LibraryPrep Kit for Illumina (NEB, Ipswich, MA, USA) according to the manufacturer’s instructions, and index codes were used to label the sequences of each sample. After cluster generation, the library preparations were sequenced on an Illumina Hiseq platform (Illumina, San Diego, CA, USA). Raw data of the performed RNA-seq have been recorded in the SRA public database (Accession number: SRP173986).

Reference genome mapping and transcriptome assembly

Raw data in fastq format were processed through in-hoseperl scripts. In this step, clean reads were obtained by removing reads with adapter contamination, reads that contained poly-N, and low-quality reads from raw data. Simultaneously, the Q20, Q30, and GC contents of the clean data were calculated. All downstream analysis was based on high-quality clean data. HiSAT2 (Pertea et al., 2016) was used to align clean reads of each sample to the sheep reference genome Oar_v3.1. StringTie (Pertea et al., 2015) was used for transcriptome assembly and reconstruction. Thus, known lncRNA and mRNA transcripts were identified, and the position of transcripts was obtained.

Identification of potential lncRNA candidates

LncRNAs were identified using the following workflow. (1) Transcripts >200-nt long with >2 exons were obtained. (2) Transcripts with a coverage less than 5 in all samples were removed. (3) The different classes of class_code annotated by “u”, “i”, and “x” were retained, which corresponded to lincRNAs, intronic lncRNAs, and anti-sense lncRNAs, respectively. (4) Used Gffcompare to compare with annotation files to screen out known mRNAs and other non-coding RNAs (e.g., rRNAs, tRNAs, snoRNAs, snRNAs). Transcripts without coding potential, as predicted by CNCI, CPC, PFAM, and CPAT, were candidate lncRNAs.

Differentially expressed gene analysis

The FPKM was used to normalize the expression levels of lncRNAs and mRNAs, which eliminated the effect of sequencing depth, gene length, and sample difference on gene expression levels (Trapnell et al., 2010). For experiments with three biological replicates, the differentially expressed lncRNAs and mRNAs were identified using the R package DEseq (Anders & Huber , 2010) after the negative binomial distribution. For biological replicates, lncRNAs and mRNAs with P <0.05 and Fold change >1.5 were considered differentially expressed between the polytocous and monotocous ewes of different estrus cycles.

GO and KEGG pathway enrichment analysis of differentially expressed genes

Gene Ontology (GO) enrichment analysis of differentially expressed genes or lncRNA target genes was implemented by the GOseq R package, in which gene length bias was corrected (Young et al., 2010). GO classifies functions into three groups: cellular components, molecular functions, and biological processes. The KEGG biological pathways database (http://www.genome.jp) is a central public database for understanding high-level functions and regulatory network research. Enrichment analysis was performed on each Pathway in KEGG using a hypergeometric test. The calculated P value and 0.05 being defined as the significant threshold, the genes were screened and enriched for the pathways.Next, the significance of the pathway enrichment analysis was corrected by FDR, and the corrected P-value (Q-value) was obtained. Differentially expressed genes were further studied using the GO and KEGG databases to study the functions of the genes and identify the pathways in which they participate. If a P value was ≤0.05, enrichment was considered significant.

Prediction and functional analysis of differentially expressed lncRNA target genes

The primary role of lncRNAs, which are a type of noncoding RNA, is to regulate their target genes by cis-regulating nearby protein-coding genes and trans-regulating distal protein-coding genes. Protein-coding genes located 50-kb upstream and downstream of a lncRNA in a genome are cis-target genes, whereas protein-coding genes with a correlation coefficient >0.9 with a lncRNA were trans-target genes (Fatica & Bozzoni, 2013). GO annotation and KEGG pathway enrichment analyses were performed on the obtained lncRNA target genes to identify the biological processes and signaling pathways enriched in the lncRNAs. Then, we predicted the functions of lncRNAs.

LncRNA–mRNA network construction

To further explore the interactions between the lncRNAs, target genes and differentially expressed genes in sheep reproduction. Based on the targeting relationship between mRNA and lncRNA, we screened networks related to uterine function and reproduction with reference to their GO and KEGG enrichment terms and keywords for classification. Visualization of gene interactions is achieved through an open software platform called Cytoscape (V3.1.1) (Saito et al., 2012).

Gene expression validation by quantitative real-time PCR

We used qRT-PCR to verify gene expression levels. We used approximately 0.1 µg of each RNA sample and reverse transcribed into cDNA using RT reagent. Real-time PCR was performed at 95 °C for 10 min, followed by 95 °C for 15 s, 60 °C for 60 s for 45 cycles, and 72 °C for 30 s. qPCR was performed on the LightCycler 480II (Roche, Basel, Sweden) using SYBR Green Real-time PCR Master Mix (TOYOBOCO, LTD, Osaka, Japan). β-Actin was used as an internal reference to normalize target gene expression. All experiments were performed in triplicate. mRNA and lncRNA primers are shown in Table 1.

Statistical analyses

All data were expressed as “means ± SD” At the time of comparison, a Student’s t-test was performed, and P < 0.05 was considered statistically significant.

Identification and characterization of lncRNA

A total of 1,489,144,532 clean single-end reads were obtained by sequencing all 12 libraries. Each library single-end reads of obtained were above 11.9 million. Reads were then aligned onto the Ovis aries reference genome using HiSAT2. Approximately 90% to 94% of the reads were successfully aligned to the Ovis aries reference genome (Table 2).

A total of 25,104 lncRNAs were identified in uterine tissues of the 12 ewes using four programs: CNCI, CPC, PFAM, and CAPT (Fig. 1A); 20,908 mRNAs and 16,016 novel transcripts were identified. Many lncRNAs have only two or three exons, whereas mRNAs contain a wide range of exons from two to thirty (Fig. 1B). Overall, the distribution of lncRNAs and protein-coding gene lengths were consistent, and the transcript lengths of lncRNAs were longer than those of mRNAs (Fig. 1C).

Gene expression levels and differential expression analysis

The Fig. 1D box plot shows that lncRNA transcript expression levels were all lower than those of mRNAs in the uterus of both polytocous and monotocous Small Tail Han ewes. Based on a fold change of >1.5 and a false discovery rate of <0.05, in the follicular phase, 242 lncRNA transcripts were up-regulated and 231 were down-regulated, and 33 mRNA transcripts were up-regulated and 133 were down-regulated in the PG ewes (Fig. 2A, Tables S1 and S2). Moreover, in the luteal phase, 330 lncRNA transcripts were up-regulated and 637 were down-regulated, and 359 mRNA transcripts were up-regulated and 146 were down-regulated in the PG ewes (Fig. 2B, Tables S3 and S4).

GO annotation and KEGG enrichment analysis of differentially expressed genes

A total of 166 follicular phase and 505 luteal phase differentially expressed mRNAs were analyzed by GO analysis (P < 0.05); all were categorized into biological processes, cellular components, and molecular function (Tables S5 and S6). Reproductive process was enriched in the top 10 terms, which indicates that some genes regulate sheep reproductive traits in the follicular phase (Fig. 3A) and luteal phase (Fig. 3B) of polytocous and monotocous sheep.

KEGG pathway analysis revealed that, in the follicular phase, the differentially expressed mRNAs were assigned to 127 pathways (Table S7); they participate in ovarian steroidogenesis, steroid hormone biosynthesis, oxytocin signaling pathway, cell cycle, retinol metabolism and other critical regulatory processes (Fig. 4A). In the luteal phase, the differentially expressed mRNAs were distributed in 240 pathways (Table S8); they participate in the prolactin signaling pathway, steroid biosynthesis, steroid hormone biosynthesis, ovarian steroidogenesis, TGF-β signaling pathway, estrogen signaling pathway and other important regulatory processes (Fig. 4B). These results indicated that the differentially expressed genes in the ovine uterus of PG and MG ewes are involved in reproduction performance.

LncRNA target genes and functional analysis

The genes transcribed within a 50-kb window upstream or downstream of the lncRNAs were considered cis-target genes, while targets in trans were predicted via calculating the expressed correlation (correlation coefficient ≥ 0.9) with lncRNAs.

In the follicular phase, three annotated lncRNAs corresponded to 11 target genes, and 470 novel lncRNAs corresponded to 1,578 target genes. GO terms of these lncRNA targets were enriched for several processes, such as reproductive process and cell proliferation (Fig. 3C, Table S9). KEGG analysis of these differentially expressed lncRNA target genes revealed that they were enriched in the foxo signaling pathway, TGF-β signaling pathway, estrogen signaling pathway, cell cycle, steroid hormone biosynthesis, ovarian steroidogenesis and retinol metabolism (Fig. 4C, Table S10).

In the luteal phase, 19 annotated lncRNAs corresponded to 157 target genes, and 948 novel lncRNAs corresponded to 3108 target genes. GO terms of these lncRNA targets were enriched for several processes, such as reproductive process, cell proliferation and metabolic process (Fig. 3D, Table S11). KEGG analysis of these differentially expressed target genes of lncRNAs revealed that they were enriched in the foxo signaling pathway, estrogen signaling pathway, TGF-β signaling pathway, ovarian steroidogenesis and VEGF signaling pathway (Fig. 4D, Table S12).

LncRNA–mRNA co-expression network analysis

In the follicular phase, lncRNA-mRNA co-expression networks were constructed using 38 differentially expressed lncRNAs and 46 target genes involved in reproductive-related pathways. As shown in Fig. 5, some differentially expressed lncRNAs were at the center of the network, such as MSTRG.134747, MSTRG.82376, MSTRG.135103, MSTRG.82370, and MSTRG.163615 (Table S13). In the luteal phase, lncRNA-mRNA co-expression networks were constructed using 174 differentially expressed lncRNAs and 164 target genes involved in reproductive-related pathways. As shown in Fig. 6, some differentially expressed lncRNAs were at the central positions of the network, such as MSTRG.134747, MSTRG.134749, MSTRG.134746, MSTRG.134751, MSTRG.170669, MSTRG.153309, MSTRG.116096, MSTRG.134753, and MSTRG.135103 (Table S14). The network model showed that each lncRNA was co-expressed with polytocous genes, which indicates mutual regulation of lncRNA and mRNA in reproduction.

RNA-Seq data validation by real-time PCR

To further validate the sequencing data, we selected eight involved in reproductive-related differentially expressed mRNAs and eight differentially expressed lncRNAs of targeted reproductive-related genes, and determined their expression levels by qRT-PCR (Fig. 7). The expression of each mRNA or lncRNA in the PG and MG ewes in the follicular and luteal phases (Figs. 7A, 7B) were consistent with those obtained by sequencing.