Promoter hypermethylation of PIWI/piRNA pathway genes associated with diminished pachytene piRNA production in bovine hybrid male sterility

ABSTRACT

Hybrid male sterility (HMS) is a postzygotic reproductive isolation mechanism that enforces speciation. A bovine example of HMS is the yattle (also called dzo), an interspecies hybrid of taurine cattle (Bos taurus) and yak (Bos grunniens). The molecular mechanisms underlying HMS of yattle are not well understood. Epigenetic modifications of DNA methylation and P-element induced wimpy testis (PIWI)-interacting RNA (piRNAs) are important regulators in spermatogenesis. In this study, we investigated DNA methylation patterns and piRNA expression in adult testes in hybrid infertile yattle bulls and fertile cattle and yak bulls using whole genome bisulphite-seq and small RNA-seq. Promoter hypermethylation in yattle were associated with DNA methylation involved in gamete generation, piRNA metabolic processes, spermatogenesis, and spermatid development (P < 2.6 × 10⁻⁵). Male infertility in yattle was associated with the promoter hypermethylation-associated silencing of PIWI/piRNA pathway genes including PIWIL1, DDX4, PLD6, MAEL, FKBP6, TDRD1 and TDRD5. The downstream effects of silencing these genes were diminished production of 29- to 31- nucleotide pachytene piRNAs in yattle testes. Hypermethylation events at transposable element loci (LINEs, SINEs, and LTRs) were found in yattle. LINE-derived prepachytene piRNAs increased and SINE-derived prepachytene piRNAs were reduced in yattle testes. Our data suggests that DNA methylation affects the PIWI/piRNA pathway and is involved in gene expression and pachytene piRNA production during spermatogenesis in bovine HMS. DNA hypermethylation and disruption of piRNA production contributed to unsuccessful germ cell development that may drive bovine HMS.

KEYWORDS:

• Hybrid male sterility
• DNA methylation
• piRNA biogenesis
• transposon

Introduction

Hybrid male sterility (HMS) is a postzygotic reproductive isolation mechanism that enforces speciation by restricting gene flow between species or subspecies. It is thought to result from the acquisition of genetic incompatibilities [1,2]. HMS is a universal phenomenon in many eukaryotic interspecies hybrids [3]. Domestic yaks (Bos grunniens) and domestic taurine cattle (Bos taurus) are closely related species that diverged 4.9 million years ago. The synteny of their genomes is as high as 94% [4]. An example of HMS in bovines is the yattle (also called dzo), the interspecies hybrid of taurine cattle and yak. F1 hybrid yattle and F2 backcross males are sterile while females are fertile. Seminiferous tubules of yattle contain few primary spermatocytes at the leptotene to pachytene stage and no spermatocytes beyond this stage [5–7]. Thus, the male infertility of yattle may be due to arrest of meiosis at about the pachytene stage [7]. Testes transcriptomics and proteomics show that downregulated genes primarily associated with meiosis and spermiogenesis are poorly expressed in sterile yattle [8,9]. However, the regulatory network that controls germline transcription in yattle is not understood.

Interspecific hybridization causes genomic stress that can lead to genetic instability in offspring. Accumulating evidence indicates that hybrids show extensive modifications of the genome and transcriptome, including novel gene expression pattern changes, chromosomal rearrangements, transposable element mobilization, miRNA and P-element induced wimpy testis (PIWI)-interacting RNA (piRNAs) deficiency, and DNA methylation changes (reviewed in [10]). Epigenetic modifications of DNA methylation, histone modification, and piRNAs are important regulators in spermatogenesis [11–13]. DNA methylation involves the transfer of a methyl group from S-adenosyl-l-methionine to a cytosine residue and generally occurs in CpG dinucleotide sites [14]. DNA methylation in promoters is usually restricted to genes in long-term, stable repressed states [15] and promoter DNA hypermethylation generally suppresses transcription [16]. The methylation state of promoters of genes involved in spermatogenesis has been studied. Hypermethylation of the promoter regions of methylenetetrahydrofolate reductase (MTHFR) and voltage dependent anion channel 2 (VDAC2) down-regulate MTHFR and VDAC2 expression, leading to nonobstructive azoospermia and potentially idiopathic infertility in humans [17,18]. Studies show that some meiosis genes have significantly reduced expression in testes of yattle and have hypermethylation in their promoters and CpG islands (CGI), including DAZL (deleted in azoospermia like), SYCP3 (synaptonemal complex protein 3), PIWIL1 (piwi-like RNA-mediated gene silencing 1), DDX4 (DEAD-box helicase 4), MEI1 (meiotic double-stranded break formation protein 1), and FKBP6 (FK506 binding protein 6) using bisulphite sequencing and RT-qPCR [19–23]. Genome-wide DNA methylation in testes of yattle and cattle has not been investigated.

