The exploration of miRNAs and mRNA profiles revealed the molecular mechanisms of cattle-yak male infertility

Cattle-yak, the first-generation offspring of cattle and yak, inherited many excellent characteristics from their parents. However, F1 male hybrid infertility restricts the utilization of heterosis greatly. In this study, we first compared the testicular tissue histological characteristics of three cattle, three yaks, and three cattle-yak. Then we explored the miRNA profiles and the target functions of nine samples with RNA-seq technology. We further analyzed the function of DE gene sets of mRNA profiles identified previously with GSEA. Testicular histology indicated that the seminiferous tubules became vacuolated and few active germ cells can be seen. RNA-seq results showed 47 up-regulated and 34 down-regulated, 16 up-regulated and 21 down-regulated miRNAs in cattle and yaks compared with cattle-yak, respectively. From the intersection of DE miRNAs, we identified that bta-miR-7 in cattle-yak is down-regulated. Target prediction indicated that the filtered genes especially MYRFL, FANCA, INSL3, USP9X, and SHF of bta-miR-7 may play crucial roles in the reproductive process. With further network analysis and GSEA, we screened such hub genes and function terms, we also found some DE gene sets that enriched in ATP binding, DNA binding, and reproduction processes. We concluded that bta-miR-7 may play an important role in influencing fecundity. Our study provides new insights for explaining the molecular mechanism of cattle-yak infertility.

Introduction

Studies have revealed that only about 2% of genome transcripts have protein-coding ability (1). However, as the most numerous components of genomic transcriptional products, the non-coding RNAs have also been confirmed that play significant roles in regulating gene expression at different levels such as transcription, RNA processing, translation, and epigenetic regulation (2–4). The non-coding RNAs can be classified into small and long ones based on length (5). microRNAs (miRNAs) are endogenous small non-coding RNA molecules 19–25 nucleotides in size that regulate posttranscriptional silencing of target genes, and single miRNAs can influence the expression patterns of many genes by targeting hundreds of mRNAs (6). In mammals, many studies have focused on miRNA and mRNA and increasing evidence had demonstrated that miRNAs can participate in multiple biological processes such as proliferation (7), differentiation (8, 9), and apoptosis (10). Moreover, miRNAs can act as epigenetic modulators to affect the protein levels of the target mRNAs without modifying the gene sequences, and the expression of miRNAs is also regulated by epigenetic machinery such as DNA methylation, and RNA modification (11). In recent years, there have also been a large number of related studies on the regulation of azoospermia and male sterility by miRNAs.

Cattle-yak, an excellent hybrid of cattle (Bos taurus) and yak (Bos grunniens), exhibits better economic traits than its parents such as meat quality, cold resistance, and growth performance. For example, the phospholipids containing long-chain polyunsaturated fatty acids have been confirmed more abundant in the muscles of plateau cattle than plain cattle, and are considered beneficial to human health (12–14). But the Haldane's rule demonstrates that the heterogametic sex of first-generation hybrids is infertile. Male infertility of F1 generation of cattle-yak resulting from spermatogenic arrest has greatly restricted the effective utilization of heterosis from the hybrids and the cultivation of new breeds (15, 16). Male cattle-yak have normal external genitalia, libido, and mating behavior, but testicular histology revealed that its seminiferous tubule walls were thinner and spermatogenesis was arrested at the stage of spermatogonia differentiation, which was hardly ever able to proceed further than meiosis I (17, 18). In epigenetics, unlike yaks, the expression levels of acetyl-histone H3 Lys9 (ACK9) in cattle-yak were abnormal and the 5-methylcytosine level was higher (19). Studies also observed that the androgens levels were decreased and androgen receptors led to higher expression of 3β-hydroxysteroid dehydrogenase(3βHSD) in the Leydig cells of cattle-yak testes (20). As a member of deleted azoospermia family genes, the Boule gene has also been proven to play important role in regulating primordial germ-cell formation (21). Entirely different from cattle, the expression level of the Boule gene was lower because the long CpG island in the promoter was hypermethylated in the testes of cattle-yak hybrids, and the methylation level was higher in the core promoter than outside (22). With the development of molecular biology in recent decades, such genes and transcription factors associated with male infertility had been identified. For example, DNA-binding protein PRDM9 is involved in double-strand breaks initiating meiotic recombination, and the position of the double-strand breaks (DSBs) is closely related to the fertility of male hybrids (23). Researchers had found that the expression level of bovine meiosis defective 1 (bMei1) participated in spermatogenesis was significantly decreased and the methylation level of promoter and gene body was significantly increased in the testis of cattle–yak compared with cattle. On the contrary, with the treatment of methyltransferase inhibitor 5-aza-2′ deoxycytidine (5-Aza-CdR), the activation of bMei1 would recover, further illustrating the association between bMei1 expression level and methylation (24). Besides, as a subfamily of the Argonaute protein, PIWI-like protein 1 (PIWIL1) plays a crucial role during germ-line development, stem cell maintenance, spermatogenesis, and transposon regulation, and the expression levels of PIWIL1 mRNA and PIWI protein in the testes of cattle-yak are significantly lower than that of cattle or yaks because of hypermethylation (25–29). In addition, Rad51 and meiotic-specific DMC1 recombinases can act synergistically in the process of homologous recombination that is important for spermatogenesis. A previous study had proved the levels of Rad51 and DMC1 are extremely decreased in the male cattle-yak testis with a corresponding higher incidence of germ cell apoptosis (30). As a family of transcription factors, heat shock factors can regulate the expression levels of heat shock proteins (HSP) and perform crucial roles during gametogenesis and development in physiological conditions, infertility of cattle-yak may be caused by the upregulation of HSP53 (31–33). With the development of high-throughput sequencing technology, many studies had preliminary performed comparative transcriptome profiles analysis in testes tissues and epididymis of cattle-yak and yak (34–37). However, few studies explored the miRNA profiles and the function of their target genes of cattle, yak, and cattle-yak comprehensively.

In our previous study, we compared the lncRNA and mRNA profiles of the testicular tissues for the differential expression analysis of cattle, yaks, and cattle-yak (38). In this study, we also used RNA-seq technology to conduct the comparison of three varieties of miRNAs profiles and their target genes. we also performed GSEA of DE gene sets and network analysis to identify the differences in gene expression patterns. We hope to find some significant biological markers and function terms associated with cattle-yak male sterility and have a deeper understanding of the molecular mechanism of this problem.

Results

Testicular tissue histological characteristics of nine samples

Compared with cattle and yaks, the testicular tissue sections of cattle-yak showed that seminiferous tubules were highly vacuolated and the walls became thinner. The reverse of cattle and yak (Figures 1A,B), the area of the longitudinal section of the testis was expanded and there are few germ cells or spermatocytes can be seen in the seminiferous tubules of cattle-yak (Figure 1C).

FIGURE 1

