Combined analysis of mRNA–miRNA from testis tissue in Tibetan sheep with different FecB genotypes

Abstract Testis size is important for identifying breeding animals with adequate sperm production. The aim of this study was to survey the expression profile of mRNA and miRNA in testis tissue from rams carrying different FecB genotypes, including the wild-type and heterozygous genotypes in Tibetan sheep. Comparative transcriptome profiles for ovine testes were established for wild-type and heterozygote Tibetan sheep by next-generation sequencing. RNA-seq results identified 3,910 (2,034 up- and 1,876 downregulated) differentially expressed (DE) genes and 243 (158 up- and 85 downregulated) DE microRNAs (miRNAs) in wild-type vs heterozygote sheep, respectively. Combined analysis of mRNA-seq and miRNA-seq revealed that 20 miRNAs interacted with 48 true DE target genes in wild-type testes compared to heterozygous genotype testes. These results provide evidence for a functional series of genes operating in Tibetan sheep testis. In addition, quantitative real-time PCR analysis showed that the expression trends of randomly selected DE genes in testis tissues from different genotypes were consistent with high-throughput sequencing results.


Introduction
In mammals, testes are male reproductive organs that produce sperm, secrete androgenic hormones, and maintain male characteristics [1]. In production, we usually use a scrotal circumference as a proxy for testicular size, which is significantly and positively correlated with testicular size, thus indirectly indicating that the trait is easy to measure and has high heritability [2]. Studies have found positive correlations between testicular size and sperm motility, production, concentration, and ejaculate volume and a negative association with abnormal sperm percentage in livestock animals, including cattle, sheep, goat, and pig [3][4][5][6]. Testis size affects the age at puberty and the postpartum interval of offspring in females [7,8]. Therefore, the selection of elite sires based on testicular size is an effective method to improve reproductive performances. Thus, the testes are one of the important indicators for predicting male reproductive ability.
Testicles descend before puberty is reached and grow in size after puberty, increasing in volume by more than 500% [9]. This growth is due to lytic function, sperm production, and interstitial fluid and Sertoli cell fluid production. Sertoli cells and germ cells rapidly proliferate at birth and undergo a radical change in morphology and function at the onset of puberty [9]. This is due to a complex biological process, involving genetic and microRNA modulation, which is essential for testis development and sperm function in adult males. A gene called INHBA has been shown to affect testis morphogenesis, testicular cell proliferation, and testis weight in mice, and has a potential role in regulating testis size [10][11][12]. Among the regulators of spermatogenesis, small RNAs, including microRNAs (miRNAs) and PIWI-interacting RNAs (piRNAs), have recently attracted great attention [13]. miRNAs are small (∼22 nucleotides) endogenous RNAs that negatively regulate gene expression by targeting the 3′-untranslated region (3′-UTR) [14] or coding region [15] of mRNAs. Studies have shown that miRNAs play an important role in testicular development and spermatogenesis in mammals. New generation sequencing tools on mouse testis tissues at 8, 16, and 24 d found key miRNAs such as miR-34c, miR-221, miR-222, miR-20, and miR-106a [16]. Research has shown that miR-34c may inhibit the proliferation of male reproductive stem cells in dairy goats, and differentially expressed (DE) mRNAs and miRNAs have been identified in pig testis tissues at different developmental stages [17]. Expression profile and biological function studies on testicular tissues of yak, cattle, and cattle-yak have found that miRNAs are involved in apoptosis and are associated with the blockade of cattleyak spermatogenesis [18]. Another well-known gene that controls reproductive traits in sheep, BMPR-IB (FecB), has been studied by many researchers, focusing on ewes [19][20][21]. Ewes carrying the heterozygous genotype tend to produce more lambs, and, conversely, individuals carrying the wild type produce fewer lambs [22]. The relationship between ram fertility and FecB has not been elucidated, and we propose that a molecular mechanism exists for the higher fertility of rams carrying the heterozygotes compared to wild-type individuals. Although these studies describe mRNA and miRNA expression profiles of testicular tissues from different livestock animals, the genetic mechanisms of the testis tissue from different genotypes on FecB locus have not been investigated. Studies on the expression profiles of mRNA and miRNA in testis tissue of rams carrying wild type and heterozygotes in Tibetan sheep are lacking.
An indigenous sheep breed (Tibetan sheep) was selected for this study. This sheep is widely distributed as livestock in the Qinghai-Tibet Plateau (QTP) and adjacent areas in China. Tibetan sheep are well adapted to the harsh natural environment and are irreplaceable livestock in the high-altitude areas of QTP. They are the main source of meat in Tibet and have been found to have a low rate of double birth. Qiao et al. found FecB mutations in the ewes of Tibetan sheep [23,24]. Using a molecular marker of the FecB gene, the core population of multiparous Tibetan sheep was established, and the technical problems of low fecundity in ewe were solved. However, ram fertility is poorly researched. Mounting evidence indicates a series of genes that are critical for testis development, especially testis size and volume. AKT was shown to act as an anti-apoptosis gene in animal growth development and plays an important role in mouse testicular spermatogenesis [25]. Moreover, PLK1 was shown to be critical for testis development, and EHD4 is required for the survival of germ cells and attainment of normal prepubertal testicular size in mice [26]. EPAS1 is required for spermatogenesis in postnatal mouse testes. Postnatal EPAS1 ablation leads to male infertility, with reduced testis size and weight [27]. KITl is crucial for the testicular environment postnatally, and KIT is important for the proliferation and survival of spermatogonia [28].
Tibetan sheep reach sexual maturity at 8 months. In our preliminary research, we made a breakthrough in identifying the presence of FecB in Tibetan sheep and using this molecular marker for selection to expand the core population of Tibetan sheep. Because the relationship between ram fertility and FecB has not been elucidated in Tibetan sheep, we propose a potential molecular mechanism for the higher fertility of rams carrying the heterozygotes compared to the wild type. FecB has been reported more frequently in ewes and rarely in rams. This study is the first to study the relationship between FecB and testicular transcriptome in rams. In this study, we used Solexa deep sequencing to investigate mRNA, miRNA profiles, and construct mRNA-miRNA networks for testes from Tibetan sheep with different genotypes. This work provides a theoretical basis for the study of ram fertility and provides insights into testicular development in ruminants.