Two studies demonstrated that the promoter hypermethylation-associated silencing of PIWI/piRNA pathway genes PIWIL1, PIWIL2, and TDRD1 is associated with male infertility [24] and human testicular tumorigenesis [25]. Promoter hypermethylation of PIWI/piRNA pathway genes diminished piRNA production and transposon loci such as LINE-1 are hypomethylated in human testicular tumorigenesis and spermatogenic disorders [24,25]. PIWI proteins are a subfamily of the Argonaute protein family that is expressed mainly in germline cells and forms specific RNA-induced silencing complexes (RISCs) with piRNAs that are termed piRISCs. PiRNAs are a class of small RNAs of 25–32 nucleotides (nt) that are highly abundant in testes [26–28] and predominantly contribute to spermatogenesis in a stage-specific manner [29,30]. During mouse spermatogenesis, two distinct populations of piRNAs associate with the PIWI proteins PIWIL1/MIWI, PIWIL2/MILI, and PIWIL4/MIWI2 at two different developmental stages: prepachytene piRNAs and pachytene piRNAs [31]. Prepachytene piRNAs are mainly derived from transposable elements. They are associated with MILI and MIWI2 in the gonocyte stage and found in leptotene spermatocyte after birth. They have 26- to 27-nt MILI-bound modal lengths [26] and a 28-nt MIWI2-bound modal lengths [32]. Pachytene piRNAs originate from piRNA clusters in various regions of the genome and bind to both MILI and MIWI in pachytene spermatocytes to the round spermatid stage. They have 30-nt MIWI-bound modal lengths and are the majority of piRNAs in adult mouse testes [26,27]. In mammalian testes, piRISCs in germ granules regulate gene silencing and maintain genome stability by directing CpG DNA methylation at transposon loci [13,32–34]. Prepachytene piRNAs direct DNA and histone methylation at transposon sequences during embryogenesis [32,35,36]. Pachytene piRNAs post-transcriptionally repress transposons late in spermatogenesis [33,37] and ensure completion of meiosis and successful spermiogenesis [34,38,39]. These results suggest that DNA methylation linked to the PIWI/piRNA pathway is crucial for regulating germline transcription and may contribute to male infertility disorders.

Whether DNA methylation and piRNA production are associated with bovine HMS is currently unclear. To address this matter, we investigated DNA methylation patterns and piRNA expression in hybrid infertile yattle bulls and fertile cattle and yak bulls using whole genome bisulphite-seq (WGBS) and small RNA-seq (sRNA-seq). Promoter hypermethylation of PIWI/piRNA pathway genes was associated with transcriptional silencing and diminished pachytene piRNA expression in yattle testes. Our data indicated that epigenetic disruption of the PIWI/piRNA pathway may drive bovine HMS.

Results

Global view of WGBS data

Yattle were interspecies hybrid of taurine cattle (B. Taurus, ♂) and yaks (B. grunniens, ♀) used in this study, and males were sterile. To investigate DNA methylation patterns in hybrid infertile bulls and fertile parent bulls, testicular tissues from one infertile yattle and one fertile cattle were used for WGBS. A total of 109.39 G raw bases were generated for the yattle sample and 120.91 G for the cattle sample. After data filtering, 358,346,324 clean reads were generated for yattle and 395,585,807 for cattle. The bisulphite sequencing conversion rates for the groups were as high as 99.9%. Mapping reads covered 77.62% (278,148,416 reads) of the cattle genome (Bos_taurus_UMD_3.1.1) for yattle and 81.65% (322,995,811) for cattle. Sequencing depths were 25.40× for yattle and 29.43× for cattle. Approximately 3% of all genomic cytosine sites were methylated in yattle and cattle. Methylation in yattle and cattle occurred in three sequence contexts: CpG, CHG (where H is A, C, or T), and CHH. Only CpG sites were used for further analysis.

Identification of DMRs between sterility and fertility testes

To highlight aberrant methylated CpG sites in yattle testes and allow downstream processing and analyses, differential methylation regions (DMRs) were determined. We identified 47,252 hypermethylated and 108,231 hypomethylated DMRs for CpG in yattle testes compared to cattle. We used the top 10,000 significant DMRs of CpG for further analysis.

The length distribution of the top 10,000 DMRs was major from 250 bp to 900 bp (Figure 1(a)). The 10,000 DMRs were annotated in 24,857 genomic regions including CGI, CGI shores, promoters, transcription start sites (TSSs), 5ʹ untranslated regions (utr5), exons, introns, utr3, transcription termination sites (TESs), and repeats (Figure 1(b)). We mapped 70.22% of hypermethylated and 55.83% of hypomethylated DMRs to genomic regions with no gene annotation (intergenic regions). Low percentages of DMRs, specifically 1.96% (446) of hypermethylated and 3.64% (78) of hypomethylated, DMRs mapped to promoter regions of known genes (Figure 1(c)). Among the top 10,000 DMRs, 9142 were hypermethylated and 858 were hypomethylated in yattle. The methylation level of DMRs was higher in yattle than in cattle as seen with heatmaps and boxplots (Figure 1(d–e)). The chromosomal distribution of hypermethylated and hypomethylated DMRs is shown using Circos in figure 1(f). DMRs were not evenly distributed in the genome. The ends of chromosomes had preferentially hypermethylated or hypomethylated DMRs in yattle. Chromosomes 1, 7, 8, 9, and X exhibited fewer hypomethylated DMRs. The top 10,000 DMRs of CpG were annotated in 9255 repeats. Mainly short interspersed elements (SINEs), long interspersed elements (LINEs), long terminal repeats (LTRs) retrotransposons were more hypermethylated in yattle than cattle testes (Figure 1(g)). LINE1, LINE (RET) and SINE (Core) mRNA were extensively expressed in all tested tissues (Figure 1(h), left). RT-qPCR showed that LINE1, LINE (RET) and SINE (Core) mRNA expression were reduced in yattle testes (Figure 1(h), right).

Figure 1. Characteristics of differential methylation regions (DMRs) between adult yattle and cattle testes.

(a) Length distribution of DMRs. (b–c) Genomic region distribution of DMR between adult yattle (CY) and cattle (C) testes. Hypermethylation, methylation levels higher in yattle than cattle. (d–e) Heatmap and boxplot of DMR methylation levels in yattle and cattle. (f) Chromosomal distribution of DMRs. (g) Methylation level of repeats in yattle and cattle testes. (h) Left: RNA samples from cattle, yak and yattle tissues showing expression patterns. Te: testes; He: heart; Sp: spleen; Lu: lung; Li: liver; Mu: muscle; Ma: marker. Right: mRNA levels of LINE (L1), LINE (RTE-BovB), and SINE (Core-RTE) in testes of cattle, yaks, and yattle.

Overlap between methylated promoters and CpG islands