Figure 1. The histological differences among the testis of cattle, yak, and cattle-yak. The cattle longitudinal section of testicular tissues (A) and yak longitudinal section of testicular tissues (B) showed that the relatively thick seminiferous tubules filled with spermatogonia, spermatocytes, and spermatids. On the contrary, seminiferous tubules of cattle-yak (C) were highly vacuolated and there are few spermatids can be seen.

Quality control of sequencing data

After removing the 3' adaptors with cutadapt and filtering the low-quality reads with trimmomatic, clean reads were finally converted to FASTA format. There were 3,192,697–6,617,422 unique clean reads generated by filtering low-quality and covering the repeated reads. The GC percentage of these reads was from 46.31 to 49.29%, and all the Q30 values of these reads were surpass 96% (Table 1).

TABLE 1

Table 1. The statistics results of clean data.

Reads mapping and filtration

Blastn was used to align the reads to the Rfam database (rRNA, snRNA, snoRNA, tRNA) and then filtered the mapped reads with gapopen =0, evalue <0.01, and mismatch ≤ 1 (Supplementary Table 1). Bowtie was used to map the exon and intron sequence with mismatch ≤ 1 and then filtered the reads that can be mapped to the exon sequence and can not be mapped to the intron sequence. Clean reads were mapped to the bovine reference genome (ARS-UCD1.3) with mismatch ≤ 1 (Table 2).

TABLE 2

Table 2. The result of reads mapped to the reference genome.

The miRNAs family and expression level analysis

After filtering the reads that did not belong to miRNAs, the mirDeep2 was used to predict the novel miRNAs and secondary structure of miRNAs. After quantification of miRNAs, we analyzed the miRNAs family to further understand the function of miRNAs to their target genes (Figure 2). Then these counts were normalized with RPM. The RPM distribution of nine samples can be seen as the RPM boxplot (Figure 3A). At the same time, we analyzed the principal component of nine samples. These nine samples can be divided into H, M, and P groups (Figure 3B). With the miRNAs expression levels of each sample, we further analyzed the correlation between three samples in one group and the distance among three groups, the samples of each group were highly relevant and the distance among three groups was obvious (Figure 4).

FIGURE 2

Figure 2. miRNAs family count top 20 and total reads top reads. (A) Top 20 counts of mature miRNAs from one miRNAs family. (B) Top 20 counts of reads of these miRNAs. NA represents that the reads are consistent with the characteristics of miRNAs but didn't belong to any miRNA families we identified.

FIGURE 3

Figure 3. The RPM distribution and principal component of nine samples. (A) The RPM boxplot of nine samples, the x-axis represents each sample, and the y-axis represents log₁₀ RPM. (B) Principal component analysis of nine samples, the x-axis represents the PCA1, and the y-axis represents the new dimension PCA2, each color and shape represents each variety.

FIGURE 4

Figure 4. The correlation and distance among nine samples. (A) Sample correlation heatmap of nine samples, each small square represents the correlation between each sample of horizontal and vertical coordinates. The larger the correlation coefficient, the closer the color is to red. (B) Sample to sample distance heatmap of nine samples, each small square represents the distance between each sample of horizontal and vertical coordinates. The farther the samples are, the lighter the color.

Differential expression analysis of miRNAs among nine samples

The R package edgeR was used to analyze the differential miRNAs expression among nine samples with CPM >1, p-value < 0.05, |log₂FoldChange| >1. And then the clustering analysis and differentially expression of H vs. P and M vs. P were then diagramed with a heatmap (Figure 5). There were 81(47 up-regulated and 34 down-regulated) and 37(16 up-regulated and 21 down-regulated) differentially expressed miRNAs in H vs. P and M vs. P, respectively (Figure 6). The list of differentially expressed miRNAs can be seen in Supplementary Table 2. Then, we consolidate the common up-regulated and down-regulated miRNAs of cattle-yak in the merged analysis of H vs. P and M vs. P (Figure 7).

FIGURE 5

Figure 5. Heatmap of clustering analysis and differentially expression of H vs. P and M vs. P. (A) Heatmap of differentially expression analysis of H vs. P. (B) Heatmap of differentially expression analysis of M vs. P. Red represents the high expression of miRNA of each comparison group, and green represents the low ones.

FIGURE 6

Figure 6. Volcano plots of differentially expressed miRNAs of H vs. P and M vs. P. (A) Volcano plot of differentially expressed miRNAs between cattle (H) and cattle-yak (P). (B) Volcano plot of differentially expressed miRNAs between yak (M) and cattle-yak (P). Each point represents each DE miRNA, red represents up-regulated and green represents down-regulated.

FIGURE 7

Figure 7. Venn plot of cattle-yak miRNAs in the merged analysis between H vs. P and M vs. P. (A) Venn plot of up-regulated miRNAs of cattle-yak in the merged analysis between H vs. P and M vs. P. (B) Venn plot of down-regulated miRNAs of cattle-yak in the merged analysis between H vs. P and M vs. P.

Target prediction, functional enrichment, and network analysis

With the list of the differentially expressed miRNAs of H vs. P and M vs. P, we used Miranda software to predict their target genes. In the H vs. P group, there were 6 up-regulated and 1 down-regulated miRNA in cattle that can target prediction successfully (Table 3).

TABLE 3

Table 3. The list of differentially expressed miRNAs and their target genes of H vs. P.

Besides, M vs. P group filtered 2 up-regulated and 3 down-regulated miRNAs in yak can target the protein-coding genes and further regulate the associated signaling pathways (Table 4). We also observed that bta-miR-7 is commonly down-regulated in cattle-yak between H vs. P and M vs. P, the predicted target genes MYRFL, FANCA, INSL3, USP9X, and SHF may play significant roles in cattle-yak male-sterility. With the predicted genes, we carried out GO and KEGG enrichment analyses to annotate their function (Supplementary Tables 3–6). We selected the top 10 GO terms of molecular function, biological process, cellular component in GO enrichment, and the up-regulated genes of cattle most engaged in the negative regulation of ubiquitin-dependent protein catabolic process (GO:2000059), positive regulation of cell migration (GO:0030335), thiol–dependent ubiquitinyl hydrolase activity (GO:0036459), and positive regulation of mitochondrial DNA replication (GO:0090297). Besides, KEGG enrichment analysis demonstrated that these target genes are most enriched in the cell cycle (bta04110), RNA transport(bta03013), and PI3K–Akt signaling pathway(bta04151), Focal adhesion (bta04510) pathways. On the other hand, the up-regulated genes in yak are most enriched in the DNA-binding transcription factor activity (GO:0003700), positive regulation of DNA demethylation (GO:1901537), protein deubiquitination (GO:0016579), zinc ion binding (GO:0008270) terms and relaxin signaling pathway (bta04926), RNA polymerase (bta03020) pathways, such the up-regulated GO terms and pathways in cattle and yak may positively participate in the reproductive processes (Figures 8A–D). With the analysis of miR-7 target genes function, we found these genes enriched in hormone activity, signaling receptor binding and activator activity, molecular function regulator process, Fanconi anemia pathway, and relaxing signaling pathway that is associated with hormone secretion and signal transduction (Figures 8E,F). At the same time, we performed Gene Set Enrichment Analysis (GSEA) to interpret the enrichment results of DEG sets. We found that the gene sets of cattle enriched in the ATP binding, reproduction, and developmental process involved in reproduction processes and focal adhesion, olfactory transduction, and MAPK signaling pathway in the H vs. P group. Likewise, the gene sets of yaks are most involved in the ATP binding, DNA binding, and reproduction processes, as well as the cytokine receptor interaction and MAPK signaling pathways (Figure 9). We also found that the positive regulation of cell migration (GO:0030335), negative regulation of cell adhesion (GO:0007162), SMAD2-SMAD3 protein complex processes, PIK3R1 gene, SMAD2 gene, relaxin signaling pathway (bta-04926) and cellular senescence (bta04218) were located in the center of GO and KEGG regulatory network of H vs. P DE miRNAs target genes (Figures 10A–C). Similarly, in the M vs. P DE target genes GO and KEGG regulatory network, membrane raft (GO:0045121) process, and focal adhesion (bta04510) pathways were in the central location. The ITGB3 gene may be the hub genes among these pathways (Figures 10D–F). All the enrichment results demonstrated that the male sterility of cattle-yak may be associated with a lack of regulation of such processes.