Animals and sample collection
Blood from 40 individual rams was randomly tested from the Tibetan sheep population for genotyping at the FecB gene locus (A746G) according to molecular detection methods for the major gene FecB in sheep [29]. In general, there are three genotypes at the FecB mutation locus: AA (also known as ++ or wild type), AG (heterozygous), and GG (homozygous). Based on the genotyping results of the FecB locus and the preciousness of heterozygote rams, which can be irreversibly damaged by slaughter, we selected six rams, including three wild types and three heterozygotes, for high-throughput sequencing. To ensure that the environmental conditions were similar throughout the experiment, all the lambs were housed in a well-ventilated room with controlled temperature and humidity. Testicular tissues were collected as described by Xu et al. [30]. All the animals and samples used in this study were collected according to the protocol approved by the institutional animal care and use committee of Qinghai University.
Ethical approval: The research related to animal use has been complied with all the relevant national regulations and institutional policies for the care and use of animals, and has been approved by the institutional animal care and use committee of Qinghai University.

RNA isolation and quantification
A total amount of 1 μg RNA per sample was used as input material for the RNA sample preparations. The total RNA and small RNA were isolated using TRIzol® reagent (TransGen Biotech, Beijing, China). RNA purity was checked using a NanoPhotometer spectrophotometer (IMPLEN, CA, USA). RNA degradation and contamination were monitored on 1% agarose gel. RNA concentration and integrity were measured using Qubit ® RNA Assay Kit in Qubit ® 2.0 Fluorometer (Life Technologies, CA, USA) and RNA Nano 6000Assy Kit of the Bioanalyzer 2100 system (Agilent Technologies, CA, USA), respectively. The experimental protocols were performed according to the manufacturer's technical instructions.

mRNA sequencing and analysis
cDNA libraries were prepared by first isolating total RNA from testis tissue. Total RNA was then reverse transcribed into cDNA, which was then amplified to create cDNA molecules. cDNA libraries were used for transcriptome sequencing. In this study, six RNA libraries were constructed and named TWT1, TWT2, TWT3, THT1, THT2, and THT3. Sequencing libraries were generated using the NEBNext Ultra RNA Library Prep Kit for Illumina (NEB, USA) following the manufacturer's recommendations. After cluster generation, the libraries were sequenced on an Illumina HiSeq Xten platform and 150 bp-paired end reads were generated. Clean reads were obtained by filtering out adapter sequences and removing low-quality reads from raw data. Simultaneously, Q30 and GC content of clean data were calculated. RNA-seq reads were aligned to the sheep reference genome (Oar_v3.1) using Tophat 2.0.10 [31] with default parameters. The gene expression level was calculated by reads per kilo-base per million reads after the read numbers mapped to each gene were counted by HTSeq v0.6.1 [32]. DE genes (DEGs) were examined using the DEGseq R package (1.12.0) [33]. A corrected P-value of 0.05 and an absolute value of log2 (fold change) of 1.5 were set as the thresholds for significant differential expression. Gene Ontology (GO) enrichment analysis of the DEGs was implemented by the GOseq R packages based on Wallenius non-central hyper-geometric distribution [34].