Although the effect of DNA methylation in intergenic regions is less clear, aberrant methylation in promoter regions is linked to altered chromosomal state and thus to gene transcription silencing [16]. In this study, 3206 DMRs (2956 hypermethylated and 250 hypomethylated) were at CGIs and 524 DMRs (446 hypermethylated and 78 hypomethylated) were at promoters. Among them, 264 hypermethylated and 28 hypomethylated gene promoters overlapped with CGIs, and were 59.2% (264/446) of the hypermethylated and 35.9% (28/78) of the hypomethylated promoters. These results indicated that more genes exhibited differential methylation in promoters that were not associated with CGIs. Therefore, aberrant promoter methylation was not restricted to CGIs.

Promoter hypermethylation of piRNA metabolism and spermatogenesis was associated with male sterility in yattle

When regions of a DMR and a specific gene functional element overlapped, the corresponding gene was selected as a DMR-related, differentially methylated gene (DMG). A total of 7710 CpG DMRs were annotated to 2848 DMGs, and more DMGs were observed in gene bodies than in promoters (Figure 2(a)). Among the DMGs, 2362 were more prone to hypermethylation in yattle than cattle (Figure 2(a)). The 419 DMGs displayed hypermethylation in their promoters (Figure 2(a)).

Figure 2. Promoter hypermethylation in yattle involved in PIWI/piRNA pathway and spermatogenesis.

(a) Venn diagram of hypermethylation and hypomethylation of promoter and gene bodies. (b) Enriched GO terms and rich factor of 419 promoter-hypermethylated genes. (c) MRNA levels of piRNA pathway genes in adult yattle, cattle, and yak testes. FC: Fold-change for yattle to cattle or yaks. Bold, hypermethylation in promoter and CGI. (d) MRNA levels of piRNA pathway genes in adult yattle, cattle, and yak testes using RT-qPCR. (e) MRNA levels for piRNA pathway genes in purified mouse spermatogenic cell populations. 2 C: most contained somatic cells; LZ: leptotene and zygotene spermatocytes; PS: pachytene spermatocytes; RS: round spermatids. Data from da Cruz et al., 2016. (f) RNA samples from cattle, yak and yattle tissues showing expression patterns. Te: testes; He: heart; Sp: spleen; Lu: lung; Li: liver; Mu: muscle; Ma: marker.

To examine the potential function of 419 promoter-hypermethylated genes, we performed functional analysis. Significantly gene ontology (GO)-enriched terms were in spermatogenesis, meiotic nuclear division, meiosis I, meiotic chromosome segregation, DNA methylation involved in gamete generation, piRNA metabolic process, spermatid development, synapsis, and homologous chromosome segregation (Figure 2(b), P < 2.3 × 10⁻⁵, FDR P < 0.04). Although only one yattle and one cattle sample were used for WGBS, and not from yak. Our data identified hypermethylated promoters and CGIs for DAZL, SYCP3, PIWIL1, DDX4 and FKBP6 in yattle testes compared to cattle. Hypermethylated promoters and CGIs for these genes were identified in yattle compared to yak and cattle using bisulphite sequencing in previous studies [19–23]. This result gave us confidence in our WGBS data. These findings suggested that hypermethylation of genes at promoters in testes of yattle was involved in spermatogenesis, meiosis I, gamete generation and piRNA metabolic process.

Promoter hypermethylation of PIWI/piRNA pathway genes in yattle were associated with transcriptional silencing

Basing on GOTERM_DIRECT analysis, promoter-hypermethylated genes were significantly enriched for DNA methylation involved in gamete generation (GO: 0043046), piRNA metabolic process (GO: 0034587), chromatoid body (GO: 0033391), and pi-body (GO: 0071546). Fold-enrichment factors were higher than 30 (Table S1). PIWIL1, PIWIL2, PIWIL4, DDX4, MAEL, PLD6, TDRD1, TDRD5, TDRD9, TDRD12, RNF17, FKBP6, MOV10L1, ASZ1, TDRKH, HENMT1, EXD1, GPAT2, PNLDC1, BTBD18, and HSP90 were involved the piRNA pathway (reviewed in [13,40]). Among these genes, well-known piRNA metabolic genes DDX4, MAEL, PLD6, PIWIL2, TDRD1, TDRD5, RNF17, FKBP6, and MOV10L1 were hypermethylated in promoters and CGIs (Table S1). PIWIL1 and TDRD9, were hypermethylated in exons or near CGIs. The mRNA levels of PIWIL1, DDX4, MAEL, PLD6, TDRD1, TDRD5, TDRD9, RNF17, FKBP6 and MOV10L1 genes were significantly lower in yattle than cattle and yaks (Figure 2(c), fold change > 2). The mRNA levels were from unpublished RNA-Seq data (Figure 2(c), n = 3 for cattle, yaks, and yattle, including same samples as the WGBS) and RT-qPCR data (Figure 2(d), n = 8 for cattle, yaks, and yattle, including same samples as the WGBS and sRNA-seq).

In mice, DDX4, RNF17, FKBP6, MOV10L1, TDRD9, MAEL genes were turned on in leptotene and zygotene spermatocytes (LZ; meiotic prophase I early stages). Peak expression of PIWIL1, PLD6, TDRD1, TDRD5 was in pachytene spermatocytes (PS; medium meiotic prophase) (Figure 2(e), data from [41]). We were unable to obtain highly purified spermatogenic cell populations to investigate stage-specific expression of PIWI/piRNA genes in bovine. RT-PCR agarose-gel results showed that PIWIL1, PIWIL4, DDX4, RNF17, FKBP6, MOV10L1, MAEL, PLD6, TDRD1, TDRD5 and TDRD9 were uniquely expressed in adult testes of cattle, yak and yattle (figure 2(f)). Overall, promoter hypermethylation of PIWI/piRNA pathway genes was associated with transcriptional silencing in yattle.