TABLE 4

Table 4. The list of differentially expressed miRNAs and their target genes of M vs. P.

FIGURE 8

Figure 8. Target genes GO and KEGG enrichment top bar plot and point plot of different groups. (A) Top 30 bar plot of enriched GO terms of DE miRNAs target genes in the H vs. P group. (B) Top 20 points plot of KEGG pathway of DE miRNAs target genes in the H vs. P group. (C) Top 30 bar plot of enriched GO terms of DE miRNAs target genes in the M vs. P group. (D) Top 20 points plot of KEGG pathway of DE miRNAs target genes in the M vs. P group. (E) Point plot of miR-7 target genes enriched GO terms. (F) Point plot of miR-7 target genes enriched pathways.

FIGURE 9

Figure 9. GSEA enrichment plot of differentially expressed genes in H vs. P and M vs. P. (A) The first five go terms of DEG in the GSEA-GO analysis of H vs. P. (B) The first five KEGG pathways of EDG in the GSE-KEGG analysis of H vs. P. (C) The first five go terms of DEG in the GSE-GO analysis of M vs. P. (D) The first five KEGG pathways of EDG in the GSEA-KEGG analysis of M vs. P. One curve for each color represents each term or pathway.

FIGURE 10

Figure 10. The network of GO enriched terms and pathways in the H vs. P and M vs. P groups. (A) The GO to GO network of up-regulated terms in cattle between the comparison of cattle and cattle-yak. (B) The pathway enrichment top network between the comparison of cattle and yak. (C) The pathway to pathway top network between of up-regulated in cattle between he comparison of cattle and cattle-yak. (D) The GO to GO network of up-regulated terms in yak between the comparison of yak and cattle-yak. (E) The pathway enrichment top network between the comparison of yak and cattle-yak. (F) The pathway to pathway top network between of up-regulated in yak between he comparison of yak and cattle-yak.

qPCR validation

To confirm the relative expression levels of DE miRNAs and their target genes were consistent with the RNA-seq results, we selected the DE miRNAs bta-miR-449a, the target genes FANCA, MYRFL, and SHF of bta-miR-7 to perform the qPCR validation. We found that the miR-449a relative expression of cattle-yak was significantly lower than cattle and yak. And the target gene expression levels of FANCA, MYRFL, and SHF were significantly up-regulated in cattle-yak, which is consistent with the lower expression of cattle-yak and negative regulation of bta-miR-7 (Figure 11).

FIGURE 11

Figure 11. The relative expression of selected DE miRNAs and target genes in cattle, yak, and cattle-yak groups with qPCR validation. (A) The relative expression of miR-449a in cattle, yak, and cattle-yak. (B) The target gene FANCA relative expression of bta-miR-7 in cattle, yak, and cattle-yak. (C) The target gene MYRFL relative expression of bta-miR-7 in cattle, yak, and cattle-yak. (D) The target gene SHF relative expression of bta-miR-7 in cattle, yak, and cattle-yak.

Discussion

As important post-transcriptional non-protein-coding regulators, miRNAs engaged in many cellular and biological processes such as proliferation, apoptosis, differentiation, and tumorigenesis. miRNAs most mediate negative post-transcriptional regulation of their targets by base pairing to sequence motifs in the 3′ UTR of mRNAs as epigenetic modifiers (39–41). It is elucidated that individual miRNAs can control the expression of more than one target mRNAs and that each mRNA may be regulated by multiple miRNAs (42). Over the past few decades, the studies that focused on the regulation of reproductive processes by miRNAs is gradually increasing. miRNAs-34/449 family members have been verified that play crucial roles in spermatogenesis, regulating spermatozoa maturation and functionality, the dysregulation of ciliogenesis in the efferent ductules is significantly impaired, leading to sperm aggregation and agglutination as well as to defective reabsorption of the seminiferous tubular fluids. Then the efferent ductules were obstructed and tubular fluids were accumulated, resulting in high hydrostatic pressure into the testis, leading to testicular dysfunction and spermatogenic failure (43). Besides, a previous study used the small RNA sequencing technology to reveal the differentially expressed miRNAs in High- and Low-motile sperm populations of cryopreserved semen and found that the DE miRNAs miR-17-5p, miR-26a-5p, miR-486-5p, miR-122-5p, miR-184, and miR-20a-5p may relate to reproductive processes (44). Thus, it can be seen that miRNAs could play great roles in regulating reproductive traits. Although there were systemically studies that related to the problem of male sterility of cattle-yak, the comprehensive exploration of comparing the miRNAs profiles with cattle-yak and their parents' generations were rarely can be seen.

In this study, we compared the DE miRNAs and DEGs among cattle, yak, and cattle-yak groups. In the H vs. P groups, we found a total of 47 up-regulated and 34 down-regulated miRNAs, and there were 7 known miRNAs (6 up-regulated and 1 down-regulated) that can be mapped to their targets. The target genes ESPL1, KDM6A, and BCORL1 may be engaged in reproductive processes. In many previous studies, ESPL1 (extra spindle pole bodies like 1) have been found involved in cellular mitosis and cell cycle process (45–47). KDM6A (lysine demethylase 6A) belongs to the KDM6 family of histone H3 lysine 27 (H3K27) demethylases, which is an X-linked protein that contains a catalytic Jumonji C (JmjC) domain can facilitate the removal of the methyl group on di- and trimethylated H3K27 (H3K27me2/3), it can not only manifest demethylase activity but also regulate gene expression and enhancer activation independently of its demethylase activity. Many studies had demonstrated the function of KDM6A in differentiation, development, regeneration, and tumor suppression (48–50). BCORL1 was also proved that can be a contributor to spermatogenesis (51). Furthermore, we merged the common differentially expressed miRNAs between the two comparison groups and found that the known miRNAs bta-miR-7 were down-regulated in cattle-yak. With the analysis of target genes, we found that the FANCA gene is associated with the processes of DNA damage repair and the variation of FANCA would induce premature ovarian insufficiency and further cause female fertility (52, 53). In this study, we quantified the expression levels of the FANCA gene in cattle, yak, and cattle-yak and detected higher relative expression levels in cattle-yak. We concluded that the FANCA overexpression would increase the DNA stability excessively and interfere with the processes of spermatocyte meiosis.