miRNA sequencing and analysis
miRNA libraries were prepared by first isolating total RNA from testis tissue. Total RNA was processed with a miRNAspecific reverse transcription step to create a miRNA library. miRNA libraries were used for miRNA expression profiling. Six small RNA libraries were constructed: SRWT1, SRWT2, SRWT3, SRHT1, SRHT2, and SRHT3. Sequencing libraries were generated using NEBNext ® Multiplex small RNA Library Prep Set for Illumina ® (NEB, USA) following the manufacturer's recommendations, and index codes were added to attribute sequences to each sample. Small RNAs ligated with 50 and 30 adapters were reverse transcribed and amplified. Library preparations were sequenced on the Illumina HiSeq Xten platform according to the manufacturer's instructions. After removing adapter sequences, reads containing poly-N, and low-quality reads, the clean reads were mapped to the Repeatmasker and Rfam database to remove tags originating from repeated sequences, rRNA, tRNA, snRNA, and snoRNA. The miRNA expression levels were estimated by transcript per million values. Differential expression analysis of miRNA was performed using the DEGseq R package [33]. P-value was adjusted using q-value, and q-value < 0.01 and |log2 (fold change)| > 1 were set as the threshold for significantly differential expression by default. The methods to construct mRNA-miRNA networks consist of three key points. First, data collection was performed by collecting DE miRNA and mRNA data from both the wildtype and heterozygote groups. Second, data analysis was performed using RNAhybrid tools [35] to compare the data and identify miRNA and mRNA sequences that are DE in the same tissues. Third, network construction was performed using Cytoscape software [36] to create a network diagram with DE miRNA and mRNA nodes.

Validation of mRNA by qRT-PCR
Total RNA was extracted using TRIzol ® reagent (TransGen Biotech, Beijing, China) as described by the manufacturer. After acquiring high-quality total RNA, mRNAs were reverse transcribed using the RevertAid™ First Strand cDNA Synthesis Kit (TransGen Biotech, Beijing, China). Quantitative real-time PCR (qRT-PCR) analyses on the mRNAs were performed using SYBR Green PCR Master Mix (ABI, USA, 4304437) in a Bio-rad CFX96 real-time system. Beta actin was used as an internal control for each mRNA. Primer sequences of all mRNAs were designed using primer 5 ( Table 1). For mRNA quantification, the reaction conditions were: 95°C for 10 min, followed by 40 cycles of 95°C for 10 s and 59°C for 50 s. Relative mRNA expression was evaluated using the 2 −△△CT method [37]. Differences were considered significant when P < 0.05 using one-way analysis of variance. All statistical analyses were performed using general linear models in R. Visualization of graphs was done with the ggplot package [38] in R.

mRNA expression analysis from testes in different genotypes
To obtain differential comparative transcriptome profiles for testes of different genotypes, six libraries, namely TWT1,   TWT2, TWT3, THT1, THT2, and THT3, were  showed multiple matches. Finally, 2,034 up-and 1,876 downregulated genes were identified in the wild-type group compared to the heterozygote group. These genes may be regulators of testicular growth and spermatogenesis. Sequencing analysis revealed that AKT, EHD4, TGFBR3, ODF3, INHA, VIM, COLIA1, and COLIA2 were DE in the testes of different genotypes. An extensive literature research revealed that these genes are closely related to testicular development, sperm morphology, and spermatogenesis. Therefore, we selected these genes for qPCR validation. To confirm the gene expression pattern identified by RNA sequencing, eight genes including four up-(AKT, EHD4, TGFBR3, ODF3) and four downregulated (INHA, VIM, COLIA1, COLIA2) genes were selected for real-time RT-PCR using testis samples from different genotypes. The results showed that the expression patterns of eight genes were consistent with those of the RNA sequencing results (Figure 1), implying that our RNA-seq results and analysis methods were highly reliable. GO enrichment analyses showed that the DEGs were significantly enriched for GO terms involving cell morphogenesis, tissue development, and spermatogenesis ( Figure 2).