Small RNA sequencing data

The disruption of genes associated with the piRNA processing machinery is directly related to spermatogenesis failure due to maturation arrest, resulting in male sterility in mouse models (reviewed in [13,40]). To determine promoter hypermethylation-associated transcriptional silencing of PIWI/piRNA pathway genes and effects on piRNA production, we analysed small RNA-seq (sRNA-seq) libraries from two testicular tissues samples from yattle, cattle, and yaks by deep sequencing. After removing low-quality reads and adaptors, we obtained an average of 19.67 million clean reads from cattle libraries, 17.68 million from yattle libraries, and 20.36 million from yak libraries (Table S2). For sRNA analysis, we used an average 19.32 million reads for cattle, 17.29 million for yattle, and 19.99 million for yaks, ranging from 18 to 35 nt (Table S2). Clean reads from 18 to 35 nt were mapped to the cattle genome (Bos_taurus_UMD_3.1.1). Mapping read coverage was above 91% for cattle and yattle and about 88% for yaks (Table S2).

For each sRNA-seq library, more than 94% reads were 21 to 32 nt. The frequency of percentages of approximately 22- and 30-nt modal lengths were different for yattle, cattle and yaks (Figure 3(a)). In cattle and yak fertile bulls, the average frequencies of percentages were 11% of total sRNA reads for 21–24 nt, 19% for 25–28, and 65% for 29–32 modal lengths. However, for hybrid infertile yattle, the average frequency percentage for 21- to 24-nt modal length was up to 51% and the 29- to 32-nt modal length decreased to 23% of total sRNA reads (Figure 3(a)).

Figure 3. Reduced piRNA levels in adult yattle testes.

(a) Length distribution and percentage of small RNA reads from adult cattle, yak and yattle testes. (b) piRNA base bias at each position. (c) PiRNAs exhibited a 5′-end U bias at position 1 in cattle, yaks and yattle adult testes. (d) Length distribution and percentage of piRNAs reads from adult cattle, yak and yattle testes. (e) Genomic annotation of piRNAs from adult cattle, yak and yattle testes. Sequence reads (25–32 nt) from piRNA reads aligned to cattle genomic sequence sets in the order: piRNA clusters, exon, repeats, intronic sequences and other. Percentage of mapped reads is shown. (f) Examples of monodirectional clusters (clusters 5129 and 2680) and a bidirectional cluster (cluster 1821). DNA methylation levels for yattle and cattle of these regions is shown.

Pachytene piRNA production was reduced in adult yattle testes

In this study, piRNA sequences were identified using a k-mer scheme [42]. We obtained an average of 0.65 million distinct piRNA species for cattle, 0.56 million for yattle, and 0.68 million for yaks. An average of 5.29 million piRNA sequence reads were obtained for cattle, 1.80 million for yattle, and 5.36 million for yaks (Table S2). Primary piRNAs show an antisense and 1 U bias, secondary piRNAs loaded onto AGO3 show 10-nt complementarity at their 5ʹ ends with PIWI-bound piRNAs and possess a 10A bias that are referred to as the ping-pong signature [27,43]. PiRNA has a strong bias for U at the first position in this study (Figure 3(b)). More than 93.0% of piRNA sequences began with U at the 5ʹ-end in cattle and yak testes (Figure 3(c)). The percent of 1 U bases decreased about 2.25% and 4.93% from two biological replicates in yattle. Meanwhile, the percent of 10A bases (Figure 3(c)) and mRNA expression of AGO3 (Figure 2(d)) were no significantly different between three groups. These results indicated that piRNA 5′ formation and the ping-pong cycle were normal in cattle, yaks and yattle.

In mice, modal lengths are 26–27 nt for MILI (PIWIL2)-bound piRNAs (prepachytene piRNAs) [26,32,44], 28 nt for MIWI2 (PIWIL4)-bound piRNAs [32], and 30 nt for MIWI (PIWIL1)-bound piRNAs (pachytene piRNAs) [31,45]. Because no datasets were available for small RNAs from different stages of bovine spermatogenesis, we used length of reads to computationally separate prepachytene from pachytene piRNAs. Pachytene piRNAs are the vast majority of piRNAs in adult mouse testes and accumulate rapidly at the pachytene stage of meiosis [31]. More than 73% of piRNAs with 29- to 31-nt modal lengths were present in adult cattle and yak testes. The relative percentages of 29- to 31-nt pachytene piRNAs were diminished in yattle, while 26- to 27-nt modal length piRNAs were increased in yattle testes (Figure 3(d)). Together, these data established that pachytene piRNAs production was diminished in adult yattle testes.

Characteristics of piRNA clusters in bovine adult testes

Previous work showed that piRNAs are primarily generated from hundreds of unique genomic loci (piRNA clusters) through a primary processing pathway [27,31,46]. PiRNA clusters contain a ten to thousands of piRNAs, varying in size from 1 to 100 kb [47,48]. Using a coverage cut-off of two and a minimum piRNA cluster length cut-off of 200, we identified 9222 cluster regions from 200 bp to 181.5 kb. Putative cluster-piRNAs mapped uniquely in more than 68% of total piRNAs in cattle and 59% in yaks (Figure 3(e)). Also, 3% piRNAs mapped to repeats of SINEs, LINEs, LTRs retrotransposons. Only 48% of putative piRNAs were located in piRNA clusters in yattle. Up to 5.48% of piRNAs were in LINEs in yattle, with 1.22% in cattle and 1.27% in yaks (Figure 3(e)).