Furthermore, we performed GO, KEGG, and GSEA analyses to annotate the function of the DE miRNAs target genes and DEGs. With GO and KGEE enrichment annotation, we found that the up-regulated miRNAs target genes in cattle and yak were most engaged in negative regulation of ubiquitin-dependent protein catabolic process, positive regulation of cell migration, positive regulation of mitochondrial DNA replication, positive regulation of DNA demethylation, protein deubiquitination processes and cell cycle (bta04110), PI3K–Akt signaling pathway(bta04151), RNA polymerase (bta03020), Focal adhesion (bta04510), and relaxin signaling pathway (bta04926). In our previous study, we have described the function of lncRNA target genes preliminarily (38). In this work, we further used gene set enrichment analysis (GSEA) to expatiate the function of DE gene sets. We performed GSE GO and GSE KEGG to annotate the function of DEGs, respectively, and found that the DE gene sets were enriched in ATP binding, the developmental process in reproduction, DNA binding, transporter activity, extracellular region processes, MAPK signaling pathway, pathways in cancer, and focal adhesion pathways. Studies had proved that both noncatalytic ATP binding and stable ADP binding can mediate chemo-mechanical transduction in axonemal dynein (54). ATP levels were highly associated with the levels of GSK3α sperm hexokinase activity (55). As an efflux pump of the ABCG subfamily of the ATP-binding cassette (ABC) transporters, ATP-binding Cassette Transporter G2(ABCG2) was also detected in the blood-testis barrier and epididymal, which was proven to play a great role in epididymal sperm maturation (56–58). Mitogen-activated protein kinases (MAPKs) are important signal transducing enzymes that are involved in many facets of cellular regulation such as gene expression, cell proliferation, cell differentiation, programmed cell death, and cell motility (59). MAPK signaling pathway can regulate dynamics of tight junctions and adherens junctions, proliferation, and meiosis of germ cells, and proliferation and lactate production of Sertoli cells (60, 61). Focal adhesion kinase is a widely expressed protein-tyrosine that plays an important role in intracellular signal transduction pathways triggered in response to cell interactions with the extracellular matrix (62). Focal adhesion pathways of focal adhesion kinase-mediated have been implicated in a diverse array of cellular processes, such as cell migration, growth factor signaling, cell cycle progression, and cell survival (63). Such the processes and pathways in cattle and yak up-regulated target genes and gene sets may positively engage in the sperm cell proliferation and seminal tubules differentiation.

The spermatogenesis contains spermatocytogenesis (mitosis), meiosis, and spermiogenesis (64). With a cycle of several mitotic divisions, spermatogonia and primary spermatocytes were produced and stem cells were renewed. In the stage of meiosis, the genetic material was duplicated and exchanged and then four haploid round spermatids were produced. Lastly, the haploid spermatids were differentiated, and mature spermatozoa were released into the lumen of seminiferous tubules. However, the male hybrid of cattle and yak cannot finish the above processes. Histological examination of seminiferous tubules revealed that gonocytes and spermatocytes were established normally, however, spermatogenesis was arrested at the meiosis phase began 10 months after birth in the hybrids (65). With the testicular tissue histological characteristics, we found that the seminiferous tubules were highly vacuolated and there were few spermatids to be seen. Previous studies also demonstrated that although cattle-yak has the same number of chromosomes (2n = 60) as that of cattle and yak, most of the primary spermatocytes of cattle-yak performed morphologically abnormal and have no XY bivalents, which hindered the subsequent processes of meiosis and spermiogenesis (18). The genes that participated in meiosis such as MEI1, SYCP3, and DAZL have been detected significantly down-regulated in adult cattle-yak testes compared to those in yak testes (24, 66, 67). Besides, the target gene USP9X identified in this study is essential for proper spermatogenesis. USP9X-null spermatogenic cells underwent apoptotic cell death at the early spermatocyte stage and then caused subsequent aberrant spermiogenesis, which resulted in complete infertility of USP9X conditional knockout male mice (68). The comprehensive regulation of such miRNAs target genes and DE genes would disturb the spermatocyte's meiosis process and further influence the reproductive traits of cattle-yak.

Materials and methods

Animals and testis sample collection

Three male cattle (H1, H2, and H3), three male cattle-yak (P1, P2, and P3), and three male yaks (M1, M2, and M3) of good health and consistent diet levels, with the same growth conditions aged at 18 months were randomly selected on a pasture of Hongyuan country, Sichuan province of China. They were then divided into three groups (H, P, and M). Testicular tissue samples were collected and stored at −80°C until total RNA was extracted when these animals were slaughtered.

Histological analysis of testicular tissues

Testicular tissues collected from nine samples were immersed in 4% phosphate-buffered formalin overnight. Then the fixed paraffin-embedded tissues were sectioned with 4μm thickness. The sections have further been stained with hematoxylin and eosin (HE) and toasted in an oven set to 70°C.

RNA extraction, library preparation, and sequencing

Total RNA was extracted from nine testicular tissue samples with a Total RNA Extractor (Trizol) (Sangon Biotech, Shanghai, China) according to the manufacturer's instructions, then Qubit2.0 (Life Technologies, Thermo Fisher Scientific, USA) and agarose gel were used to detect the concentration and integrity of total RNA, respectively. After that, the RNA Mix with 3' Adapter was incubated 2 min at 70°C, and then T4 RNA Ligase 2(New England Biolabs, USA) was used to connect the 3' Adapter at 22°C for 1 h. Similarly, T4 RNA Ligase 1(New England Biolabs, USA) was used to connect the 5' Adapter at 20°C for 1 h. Then the connected product was used to reverse transcription with 5 × First Strand Buffer and PCR amplification (Table 5). The amplified cDNA library was used to detect the integrity with 12%-page gel, and was sequenced on the Illumina platform (Illumina Hiseq X Ten).

TABLE 5

Table 5. The information on the primers used in qRT-PCR validation.

Quality control, reads mapping, and sequence analysis