miRNA expression analysis from testes of different genotypes
The Q30 values of the wild type and heterozygotes were more than 97.17%. The read length demonstrated a bimodal with 21-23 and 26-32 nt in six libraries (Figure 3), which was consistent with the classical fragments of miRNAs. These results showed that the sequencing quality was feasible. After obtaining clean reads, the clean reads were categorized using some database. Subsequently, unannotated clean reads of sRNA were mapped to the sheep reference genome (Ovis_aries3.1). Approximately 58.74-65.55% of clean reads were mapped to the sheep reference genome (Ovis_aries3.1) ( Table 3). The results indicated that the diversity and distribution of miRNAs in testes samples from different genotypes, including wild type and heterozygous, are relatively abundant. These results suggest that miRNAs are involved in testicular development. Finally, 243 DE miRNAs, including 158 up-and 85 downregulated miRNAs, were found in the wild type compared to the heterozygote group.

Integrated function analysis of genes and miRNA expression
Compared to the heterozygote group, 48 TDETGs were observed in the wild type, including 16 upregulated TDETGs corresponding to 8 downregulated miRNAs and 32 downregulated TDETGs corresponding to 12 upregulated miRNAs ( Table 4). miRNAs and protein-coding genes are involved in regulatory networks to regulate spermatogenesis and testes development. In terms of miRNAs, which can positively and negatively regulate target genes, both types of networks were constructed. Type I belongs to the regulatory network of upregulated mRNAs interacting with downregulated miRNAs in wild type compared to heterozygote group. Type II belongs to the regulatory network of downregulated mRNAs interacting with upregulated miRNAs in the wild type compared to the heterozygote group (Figure 4). A total of three miRNAs and ten genes were present in the type I network, whereas two miRNAs and ten genes were present in the type II network. We found that oar-miR-27a, miR-125b, miR-127, EHD4, TGFBR3, and PDGF were all active in the regulatory network.