To increase confidence in piRNA clusters, the cut-off for piRNA clusters was set to transcripts per kilobase of exon model per million mapped reads (TPM) ≥50 for two biological replicates of cattle and yak samples. We identified 288 high-confidence piRNA clusters. The piRNA clusters showed a nonuniform distribution in chromosomes and strong strand asymmetry. Most clusters encoding piRNAs were mainly on only one strand; there were named monodirectional clusters. In some clusters, one strand changed abruptly to another strand; there were named bidirectional clusters [27,30]. For example, cluster 2680 (chr8:17538987–17667396) and cluster 5129 (chr15:76578895–76698234) were found mainly on one strand, while cluster 1821 (chr5:106856947–106973147) was a bidirectional cluster. We found no overlap between DMRs and piRNA cluster regions (figure 3(f)).

Top pachytene piRNA clusters were diminished in adult yattle testes

We used read lengths to computationally separate prepachytene from pachytene piRNA clusters. Clusters were considered prepachytene if the fraction of 26-nt reads was larger than the fraction of 30-nt reads, and vice versa for pachytene clusters [49]. Using only this parameter, we found among 288 clusters, 32 with the length of presumed prepachytene piRNA clusters (26–28 nt), and 256 (88.9%) with the length of presumed pachytene piRNA clusters (29–31 nt) mapped to intergenic regions. As a pachytene piRNA cluster example, the top 50 high-expression piRNA from cluster 5129 were diminished in adult yattle testes. A few piRNAs of cluster 5129 were mapped to repeat sequences (Figure 4(a)). Cluster 5129 was in an intergenic region and repeat sequences showed scattered distribution in this genomic region (Figure 4(b)). These results were consistent with the finding that pachytene piRNAs are poor in transposon sequences [26,27]. In contrast, the top 50 high-expression piRNAs from prepachytene piRNA cluster 3647 were increased in adult yattle testes (Figure 4(c)). The genomic region of cluster 3647 was highly enriched in LINE sequences (Figure 4(d)). By read numbers, pachytene clusters had about 10 times more piRNA reads than prepachytene reads.

Figure 4. Cluster-derived piRNAs were diminished in yattle testes.

(a and c) Comparison of top 50 high-expression piRNAs in clusters 5129 and 3647 in adult cattle and yattle testes. (b and d) Characteristics of piRNA clusters 5129 and 3647 in piRBase. (e) Number of piRNA reads mapped to top 60 high-expression piRNA clusters from small RNA libraries of adult cattle, yak and yattle testes. (f) Comparison of fold-change in top 60 high-expression piRNA-producing clusters in cattle vs. yattle and yaks vs. yattle.

An average of 66.6% piRNAs reads were derived from the 288 high-confidence piRNA clusters in adult cattle and an average of 57.5% from yak testes. Yattle testes had a specific reduced percentage of 40.6% of these cluster-piRNAs. In cattle and yak testes, the top 60 high-confidence piRNA clusters were the majority (>97%) of the 288 cluster-derived piRNAs and were diminished in yattle (Figure 4(e), Table S2). In yattle testes, piRNAs from these top 25 piRNA-producing clusters were uniformly reduced by an average of six-fold (Figure 4(f)).

LINE-derived prepachytene piRNAs were increased in yattle testes

PIWI/piRNA complexes recruit methylation factors to silence transposons in male germline stem cells. Thus, the piRNA-related machinery is directly involved in regulation of repetitive elements [32,35,36]. We annotated transposable elements (TEs) in the cattle reference genome. From sRNA libraries, an average 2.32 million reads (12.92% of total mapped reads) from cattle, 1.02 million (6.55%) from yattle, and 2.15 million (12.24%) from yaks were aligned to repeat reads (satellite, SINE, LTR, LINE, and DNA). Mapped SINE, LTR, and LINE reads to total sRNA reads were significantly lower in yattle than in cattle and yaks (Figure 5(a), Table S2). The percent of piRNA originating from LINEs was increased in yattle (Figure 3(d)). We classified repeat sequences of noncluster regions by length into 26–28 nt and 29–32 nt. The percent of 26-to 28-nt LINE-derived piRNA reads to total piRNA reads was increased in yattle (Figure 5(b)). Among the top 90 high expression LINE-derived piRNAs, expression of 26- to 28-nt piRNAs was higher in yattle than cattle (Figure 5(c)). Expression of 29- to 32-nt LINE-derived and SINE-derived piRNAs was significantly reduced in yattle (Figure 5(d)). We found no overlap between DMRs and LINE-derived and SINE-derived piRNA regions (data not shown).

Figure 5. Transposable element-driving piRNAs in adult cattle, yak and yattle testes.

(a–b) Percent of annotation repeats (LINEs, SINEs and LTRs) in small RNA libraries (a) or total piRNA reads (b) of adult cattle, yak and yattle testes. (c–d) Number of top 90 high-expression LINE-driving (c) or top 65 high-expression SINE-driving (d) piRNA reads in small RNA libraries of adult cattle, yak and yattle testes.

Discussion

Understanding the process of speciation is a central problem in biology and animal genetics. Interspecies HMS is studied for speciation because sterility is a postzygotic isolation barrier to free gene flow between organisms; thus, effectively isolating distinct species [1,2]. A related question in HMS is why genes that cause fertility in a pure-species genetic background fail to produce fertility in a hybrid genetic background [1]. To address this matter, gene expression arrays have been used to establish associations between misregulation of gene expression and sterility in hybrids between closely related species. In Drosophila, sterile hybrids are more variable in gene expression than fertile males and most misregulated genes are underexpressed [50,51]. Interspecies sterile Drosophila hybrids have mostly postmeiotic problems during spermatogenesis. This connection between sterility phenotype and tissue-specific misregulation of gene expression suggests that downregulation of postmeiotic spermatogenesis genes in sterile hybrids is a causative factor of HMS [1,51].