FastQC was used for quality checks on raw data of nine samples. Cutadapt and trimmomatic (69) were then used to remove the adaptors and low-quality reads, respectively. Blastn was used to align to Rfam (70) database and then the reads that mapped to the bovine rRNA, sRNA, snRNA, and snoRNA were discarded. Bowtie (71) was used to map to the bovine reference genome(ARS-UCD1.3), exon, and intron sequence, the reads that can be mapped to the bovine genome and the intron region were retained. Novel miRNAs prediction and known miRNAs quantification were separately run with the miRDeep2 module and quantifier module of the miRDeep2 software (72).

Differential expression analysis

The R package edgeR (73) was utilized to compare the expression level of miRNAs between H_vs_P and M_vs_P. Specifically, the lowly expressed miRNAs were filtered with the CPM normalization method, and the DE miRNAs among cattle, yak, and cattle-yak were screened with P-value < 0.05, |fold change| >2. We also carried on the Venn analysis to merge the common DE miRNAs between H_vs_P and M_vs_P for further analysis.

Target prediction and function annotation

With the filtered DE miRNAs of H_vs_P and M_vs_P, miranda was used to predict their target genes, and then ClusterProfiler (74) was used to analyze their function with GO and KEGG enrichment. In our previous study, we have already preliminarily explored the function of the novel lncRNA target genes and differentially expressed genes (DEGs) (38). To further understand the DEGs expression patterns and their function, we utilized the gene set enrichment analysis (GSEA) to demo some meaningful processes and pathways.

Net work analysis of DE miRNAs target genes

To understand the interactions among function processes and genes with pathways intuitively, igraph was used to draw the net plots and describe the correlations among pathways. Especially, the GO to GO, gene to pathways, and pathway to pathway plots were diagramed, respectively.

qPCR validation and statistical analysis

To verify the accuracy of sequencing results, miRNAs and several target genes were randomly selected for qPCR validation. The identified DE miRNAs miR-449a and the target genes FANCA, MYRFL, and SHF of miR-7 were selected for further qRT-PCR experiment with specific primers and conditions (Tables 5, 6). The 2^−ΔΔCT method was used to quantify the relative expression of the miRNAs and the target genes. Student's t-test with Graph Pad Prism 6 software was used to compare the differences among cattle, yak, and cattle-yak samples. ^*p-value < 0.05 was considered statistically significant and ^**p-value < 0.01 was considered as extremely statistically significant.

TABLE 6

Table 6. PCR reaction system and conditions.

Conclusions

In this study, we used high-throughput sequencing technology to explore miRNAs testicular tissue profiles of cattle, yaks, and cattle-yak. Finally, 47 up-regulated miRNAs and 34 down-regulated miRNAs were identified in cattle between the comparison of cattle and cattle-yak; 16 up-regulated and 21 down-regulated miRNAs were identified in yak samples between the comparison of yak and cattle-yak. Subsequently, we merged the same differentially expressed miRNAs between H vs. P and M vs. P. As a result, the known miRNA bta-miR-7 was commonly down-regulated in cattle-yak in the two comparison groups, revealing that miR-7 and its target genes MYRFL, FANCA, INSL3, and USP9X may play an extremely important role in reproductive processes. We also conducted GO, KEGG, and GSEA to annotate the function of differentially expressed genes and target candidates identified in nine samples, the results showed that the pathways and terms were associated with spermatogenesis, and reproduction processes. We also found that the membrane raft (GO:0045121) process, focal adhesion (bta04510) pathways, and ITGB3 gene were in the hub locations of the network. Our research is important to explore the molecular mechanism of cattle-yak male sterility and provide effective strategies to resolve it.

Data availability statement

The datasets analyzed for this study have been deposited in the Genome Sequence Archive (Genomics, Proteomics & Bioinformatics 2021) in National Genomics Data Center (Nucleic Acids Res 2022), China National Center for Bioinformation / Beijing Institute of Genomics, Chinese Academy of Sciences (GSA: CRA006824) that are publicly accessible at https://ngdc.cncb.ac.cn/gsa (75, 76).

Ethics statement

The animal study was reviewed and approved by Institutional Animal Care and Use Committee at the College of Animal Science and Technology, Sichuan Agricultural University (Permit Number: DKY2020050). Written informed consent was obtained from the owners for the participation of their animals in this study.

Author contributions

SZ, WS, and XJ performed the data analysis and drafted the manuscript. SZ, XJ, and YL collected the testicular samples involved in this study. YL and S-YC contributed to laboratory experiments. SL and S-YC revised the manuscript. XJ and JW designed the study and also revised the manuscript. All authors read and approved the final manuscript.

Funding

This research was funded by the National Key R&D Program of China (2021YFD1200403), the Sichuan innovation team of the national modern agricultural industry technology system (SCCXTD-2022-13), and the Key R&D Program of Sichuan Province (2021YFYZ0007).

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher's note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fvets.2022.974703/full#supplementary-material

References

1. Boon RA, Jaé N, Holdt L, Dimmeler S. Long noncoding rnas: from clinical genetics to therapeutic targets? J Am Coll Cardiol. (2016) 67:1214–26. doi: 10.1016/j.jacc.2015.12.051

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Cech Thomas R, Steitz Joan A. The noncoding RNA revolution—trashing old rules to forge new ones. Cell. (2014) 157:77–94. doi: 10.1016/j.cell.2014.03.008

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Lv Y, Huang S. Role of non-coding RNA in pancreatic cancer (Review). Oncol Lett. (2019) 18:3963–73. doi: 10.3892/ol.2019.10758

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Lv Y, Wang Z, Zhao K, Zhang G, Huang S, Zhao Y. Role of noncoding RNAs in cholangiocarcinoma (Review). Int J Oncol. (2020) 57:7–20. doi: 10.3892/ijo.2020.5047

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Pagani M, Rossetti G, Panzeri I, de Candia P, Bonnal RJ, Rossi RL, et al. Role of microRNAs and long-non-coding RNAs in CD4(+) T-cell differentiation. Immunol Rev. (2013) 253:82–96. doi: 10.1111/imr.12055

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Lu TX, Rothenberg ME. MicroRNA. J Allergy Clin Immunol. (2018) 141:1202–7. doi: 10.1016/j.jaci.2017.08.034

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Kim DY, Sung JH. Regulatory role of microRNAs in the proliferation and differentiation of adipose-derived stem cells. Histol Histopathol. (2017) 32:1–10. doi: 10.14670/HH-11-798

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Luo G, Hu S, Lai T, Wang J, Wang L, Lai S. MiR-9-5p promotes rabbit preadipocyte differentiation by suppressing leptin gene expression. Lipids Health Dis. (2020) 19:126. doi: 10.1186/s12944-020-01294-8

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Raza SHA, Kaster N, Khan R, Abdelnour SA, El-Hack MEA, Khafaga AF, et al. The role of microRNAs in muscle tissue development in beef cattle. Genes. (2020) 11:295. doi: 10.3390/genes11030295

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Guo L, Huang Q, Zhao J, Liu H, Lu W, Wang J. microRNA-10b promotes the apoptosis of bovine ovarian granulosa cells by targeting plasminogen activator inhibitor-1. Theriogenology. (2021) 176:206–16. doi: 10.1016/j.theriogenology.2021.09.035

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Yao Q, Chen Y, Zhou X. The roles of microRNAs in epigenetic regulation. Curr Opin Chem Biol. (2019) 51:11–7. doi: 10.1016/j.cbpa.2019.01.024

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Gu X, Sun W, Yi K, Yang L, Chi F, Luo Z, et al. Comparison of muscle lipidomes between cattle-yak, yak, and cattle using UPLC–MS/MS. J Food Compos Anal. (2021) 103:104113. doi: 10.1016/j.jfca.2021.104113