Discussion
Xizang sheep, also known as Tibetan sheep, are widespread in QTP and adjacent areas. Using molecular markers of the FecB gene, the core population of multiparous Tibetan sheep was established and the technical problems of low fecundity in ewe were solved. However, ram fertility is poorly researched. Therefore, we investigated testicular genetics in wild-type rams and heterozygote rams to find effective molecular markers for measuring ram fertility. Sequencing results revealed a big difference in the wild type compared to the heterozygote group. In general, clean reads were assigned to three major regions, including exonic, intergenic, and intronic. Approximately 65.68-74.72% of clean reads were mapped to the sheep reference genome, of which 63.27-73.03% were uniquely mapped. The efficiency of mapping to the reference genome was low and we give several possible reasons for this. First, the Tibetan sheep selected are breed specific and may differ from the reference genome. Second, the sheep reference genome does not contain the Y chromosome, and we sequenced testicular tissue. A total of 2,034 Figure 1: Differentially expressed mRNAs were examined in different genotype testis tissue. One asterisk (P value <0.05) indicates a significant difference between two groups, two asterisks (P value <0.01) indicate a possible highly significant difference between two groups, and three asterisks (P value <0.001) indicate a highly significant difference between two groups. up-and 1,876 downregulated genes were identified in the wild type compared to the heterozygote group. These genes may be important for testicular growth and spermatogenesis.
The length distribution of small RNAs was concentrated in the two peaks of 19-24 and 26-32 nt. Notably, the clean reads of the wild-type group were mainly   miRNA based. Interestingly, the miRNA level gradually decreased, and the level of piRNA gradually increased in the heterozygote testes; these data were coherent with previous descriptions. These results indicate that the regulation of spermatogenesis is complex. Consistently, Hirano et al. [39] performed small RNA-seq analysis in the marmoset and found that piRNA is dominant in adult testes. Li et al. used high-throughput sequencing to compare miRNA expression in immature and mature pigs and reported that the proportion of miRNAs was higher before sexual maturity in large white pigs, but not in Luchuan pigs [40]. Lian et al. found a significant difference in the proportion of small RNAs before and after sexual maturation [41]. Similar to the results of this experiment, the proportion of miRNAs in the wild-type group was significantly higher than that in the heterozygote group. These results suggest that both miRNA and piRNA may be involved in testicular development during sexual maturity. TDETGs were used in this sequencing study to identify potential regulatory factors that influence testicular development. A total of 20 potential miRNAs were identified, including 8 down-and 12 upregulated miRNAs. These miRNAs include five known miRNAs (oar-miR-27a, oar-let-7b, oar-miR-125b, oar-miR-99a, and oar-miR-127), and newly discovered miRNAs. The functions of some known miRNAs have been reported. miR-27a is a differentiation marker for spermatogonial stem cells and is highly expressed in dairy goat testis cells [42]. In the mRNA-miRNA regulatory network constructed in this study, we screened miR-27a, suggesting that miR-27a may have potential regulatory effects on testicular development and differentiation of spermatogonial stem cells in sheep. Thus, these studies show that miRNA function is conserved in the vast majority of animals. In addition, studies in pigs have shown that miR-27a can be used as an internal reference gene for assessing the quality of cryopreserved sperm [43]. Studies have shown that miR-125b plays a crucial role in regulating zygote genome activation in oocytes and embryos [44]. miR-125b can target the degradation of the PAP gene, leading to an increase in testosterone concentration, a decrease in sperm count, and an increase in abnormal sperm ratio in male mice [45]. Related studies have shown that miR-125b is highly expressed in germ cells and directly regulates the secretion of testosterone by targeting PAP, increasing the mitochondrial DNA copy number in sperm cells, and ultimately affecting semen quality [45]. This indicates that miR-125b has a positive effect on reproductive performance and could be used as a drug and diagnostic target for male sterility. In the mRNA-miRNA regulatory network constructed in this study, we screened miR-125b, suggesting that miR-125b plays a key role in the development of sheep testes, especially hormone secretion and energy metabolism. miR-125b was found to play an important role in sheep reproduction. Although other miRNAs, such as miR-99 and miR-299, have been identified, their roles in reproduction are unknown [46,47]. Studies have shown that they can regulate stem cell self-renewal, ovarian cancer cells, and glioma, suggesting that their roles in reproduction should not be neglected.
Furthermore, we also found that the miRNA target genes were involved in reproduction. For example, the target gene EHD4 of miR-27a, which was originally described as a protein expressed in the extracellular matrix, has been shown to be a cytosolic and endosomal vesicle-associated protein [48]. George et al. [26] found that although Ehd4 −/− mice were fertile male, their 31 days postpartum testicular weight decreased by 50%, and apoptosis gradually increased. Other defects, including reduced seminiferous tubule diameter, abnormal regulation of spermatogenic epithelium, and abnormal head of slender sperm cells, are also present in Ehd4 −/− mice, which eventually result in lower sperm counts and lower fertility. In adult testes, EHD4 is highly expressed in primary spermatocytes and EHD4 deletion alters the levels of other EHD proteins in an age-dependent manner. The results suggest that EHD4 plays a role in the normal development of mitosis and late germ cells, and that EHD protein-mediated endocytic cycling is important for germ cell development and testicular function [26]. The target gene TGFBR3 is regulated by miR-127 and miR-27a. Beta glycans (TGFBR3) belong to one of the ligands of the transforming growth factor beta (TGFβ) superfamily. Sarraj et al. [49] found that Tgfbr3 −/− mice developed severe fine cords at 12.5-13.5 days after conception and their fetal stromal cell function was affected; their experimental data show that TGFBR3 is necessary for the normal formation of the spermatic cord, the development of fetal stromal cells, and the establishment of endocrine function. TGFBR3 was also screened in the network constructed in this study, suggesting that this gene is important for fine testicular formation, interstitial cell development, and endocrine function in sheep. The findings confirmed that miR-27a and miR-127 potentially control testis size. Testicular development is controlled by complex levels of gene regulatory proteins, growth factors, cell adhesion molecules, signaling molecules, and hormones that interact with each other through controlled relationships, usually in a short period of time.
Platelet-derived growth factor (PDGF) is a key regulator of connective tissue cells during embryogenesis and pathogenesis, including PDGFA, PDGFB, PDGFRA, and PDGFRB. The results of this study suggested that PDGFRB is a target of miR-27a. Mice with mutations in PDGF had normal testicular development before birth, normal fetal cytoplasmic cell numbers, and normal masculine functions. However, after birth, the gene mutation caused testicular shrinkage, loss of mesenchymal cells, decreased testosterone, and impaired spermatogenesis [50]. Another interesting finding is regulation of the target genes COLIA1 and COLIA2 by miR-125. There is evidence that COLIA1 and COLIA2 are involved in type A spermatogonia and play a role in mediating the isolation and migration of germ cells during spermatogenesis [51]. The aforementioned evidence revealed that these genes may play a key role in regulating testicular development.

Conclusions
In this study, 20 potential miRNAs, including 8 downand 12 upregulated miRNAs were identified by constructing mRNA-miRNA regulatory networks. The miRNA target genes, namely, TGFBR3, PDGFRB, Col1a2, and EHD4, are related to reproduction. Therefore, these miRNA molecules can be used as potential candidate small molecules to regulate the testicular size and provide a theoretical reference for the genetic regulation of testicular development in sheep.