We hypothesized that gene repression by hypermethylation of male meiosis specific genes, particularly the PIWI/piRNA pathway genes for pachytene piRNA production PIWIL1, DDX4, PLD6, MAEL, FKBP6, RNF17, MOV10L1, TDRD1, TDRD5, and TDRD9, probably drive bovine HMS. WGBS and GOTERM_DIRECT biological process analysis showed promoter hypermethylated genes in sterile yattle were significantly enriched in ‘DNA methylation involved in gamete generation (GO: 0043046)’ and ‘piRNA metabolic process (GO: 0034587)’ (Figure 2 (b), Table S1). Therefore, abnormal DNA methylation impacted spermatogenesis failure in yattle. Gene-promoter DNA hypermethylation generally suppresses transcription [16]. In several mutant mouse models, PIWIL1 [27], DDX4 [52,53], MAEL [54], PLD6 [55], TDRD1 [56,57], TDRD5 [58], TDRD9 [52], RNF17 [59], FKBP6 [60], MOV10L1 [61,62], BTBD18 [63] and TDRKH [64] are important in pachytene piRNA production and spermatogenesis. Pachytene piRNAs accumulate rapidly at the pachytene stage of meiosis [31]. Pachytene piRNAs may promote spermatogenesis by regulating mRNAs and long noncoding RNAs in mouse testes and ensure completion of meiosis and successful spermiogenesis [34,38,39]. Indeed, promoter and CGI hypermethylation associated with silencing PIWI/piRNA pathway genes PIWIL1, DDX4, PLD6, MAEL, FKBP6, RNF17, MOV10L1, TDRD1, TDRD5, and TDRD9 (Figure 2(c)), resulted in diminished production of pachytene piRNAs in yattle testes (Figures 3(d) and 4(e). This molecular connection with histologic observations of seminiferous tubules indicated that male infertility of yattle may be due to meiosis arrest at about the pachytene stage [7]. In contrast, PIWIL2, PIWIL4, EXD1, TDRD12 are involved in prepachytene piRNA production [65–68]. PNLDC1 is essential for piRNA 3'-end trimming and transposon silencing during spermatogenesis in mice [69]. HENMT1 is required for piRNA 3'-methylation and stability [70]. Corresponding to the functions of these genes, mRNA levels of BTBD18, TDRKH, PIWIL2, PIWIL4, EXD1, TDRD12, PNLDC1, HENMT1 and piRNA length were not significantly different between yattle, cattle and yaks (Figures 2(c) and 3(a)). These results indicated that promoter hypermethylation of specific pachytene piRNA-biogenesis genes had a direct impact on pachytene piRNA production. In addition, all loss-of-function mutations of PIWI/piRNA pathway genes reported so far in mice lead to male sterility, while females are fertile. These results indicate the importance of this pathway in male germline development in mammals (reviewed in [13]). PIWIL1, FKBP6, DDX4, MAEL, TDRD5 and TDRD9 were uniquely expressed in adult testes of cattle, yak and yattle (Figure 2(f)). This result may explain why males were sterile while females were fertile in yattle. Finally, PIWI/piRNA pathway genes predominantly contribute to spermatogenesis in a stage-specific manner [29,30]. Transcriptome analysis of highly purified mouse spermatogenic cell populations showed that DDX4, RNF17, FKBP6, MOV10L1, TDRD9, MAEL genes were turned on in leptotene and zygotene spermatocytes (meiotic prophase I early stages), while the expression peak of PIWIL1, PLD6, TDRD1, TDRD5 was in pachytene spermatocytes (medium meiotic prophase) [41] (Figure 2(e)). Although no bovine datasets of gene expression from different stages of spermatogenesis were available, higher expression is seen for PIWIL1 in leptotene and zygotene spermatocytes and PIWIL2 in spermatogonial stem cells and leptotene spermatocytes during human spermatogenesis [71]. Thus, PIWI/piRNA pathway genes may conservatively contribute to bovine spermatogenesis in a stage-specific manner. Overall, promoter hypermethylation of male-specific and pachytene piRNA biogenesis-specific genes have a direct impact on pachytene piRNA production, which may have contributed to pachytene arrest in yattle. Studies report that promoter hypermethylation-associated silencing of PIWI/piRNA pathway genes PIWIL1, PIWIL2, and TDRD1 and piRNA production are associated with male infertility [24] and human testicular tumorigenesis [25]. Our study showed that DNA hypermethylation of PIWI/piRNA pathway genes was associated with diminished pachytene piRNA production in bovine HMS.

The Dobzhansky-Muller incompatibility model explains HMS as the result of incompatible epistatic interactions between divergent genomes. This model of genetic incompatibilities has long provided a useful theoretical framework for speciation genetics. However, only one mammalian gene, PRDM9, is responsible for sterility phenotypes in male hybrids of certain mouse subspecies [72,73]. In addition, interspecific hybridization is a genomic stress condition that leads to activation of TEs in animals. TE insertions can disrupt genes or regulatory sequences and can cause chromosomal rearrangements, threatening host genome integrity [74]. To mitigate these deleterious effects, mechanisms of TE control are especially important in the germline. In mice, CpG DNA methylation and PIWI/piRNA complexes safeguard the male germline genome from TEs [75]. From gonocytes to leptotene spermatocytes, MILI and MIWI2 are associated with transposon sequence-rich TE piRNAs, and guide specific de novo DNA methylation to suppress TEs and maintains germline genome integrity [32,36,48]. In contrast, MIWI and MILI are highly expressed in pachytene spermatocytes and associated with pachytene piRNA species [76]; both are essential for adult spermatogenesis in mice [77,78]. MIWI and MILI carry out LINE1 transposon silencing functions in pachytene spermatocytes, likely through slicer activity in a piRNA-independent manner [33,37]. In our study, we found that 26- to 27-nt modal length piRNAs (Figure 2(d)) and 26–28 nt LINE-driver piRNAs expression (Figures 4(c) and 5(c)) were increased in yattle. Meanwhile, the LINEs and SINEs were hypermethylated (Figure 1(g)) and the LINE and SINE mRNA expression (Figure 1(g)) were reduced in yattle. We found no overlap between top 10,000 DMRs and 288 high-confidence piRNA cluster regions (Figure 3(f)). PIWIL1 expression and pachytene piRNAs production were diminished in yattle testes. These results suggested higher TE piRNAs expression may affect TE DNA methylation and lead to suppress TEs expression in yattle. However, the regulation mechanisms between TE piRNAs and DNA methylation could not be confirmed in the whole testis data.

Conclusions

In this work, we used WGBS and sRNA-seq to investigate DNA methylation and piRNA production differences in the testes of yattle, cattle and yaks. Promoter hypermethylation genes were significantly enriched in DNA methylation involved in gamete generation, piRNA metabolic process, chromatoid body, and pi-body. We confirmed that promoter hypermethylation-associated silencing of PIWI/piRNA pathway genes had the potential to reduce pachytene piRNA production in yattle. Our data indicated that epigenetic disruption of the entire PIWI/piRNA pathway was involved with bovine HMS. Gene repression by hypermethylation of male meiosis-specific genes and PIWI/piRNA pathway genes is probably a driver of bovine HMS.

Materials and methods

Sample collection

Testes tissue samples from taurine cattle (B. taurus, n = 8), yak (B. grunniens, n = 8), and yattle (B. taurus ♂ × B. grunniens ♀, n = 8) were used to obtain genomic DNA and total RNA. Detail information is in our previous study [7]. Yaks and yattle bulls were from the Hongyuan Bovine Breeding Farm in Sichuan, China. Hybrid taurine cattle were from a beef farm in Chongqing, China. The genetic relationship between individuals was unclear. Animals were slaughtered after sexual maturity (4.5 years for yaks and yattle and 2.5 years for taurine cattle). Fat and fascia surrounding testes were resected and two cross-sections of testicular tissue from the middle of the testes were obtained by fine-scale dissection. Crosscut slice was snap frozen in liquid nitrogen for transport and stored at −80°C. All experimental procedures were approved by the Institutional Animal Care and Use Committee of Southwest University (Chongqing, China).

Whole-genome bisulphite sequencing

We analysed one cattle (fertile) and one yattle (infertile) testicular DNA sample for genome-wide DNA methylation by WGBS. DNA extraction, library preparation and sequencing was performed at Novogene (Beijing, China). DNA was purified with tissue DNA extraction kits (Tiangen, cat: DP304). DNA concentrations were determined by a Qubit 2.0 fluorometer (Life Technologies, USA). Libraries were prepared from 5.2 mg genomic DNA spiked with 26 ng lambda DNA fragmented by sonication to 200–300bp with a Covaris S220. Libraries were prepared using Illumina TruSeq DNA methylation kits (Illumina, USA). DNA was terminal-repaired and A-ligated, followed by ligation of methylated adapters according to the Illumina TruSeq protocol. Libraries were bisulphite converted twice using EZ DNA methylation gold kits (Zymo Research, cat: D5005). Library concentration was quantified by a Qubit 2.0 flurometer (Life Technologies) and quantitative PCR. Insert sizes were determined on an Agilent Bioanalyzer 2100 system. Libraries were analysed in separate lanes on a HiSeq 4000 (Illumina, USA) in paired-end mode using 150-bp sequence length. Raw sequencing data were deposited to NCBI in SRA.

Processing WGBS data

Sequence data processing was according to procedures of Novogene (Beijing). Following sequencing, we used FastQC (fastqc_v0.11.5, Babraham Institute, UK) to perform basic statistics on raw read quality. We obtained high-quality reads after trimming low-quality bases from the 3′ end and adapter sequences from the 5′ end using Trimmomatic 0.36 software and parameters SLIDINGWINDOW: 4:15; LEADING:3, TRAILING:3; ILLUMINACLIP: adapter.fa: 2:30:10; and MINLEN:36. Resulting clean reads were aligned to the reference B. taurus cattle genome (UMD_3.1.1,) using Bismark alignment software (version 0.16.3) [79], using the stringent option -X 70 with slightly relaxed parameters for mapping (score-min L, 0, −0.2). Initial data analysis was by Bowtie2 [80]. Methylation level (ML) for each C site was defined as: ML = mC/(mC+umC). This study focused on CpG sites.

DMRs were identified based on a P-value cut-off of 0.05 using DDS software [81]. A sliding-window approach (selecting every 200 bp as a window with 100-bp as step lengths) was used to identify DMRs in the genome. Regions with length greater than 50 bp, corrected P-values less than 1e-05, and read coverage greater than 3 were regarded as DMRs. According to DMR distribution through the genome, DMRs were annotated as promoter, exon, intron, CGI, CGI shore, repeat, TSS, TES. We defined genes related to DMRs as having gene body regions (from TSS to TES) or promoter regions (2kb upstream from the TSS) that overlapped DMRs. When a region where a DMR and a specific gene function element overlapped, the corresponding gene was selected as the DMR-related DMG.

GO enrichment and KEGG pathway analyses were conducted using the DAVID 6.8 (https://david.ncifcrf.gov/) for DMGs to investigate biological processes and functions [82]. GO terms and KEGG were corrected using the Benjamini method. Corrected P-values less than 0.05 were considered significantly enriched for DMGs. Graphs for mRNA expression from mRNA-seq were produced with GraphPad Prism 7.0.

Small RNA sequencing

We analysed testicular RNA sample from fertile cattle (n = 2) and yaks (n = 2), and infertile yattle (n = 2) for miRNA and piRNA using small RNA-seq. Total RNA extraction, library preparation and sequencing were performed at Novogene (Beijing, China). Total RNA from testes were extracted using TRIzol reagent (Invitrogen, cat: 15596026). The quantity and purity of total RNA were analysed by Nanodrop one (ThermoFisher Scientific, USA) and agarose gel electrophoresis. RNA concentration was measured using Qubit RNA Assay Kits in a Qubit 2.0 fluorometer (Life Technologies). RNA integrity was assessed using RNA Nano 6000 Assay Kits for the Agilent Bioanalyzer 2100 system (Agilent Technologies, USA) with RNA integrity number >8.0. Total RNA at 3 μg per sample was used as input for small RNA libraries. Sequencing libraries were generated using NEBNext Multiplex Small RNA Library Prep Sets for Illumina (NEB, cat: E7300 S). Libraries were purified on 8% polyacrylamide gels (100 V, 80 min) and quantified using an Agilent Bioanalyzer 2100 system using DNA High Sensitivity Chips. Library preparations were sequenced on an Illumina Hiseq 2500 platform, and 125-bp single-end reads generated. Raw data were deposited to NCBI in SRA.

Processing of sRNA sequencing data

Sequence data processing was according to procedures at Novogene (Beijing). Clean reads were obtained by removing reads containing poly-N; with 5′-adapter contaminants; without 3′-adapters or insert tags; containing poly A, T, G or C; or having low quality in raw data. Simultaneously, Q20, Q30, and GC-content of raw data were calculated. A length range of 18–35 nt for clean reads was selected for all downstream analyses.

Small RNA tags were mapped to the reference sequence B. taurus cattle genome (UMD_3.1.1) using Bowtie (v0.12.9, -v 0 v 0 -k) [83] without mismatches to analyse expression and distribution. Small RNAs were mapped to single annotations by performing annotations in the order: miRNAs > rRNAs > tRNAs > snRNAs > snoRNAs > repeats > genes > novel miRNAs. Mapped small RNA tags were used to search for known miRNAs using miRBase21.0 as a reference.

Processing of piRNA and piRNA clusters

PiRNA and piRNA clusters were analysed based on ‘repeat’ and ‘other’ sRNA data of 25–32 nt. PiRNA sequences were identified using a k-mer scheme [42]. To obtain the first base preference map of piRNAs of different lengths, AUCG statistics were performed on identified piRNA sequences. Identified piRNA sequences were aligned with the reference genome to obtain the distribution of piRNAs on the chromosome using Bowtie (-v 0 -k 1). The coverage of piRNA on the reference sequence was obtained using samtools depth. PiRNA-generating genes were obtained from the sam file based on alignments of piRNAs mapped to exons and with coverage ≥ 5. We annotated TEs in the cattle reference genome using RepeatMasker (version open-4.0.1).

Parameters for identifying piRNA clusters were coverage ≥ 2 piRNAs, minimum piRNA cluster length 200 nt, and threshold distance of interval length 10,000 nt. Neighbouring genes within piRNA clusters were detected based on the range of each piRNA cluster, considering sequences 2000 bp upstream and downstream. Differential expression of piRNA clusters among yattle, cattle and yaks testes were measured with DESeq2 [84], setting padj < 0.05 as the cut-off. The piRNA name was identified in piRBase using the ‘Seq To Name’ tool [85]. Graphs for first-base bias and piRNA cluster expression differences were produced with GraphPad Prism 7.0.

mRNA expression arrays

mRNA expression profiles of target genes in different tissues was performed from testes, muscle, liver, spleen, heart, and lung tissues using reverse transcription PCR (RT-PCR) as described in previous studies [7,86]. The mRNA expression of the target genes was determined by reverse transcription quantitative PCR (RT-qPCR) in cattle, yaks, and yattle testes as described in previous studies [7,86]. Primers for all these genes were designed using Primer3 software and synthetized by Genewiz (Suzhou, China, Table S3). Relative mRNA expression was analysed with the system’s relative quantification software based on the 2^{delta-delta-Ct} method (Bio-Rad). Graphs for mRNA expression differences were produced with GraphPad Prism 7.0.

Author summary

Domestic yaks (Bos grunniens) and domestic taurine cattle (Bos taurus) are closely related species. The interspecies F1 hybrid yattle and F2 backcross males are sterile. Previous studies showed that the male infertility of yattle may be due to arrest of meiosis at about the pachytene stage and meiosis gene expression is dysregulated. Epigenetic modifications of DNA methylation, histone modification, and piRNAs are important regulators in spermatogenesis. In this study, we investigated DNA methylation patterns and piRNA expression in adult testes in hybrid infertile yattle bulls and fertile cattle and yak bulls using whole genome bisulphite-seq and small RNA-seq. We confirmed that promoter DNA hypermethylation-associated silencing of PIWI/piRNA pathway genes had the potential to reduce pachytene piRNA production during spermatogenesis in yattle. Our data indicated that epigenetic disruption of the entire PIWI/piRNA pathway was involved with yattle male sterility. DNA hypermethylation of male meiosis-specific genes and disruption of piRNA production contributed to unsuccessful germ cell development that may drive bovine hybrid male sterility.

Authors’ contributions

GZ and FZ conceived and supervised the project. GZ, LW, HC, JG, YW, and JZ performed most experiments. JG, JZ, ZL, and WH provided support or technical assistance. GZ, LW, YW, and ZL contributed to bioinformatics data. GZ wrote the manuscript with major inputs from LW, YW, JZ, and FZ. All authors approved the final manuscript.

Ethics declaration

All experimental procedures were approved by the Institutional Animal Care and Use Committee of Southwest University (Chongqing, China).

Disclosure statement

No potential conflict of interest was reported by the authors.

Supplemental material

Supplemental data for this article can be accessed here.

Additional information

Funding

This work was financially supported by the National Natural Science Foundation of China [grant no. 31802046] and the Chongqing Research Program of Basic Research and Frontier Technology [cstc2017jcyjBX0013 and cstc2018jcyjAX0524] and Fundamental Research Funds for the Central Universities [XDJK2019RC001]. The funders had no role in study design, decision to publish, or preparation of the manuscript.