CrossRef Full Text | Google Scholar

13. Sun N, Chen J, Wang D, Lin S. Advance in food-derived phospholipids: Sources, molecular species and structure as well as their biological activities. Trends Food Sci Technol. (2018) 80:199–211. doi: 10.1016/j.tifs.2018.08.010

CrossRef Full Text | Google Scholar

14. Vahmani P, Ponnampalam EN, Kraft J, Mapiye C, Bermingham EN, Watkins PJ, et al. Bioactivity and health effects of ruminant meat lipids. Invited Rev Meat Sci. (2020) 165:108114. doi: 10.1016/j.meatsci.2020.108114

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Xu C, Wu S, Zhao W, Mipam T, Liu J, Liu W, et al. Differentially expressed microRNAs between cattleyak and yak testis. Sci Rep. (2018) 8:592. doi: 10.1038/s41598-017-18607-0

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Delph LF, Demuth JP. Haldane's Rule: Genetic Bases and Their Empirical Support. J Heredity. (2016) 107:383–91. doi: 10.1093/jhered/esw026

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Shah MA, Xu C, Wu S, Zhao W, Luo H, Yi C, et al. Isolation and characterization of spermatogenic cells from cattle, yak and cattleyak. Anim Reprod Sci. (2018) 193:182–90. doi: 10.1016/j.anireprosci.2018.04.067

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Cai X, Yu S, Mipam T, Yang F, Zhao W, Liu W, et al. Comparative analysis of testis transcriptomes associated with male infertility in cattleyak. Theriogenology. (2017) 88:28–42. doi: 10.1016/j.theriogenology.2016.09.047

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Phakdeedindan P, Wittayarat M, Tharasanit T, Techakumphu M, Shimazaki M, Sambuu R, et al. Aberrant levels of DNA methylation and H3K9 acetylation in the testicular cells of crossbred cattle-yak showing infertility. Reprod Domest Anim. (2021) 57:304–13. doi: 10.1111/rda.14061

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Sato Y, Kuriwaki R, Hagino S, Shimazaki M, Sambuu R, Hirata M, et al. Abnormal functions of Leydig cells in crossbred cattle-yak showing infertility. Reprod Domest Anim. (2020) 55:209–16. doi: 10.1111/rda.13609

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Kee K, Angeles VT, Flores M, Nguyen HN, Reijo Pera RA. Human DAZL, DAZ and BOULE genes modulate primordial germ-cell and haploid gamete formation. Nature. (2009) 462:222–5. doi: 10.1038/nature08562

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Yao W, Li Y, Li B, Luo H, Xu H, Pan Z, et al. Epigenetic regulation of bovine spermatogenic cell-specific gene boule. PLoS ONE. (2015) 10:e0128250-e. doi: 10.1371/journal.pone.0128250

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Davies B, Hatton E, Altemose N, Hussin JG, Pratto F, Zhang G, et al. Re-engineering the zinc fingers of PRDM9 reverses hybrid sterility in mice. Nature. (2016) 530:171–6. doi: 10.1038/nature16931

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Li B, Wu W, Luo H, Liu Z, Liu H, Li Q, et al. Molecular characterization and epigenetic regulation of Mei1 in cattle and cattle–yak. Gene. (2015) 573:50–6. doi: 10.1016/j.gene.2015.07.021

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Gu Y, Li Q, Pan Z, Li M, Luo H, Xie Z. Molecular cloning, gene expression and methylation status analysis of PIWIL1 in cattle-yaks and the parental generation. Anim Reprod Sci. (2013) 140:131–7. doi: 10.1016/j.anireprosci.2013.05.010

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Girard A, Sachidanandam R, Hannon GJ, Carmell MA, A. germline-specific class of small RNAs binds mammalian Piwi proteins. Nature. (2006) 442:199–202. doi: 10.1038/nature04917

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Megosh HB, Cox DN, Campbell C, Lin H. The role of PIWI and the miRNA machinery in Drosophila germline determination. Curr Biol. (2006) 16:1884–94. doi: 10.1016/j.cub.2006.08.051

PubMed Abstract | CrossRef Full Text | Google Scholar

28. Brennecke J, Malone CD, Aravin AA, Sachidanandam R, Stark A, Hannon GJ. An epigenetic role for maternally inherited piRNAs in transposon silencing. Science. (2008) 322:1387–92. doi: 10.1126/science.1165171

PubMed Abstract | CrossRef Full Text | Google Scholar

29. Deng W, Lin H. miwi, a murine homolog of piwi, encodes a cytoplasmic protein essential for spermatogenesis. Dev Cell. (2002) 2:819–30. doi: 10.1016/S1534-5807(02)00165-X

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Robert N, Yan C, Si-Jiu Y, Bo L, He H, Pengfei Z, et al. Expression of Rad51 and the histo-morphological evaluation of testis of the sterile male cattle-yak. Theriogenology. (2021) 172:239–54. doi: 10.1016/j.theriogenology.2021.06.018

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Abane R. Mezger V. Roles of heat shock factors in gametogenesis and development. FEBS J. (2010) 277:4150–72. doi: 10.1111/j.1742-4658.2010.07830.x

PubMed Abstract | CrossRef Full Text | Google Scholar

32. Dayalan Naidu S. Dinkova-Kostova AT. Regulation of the mammalian heat shock factor 1. FEBS J. (2017) 284:1606–27. doi: 10.1111/febs.13999

PubMed Abstract | CrossRef Full Text | Google Scholar

33. Liu P, Yu S, Cui Y, He J, Zhang Q, Sun J, et al. Regulation by Hsp27/P53 in testis development and sperm apoptosis of male cattle (cattle-yak and yak). J Cell Physiol. (2018) 234:650–60. doi: 10.1002/jcp.26822

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Zhao W, Mengal K, Yuan M, Quansah E, Li P, Wu S, et al. Comparative RNA-Seq analysis of differentially expressed genes in the epididymides of yak and cattleyak. Curr Genomics. (2019) 20:293–305. doi: 10.2174/1389202920666190809092819

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Zhang F, Hanif Q, Luo X, Jin X, Zhang J, He Z, et al. Muscle transcriptome analysis reveal candidate genes and pathways related to fat and lipid metabolism in Yunling cattle. Anim Biotechnol. 2021:1-8. doi: 10.1080/10495398.2021.2009846

PubMed Abstract | CrossRef Full Text | Google Scholar

36. Wu S, Mipam T, Xu C, Zhao W, Shah MA Yi C, et al. Testis transcriptome profiling identified genes involved in spermatogenic arrest of cattleyak. PLoS ONE. (2020) 15:e0229503. doi: 10.1371/journal.pone.0229503

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Xu C, Shah MA, Mipam T, Wu S, Yi C, Luo H, et al. Bovid microRNAs involved in the process of spermatogonia differentiation into spermatocytes. Int J Biol Sci. (2020) 16:239–50. doi: 10.7150/ijbs.38232

PubMed Abstract | CrossRef Full Text | Google Scholar

38. Zhao S, Chen T, Luo X, Chen S, Wang J, Lai S, et al. Identification of Novel lncRNA and Differentially Expressed Genes (DEGs) of Testicular Tissues among Cattle, Yak, and Cattle-Yak Associated with Male Infertility. Animals. (2021) 11:2420. doi: 10.3390/ani11082420

PubMed Abstract | CrossRef Full Text | Google Scholar

39. Lai EC. Micro RNAs are complementary to 3′ UTR sequence motifs that mediate negative post-transcriptional regulation. Nat Genet. (2002) 30:363–4. doi: 10.1038/ng865

PubMed Abstract | CrossRef Full Text | Google Scholar

40. Xie X, Lu J, Kulbokas EJ, Golub TR, Mootha V, Lindblad-Toh K, et al. Systematic discovery of regulatory motifs in human promoters and 3' UTRs by comparison of several mammals. Nature. (2005) 434:338–45. doi: 10.1038/nature03441

PubMed Abstract | CrossRef Full Text | Google Scholar

41. Lin Q, Ma L, Liu Z, Yang Z, Wang J, Liu J, et al. Targeting microRNAs: a new action mechanism of natural compounds. Oncotarget. (2016) 8:15961–70. doi: 10.18632/oncotarget.14392

PubMed Abstract | CrossRef Full Text | Google Scholar

42. Cai Y, Yu X, Hu S, Yu J, A. brief review on the mechanisms of miRNA regulation. Genom Proteom Bioinform. (2009) 7:147–54. doi: 10.1016/S1672-0229(08)60044-3

PubMed Abstract | CrossRef Full Text | Google Scholar

43. Pantos K, Grigoriadis S, Tomara P, Louka I, Maziotis E, Pantou A, et al. Investigating the role of the microRNA-34/449 family in male infertility: a critical analysis and review of the literature. Front Endocrinol (Lausanne). (2021) 12:709943. doi: 10.3389/fendo.2021.709943

PubMed Abstract | CrossRef Full Text | Google Scholar

44. Capra E, Turri F, Lazzari B, Cremonesi P, Gliozzi TM, Fojadelli I, et al. Small RNA sequencing of cryopreserved semen from single bull revealed altered miRNAs and piRNAs expression between High- and Low-motile sperm populations. BMC Genomics. (2017) 18:14. doi: 10.1186/s12864-016-3394-7

PubMed Abstract | CrossRef Full Text | Google Scholar

45. He Y, Duan L, Wu H, Chen S, Lu T, Li T, et al. Integrated transcriptome analysis reveals the impact of photodynamic therapy on cerebrovascular endothelial cells. Front Oncol. (2021) 11:731414. doi: 10.3389/fonc.2021.731414

PubMed Abstract | CrossRef Full Text | Google Scholar

46. Suh YJ, Choe JY, Park HJ. Malignancy in pheochromocytoma or paraganglioma: integrative analysis of 176 cases in TCGA. Endocr Pathol. (2017) 28:159–64. doi: 10.1007/s12022-017-9479-2

PubMed Abstract | CrossRef Full Text | Google Scholar

47. Wen DY, Lin P, Pang YY, Chen G, He Y, Dang YW, et al. Expression of the long intergenic non-protein coding RNA 665 (LINC00665) gene and the cell cycle in hepatocellular carcinoma using the cancer genome atlas, the gene expression omnibus, and quantitative real-time polymerase chain reaction. Med Sci Monit. (2018) 24:2786–808. doi: 10.12659/MSM.907389

PubMed Abstract | CrossRef Full Text | Google Scholar

48. Tran N, Broun A, Ge K. Lysine demethylase KDM6A in differentiation, development, and cancer. Mol Cell Biol. (2020) 40:e00341–20. doi: 10.1128/MCB.00341-20

PubMed Abstract | CrossRef Full Text | Google Scholar

49. Wang C, Lee JE, Cho YW, Xiao Y, Jin Q, Liu C, et al. UTX regulates mesoderm differentiation of embryonic stem cells independent of H3K27 demethylase activity. Proc Natl Acad Sci U S A. (2012) 109:15324–9. doi: 10.1073/pnas.1204166109

PubMed Abstract | CrossRef Full Text | Google Scholar

50. Welstead GG, Creyghton MP, Bilodeau S, Cheng AW, Markoulaki S, Young RA, et al. X-linked H3K27me3 demethylase Utx is required for embryonic development in a sex-specific manner. Proc Natl Acad Sci U S A. (2012) 109:13004–9. doi: 10.1073/pnas.1210787109

PubMed Abstract | CrossRef Full Text | Google Scholar

51. Lu C, Zhang Y, Qin Y, Xu Q, Zhou R, Cui Y, et al. Human X chromosome exome sequencing identifies BCORL1as contributor to spermatogenesis. J Med Genet. (2021) 58:56–65. doi: 10.1136/jmedgenet-2019-106598

PubMed Abstract | CrossRef Full Text | Google Scholar

52. Yang X, Zhang X, Jiao J, Zhang F, Pan Y, Wang Q, et al. Rare variants in FANCA induce premature ovarian insufficiency. Hum Genet. (2019) 138:1227–36. doi: 10.1007/s00439-019-02059-9

PubMed Abstract | CrossRef Full Text | Google Scholar

53. Larder R, Chang L, Clinton M, Brown P. Gonadotropin-releasing hormone regulates expression of the DNA damage repair gene, fanconi anemia A, in pituitary gonadotroph cells1. Biol Reprod. (2004) 71:828–36. doi: 10.1095/biolreprod.104.030569

PubMed Abstract | CrossRef Full Text | Google Scholar

54. Inoue Y, Shingyoji C. The roles of noncatalytic ATP binding and ADP binding in the regulation of dynein motile activity in flagella. Cell Motil Cytoskeleton. (2007) 64:690–704. doi: 10.1002/cm.20216

PubMed Abstract | CrossRef Full Text | Google Scholar

55. Bhattacharjee R, Goswami S, Dey S, Gangoda M, Brothag C, Eisa A, et al. Isoform-specific requirement for GSK3α in sperm for male fertility†. Biol Reprod. (2018) 99:384–94. doi: 10.1093/biolre/ioy020

PubMed Abstract | CrossRef Full Text | Google Scholar

56. Suzuki M, Suzuki H, Sugimoto Y, Sugiyama Y. ABCG2 transports sulfated conjugates of steroids and xenobiotics. J Biol Chem. (2003) 278:22644–9. doi: 10.1074/jbc.M212399200

PubMed Abstract | CrossRef Full Text | Google Scholar

57. Bart J, Hollema H, Groen HJ, de Vries EG, Hendrikse NH, Sleijfer DT, et al. The distribution of drug-efflux pumps, P-gp, BCRP, MRP1 and MRP2, in the normal blood-testis barrier and in primary testicular tumours. Eur J cancer (Oxford, England: 1990). (2004) 40:2064–70. doi: 10.1016/j.ejca.2004.05.010

PubMed Abstract | CrossRef Full Text | Google Scholar

58. Caballero J, Frenette G, D'Amours O, Dufour M, Oko R, Sullivan R. ATP-binding cassette transporter G2 activity in the bovine spermatozoa is modulated along the epididymal duct and at ejaculation1. Biol Reprod. (2012) 86:181. doi: 10.1095/biolreprod.111.097477

PubMed Abstract | CrossRef Full Text | Google Scholar

59. Chang L, Karin M. Mammalian MAP kinase signalling cascades. Nature. (2001) 410:37–40. doi: 10.1038/35065000

PubMed Abstract | CrossRef Full Text | Google Scholar

60. Ni F-D, Hao S-L, Yang W-X. Multiple signaling pathways in Sertoli cells: recent findings in spermatogenesis. Cell Death Dis. (2019) 10:541. doi: 10.1038/s41419-019-1782-z

PubMed Abstract | CrossRef Full Text | Google Scholar

61. Qu F, Cui Y, Zeng J, Zhang M, Qiu S, Huang X, et al. Acupuncture induces adenosine in fibroblasts through energy metabolism and promotes proliferation by activating MAPK signaling pathway via adenosine(3) receptor. J Cell Physiol. (2020) 235:2441–51. doi: 10.1002/jcp.29148

PubMed Abstract | CrossRef Full Text | Google Scholar

62. Hanks SK, Calalb MB, Harper MC, Patel SK. Focal adhesion protein-tyrosine kinase phosphorylated in response to cell attachment to fibronectin. Proc Natl Acad Sci U S A. (1992) 89:8487–91. doi: 10.1073/pnas.89.18.8487

PubMed Abstract | CrossRef Full Text | Google Scholar

63. Parsons JT. Focal adhesion kinase: the first ten years. J Cell Sci. (2003) 116:1409–16. doi: 10.1242/jcs.00373

PubMed Abstract | CrossRef Full Text | Google Scholar

64. Staub C, Johnson L. Review: Spermatogenesis in the bull. Animal. (2018) 12:s27–35. doi: 10.1017/S1751731118000435

PubMed Abstract | CrossRef Full Text | Google Scholar

65. Li Y-C, Wang G-W, Xu S-R, Zhang X-N, Yang Q-E. The expression of histone methyltransferases and distribution of selected histone methylations in testes of yak and cattle-yak hybrid. Theriogenology. (2020) 144:164–73. doi: 10.1016/j.theriogenology.2020.01.001

PubMed Abstract | CrossRef Full Text | Google Scholar

66. Wang S, Pan Z, Zhang Q, Xie Z, Liu H. Li Q. Differential mRNA expression and promoter methylation status of SYCP3 gene in testes of yaks and cattle-yaks. Reprod Domest Anim. (2012) 47:455–62. doi: 10.1111/j.1439-0531.2011.01902.x

PubMed Abstract | CrossRef Full Text | Google Scholar

67. Liu Z, Li Q, Pan Z, Qu X, Zhang C, Xie Z. Comparative analysis on mRNA expression level and methylation status of DAZL gene between cattle-yaks and their parents. Anim Reprod Sci. (2011) 126:258–64. doi: 10.1016/j.anireprosci.2011.05.013

PubMed Abstract | CrossRef Full Text | Google Scholar

68. Kishi K, Uchida A, Takase HM, Suzuki H, Kurohmaru M, Tsunekawa N, et al. Spermatogonial deubiquitinase USP9X is essential for proper spermatogenesis in mice. Reprod (Cambridge, England). (2017) 154:135–43. doi: 10.1530/REP-17-0184

PubMed Abstract | CrossRef Full Text | Google Scholar

69. Bolger AM, Lohse M, Usadel B. Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics (Oxford, England). (2014) 30:2114–20. doi: 10.1093/bioinformatics/btu170

PubMed Abstract | CrossRef Full Text | Google Scholar

70. Griffiths-Jones S, Bateman A, Marshall M, Khanna A, Eddy SR. Rfam: an RNA family database. Nucleic Acids Res. (2003) 31:439–41. doi: 10.1093/nar/gkg006

PubMed Abstract | CrossRef Full Text | Google Scholar

71. Langmead B, Trapnell C, Pop M, Salzberg SL. Ultrafast and memory-efficient alignment of short DNA sequences to the human genome. Genome Biol. (2009) 10:R25. doi: 10.1186/gb-2009-10-3-r25

PubMed Abstract | CrossRef Full Text | Google Scholar

72. Friedländer MR, Mackowiak SD Li N, Chen W, Rajewsky N. miRDeep2 accurately identifies known and hundreds of novel microRNA genes in seven animal clades. Nucleic Acids Res. (2011) 40:37–52. doi: 10.1093/nar/gkr688

PubMed Abstract | CrossRef Full Text | Google Scholar

73. Robinson MD, McCarthy DJ, Smyth GK. edgeR: A Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics. (2009) 26:139–40. doi: 10.1093/bioinformatics/btp616

PubMed Abstract | CrossRef Full Text | Google Scholar

74. Yu G, Wang L-G, Han Y, He Q-Y. clusterProfiler: an R package for comparing biological themes among gene clusters. OMICS. (2012) 16:284–7. doi: 10.1089/omi.2011.0118

PubMed Abstract | CrossRef Full Text | Google Scholar

75. Chen T, Chen X, Zhang S, Zhu J, Tang B, Wang A, et al. The genome sequence archive family: toward explosive data growth and diverse data types. Genom Proteom Bioinform. (2021) 19:578–83. doi: 10.1101/2021.06.29.449849

PubMed Abstract | CrossRef Full Text | Google Scholar

76. CNCB-NGDC Members and Partners. Database resources of the national genomics data center, china national center for bioinformation in 2022. Nucleic Acids Res. (2022) 50:D27–D38. doi: 10.1093/nar/gkab951

PubMed Abstract | CrossRef Full Text | Google Scholar