Introduction

Spermatogenesis is highly conserved among all vertebrates. Although it generally consists of three phases (mitotic, meiotic, and haploid), many characteristics appear to be species-specific, e.g., the duration of each of the three phases, the seminiferous epithelial organization, and the shape and length of spermatozoa, likely reflecting the adaptive changes during evolution (1–3). Several cellular events are unique to the male germ cells, e.g., postnatal formation of the adult male germline stem cells (i.e., spermatogonial stem cells), pubertal onset of meiosis, and haploid male germ cell differentiation (i.e., spermiogenesis) (4). Unique cellular processes are often accompanied by a more complex yet unique transcriptome, which may explain why the testis expresses more genes than any other organs in the body, with the possible exception of the brain (5). Regulation of gene expression during spermatogenesis occurs at both transcriptional and post-transcriptional levels (6). As a post-transcriptional regulator, miRNAs are abundantly expressed in the testis and are required for spermatogenesis (7–11). miRNAs typically function at post-transcriptional levels by binding the complementary sequences in the untranslated regions (UTRs) of mRNAs – particularly in the 3’UTRs through the “seed sequence” (2^nd-7^th nucleotides)(12). Numerous miRNAs are subject to rapid evolution, probably in response to the accelerated rate of divergence of UTRs compared to the exonic sequences (13). Divergence of genomic sequences can be mediated by transposable elements (TEs), which are known as building blocks of the genome and mostly map to UTRs and intronic regions of protein-encoding genes (14). Each miRNA can bind numerous target mRNAs, and one mRNA can be targeted by multiple different miRNAs. This “one-to-multiple” relationship between miRNAs and mRNAs amplifies their potential to coordinate gene expression in the cell (12). Moreover, miRNA genes often exist in clusters, which are transcribed as a unit followed by nuclear and cytoplasmic cleavage events to generate individual miRNAs (12).

Multiple clusters of miRNA genes containing the same seed sequences are categorized as a miRNA family, and miRNAs within the same family likely evolved from a common ancestor sequence (15). Of great interest, many of the testis-enriched miRNA clusters map to the X chromosome (16). Sex-linked genes are generally subject to the male germline-specific phenomenon called meiotic sex chromosome inactivation (MSCI), which silences transcription during most, if not all, of meiosis (16). Indeed, prior to 2009, there were no confirmed reports of any sex-linked genes escaping the repressive effects of MSCI. Surprisingly, however, we found that many X-linked miRNA genes do escape MSCI, suggesting that they may contribute to particularly important functions during spermatogenesis (16). This notion has since been tested by the generation of knockouts (KOs) of individual miRNA genes or individual clusters of miRNA genes normally expressed during spermatogenesis. However, these KOs resulted in minimal, if any, phenotypic effects and did not appear to impede normal spermatogenesis or male fertility (15). This left the field facing multiple unanswered questions, including 1) why do so many X-linked miRNAs express uniquely or preferentially during spermatogenesis and escape MSCI, 2) what’s their origin, and 3) how and why did they evolve rapidly?

To address these questions and better understand the functional role played by these X-linked miRNAs, we investigated the evolutionary history of this unique miRNA family and also generated KOs of individual, paired, triple, quadruple, or quintuple sets of miRNA clusters within this family and tested the effects on male fertility, initially by standard monandrous mating assays. Consistent with previous efforts to inactivate miRNA genes in C. elegans and mice (15, 17–19), KOs of either individual members or individual clusters of the miR-506 family induced no discernable phenotypes and did not impact male fertility. This may reflect the level of functional redundancy inherent within the members and clusters of this miRNA family. It was only when four or more clusters of the miR-506 family were ablated that relevant phenotype became detectable, which was manifested in the form of reduced litter size despite normal sperm counts, motility, and morphology.

Interestingly, the most common male fertility testing for lab rodents is based on a monandrous mating scheme, i.e., one fertility-proven female mated with one male. However, there are additional aspects of male fertility in many mammalian species, particularly those that are normally litter-bearing. In the wild, litter-bearing females often mate with multiple different males such that a single litter may include pups sired by more than one male (20, 21). This polyandrous mating introduces the potential for additional aspects of male reproductive fitness to accrue, one of which involves sperm competition (20, 21). Sperm competition can occur when sperm from more than one male are present in the female reproductive tract simultaneously, such that they then compete to fertilize each oocyte (22). Sperm competition has been now recognized as a significant evolutionary force directly impacting male reproductive success (22). Using experiments that mimic polyandrous mating, we found that the quinKO male mice indeed displayed compromised sperm competition. Hence, the X-linked miR-506 family miRNAs appear to function to finetune spermatogenesis to enhance sperm competition and, consequently, male reproductive fitness.

Results

X-linked miR-506 family miRNAs flanked by two highly conserved protein-coding genes Slitrk2 and Fmr1 rapidly evolved across species

X-linked genes are generally more divergent between species than autosomal ones, a phenomenon known as the “faster-X effect” (23). However, despite a high degree of conservation of two protein-coding genes, Slitrk2 and Fmr1, on the X chromosome across species, the miRNA genes located between these two loci are divergent among clades across the eutherian mammals (15, 24). Through tracing the evolution of this genomic region, we found that Slitrk2 and Fmr1 mapped to chromosome 4 (syntenic with mammalian X chromosome) in zebrafish and birds, but to the X chromosome in most mammals, with divergence and multiplication of numerous miRNA genes that belong to the miR-506 family in between (Fig. 1A and Table S1). By mapping these miRNAs of various species using the UCSC genome browser (25), we found that all members of the miR-506 family are located in a region flanked by Slitrk2 and Fmr1 (Fig. 1A). Consistent with previous reports (15, 24), Slitrk2 and Fmr1 are usually on the positive strand, whereas the miR-506 family miRNAs are in the reverse orientation (Fig. 1A). Based on the location of these miRNAs, we named the miRNAs proximal to Slitrk2 (miR-892c ∼ miR-891a in humans) and Fmr1 (miR-513c ∼ miR-514a3 in humans) SmiRs (Slitrk2-proximal miRNAs) and FmiRs (Fmr1-proximal miRNAs), respectively.

Figure 1.

To evaluate the sequence conservation of these miRNAs across species, we adopted the Multiz Alignment and Conservation pipeline, which utilizes PhastCons and PhyloP algorithms (25), to search miRNA datasets from 100 different species using the human genome as a reference (Fig. 1B and Fig. S1). The mean values of phyloP and phastCons of the Fmr1 and Slitrk2 coding sequences (CDS) were ∼4.9 and ∼0.97, respectively (Fig. 1C and D), indicating that these regions are highly conserved. In contrast, the mean values of phyloP and phastCons of SmiRs were ∼1.0 and ∼0.002, respectively, and those of FmiRs are ∼0.03 and ∼0.11, respectively (Fig. 1 C and D). The phyloP and phastCons values of the CDS, FmiRs, and SmiRs are significantly different from each other (adjusted p-value < 0.05, Kruskal-Wallis test) (Fig. 1 C and D), indicating that FmiRs and SmiRs are highly divergent, and SmiRs are more divergent than FmiRs.

We then assessed the genomic sequence similarity among various species using D-GENIES-based dot plot analyses (26). Although the Slitrk2-Fmr1 genomic regions were highly variable among different species, the sequences within some clades shared a high degree of similarities, e.g., primates (rhesus monkeys, chimpanzees, and humans), cetartiodactyla (sheep and cows), rodentia (e.g., mice and rats) and carnivora (e.g., dogs and cats) (Fig. S2A). Although the miR-506 family miRNAs were highly divergent, some orthologs displayed a higher degree of sequence conservation (Fig. S2 B and C), e.g., the SmiRs within primates were similar (Fig. S2B); miR-891a and miR-891b in rhesus monkeys were similar to miR-891b and miR-891a in humans and chimpanzees, respectively, and miR-892b in chimpanzees was homologous to miR-892c in humans and rhesus monkeys in terms of sequence and location (Figs. S2B). The FmiRs, including miR-201 (assigned as Mir-506-P1 (paralogue 1) in MirGeneDB (27)), miR-547 (Mir-506-P2), and miR-509 (Mir-506-P7) in mice and rats, are orthologues of miR-506, miR-507 and miR-509 in humans, respectively (Fig. S2C, and Table S1). Of interest, although the miR-506 family miRNAs are highly divergent, the seed sequences of some miRNAs, such as Mir-506-P6 and Mir-506-P7, remain conserved (Fig. S2C), and these miRNAs represent the dominant mature miRNAs (28). It is noteworthy that the majority of the substitutions among the miR-506 family are U↣C and A↣G (Fig. S2 B and C). Furthermore, we analyzed the conservation of the miR-506 family in modern humans using data from the 1000 Genomes Project (1kGP) (Fig. 1 E and F) (29). We compared SmiRs, FmiRs, and all miRNAs with pachytene piRNAs, which are known to be highly divergent in modern humans but barely exert any biological functions (30). Of interest, the derived allele frequency (DAF) and mean nucleotide diversity (MND) of FmiRs and all miRNAs were significantly smaller than that of the pachytene piRNAs, whereas SmiRs were significantly smaller than all miRNAs (Fig. 1 E and F) (adjusted p-value < 0.05, Kruskal-Wallis test), suggesting that the miR-506 family miRNAs are more conserved than pachytene piRNAs in modern humans. Taken together, these data indicate that the X-linked miR-506 family, although rapidly evolving as a whole, contains both divergent and conserved miRNAs, suggesting both conserved and novel functions across species.

X-linked miR-506 family miRNAs are derived from MER91C DNA transposons

To visualize the family history of these miRNAs, we built a phylogram for the miR-506 family (Fig. S3). The phylogram suggests that these miRNAs shared a common ancestor and that the FmiRs emerged earlier than the SmiRs, which is also supported by the fact that some FmiRs exist in green sea turtles (Figs. S1 and S3). These data suggest that the miR-506 family miRNAs arose much earlier than previously thought (25). The two subfamilies, FmiRs and SmiRs, despite their common ancestors, may have evolved at different paces and thus, might be functionally divergent.

Studies have shown that transposable elements (TEs) drive evolution through transpositions (31). CRISPR-Cas9/Cas12a genome editing can induce irreversible small indels at the cutting sites (also called “scars”)(32), which have been used for lineage tracing (33). Inspired by this strategy, we attempted to trace the evolution of the miR-506 family miRNAs by searching the transposon database for the transpositional “scars” (partial TE sequences) after transposition. To search the potential TE sources of the miR-506 family miRNAs, we downloaded all transposons in the human, horse, dog, and guinea pig genomes and aligned them to their corresponding miR-506 family miRNAs using BLAST (Basic Local Alignment Search Tool) (34). The nonautonomous MER91C DNA transposon (∼100-150 million years) (35, 36) was the only transposon that aligned to FmiRs of the miR-506 family (> 94% identical matches) in all four species (Table S2).

Given that the FmiRs (e.g., miRs-506∼509) emerged much earlier than the SmiRs (fig. S1 and fig. S3) and that miR-513 (belonging to FmiRs) and SmiRs (including miR-891a and miR-891b) share a common ancestor (Fig. S3), we reasoned that the X-linked miR-506 family might be derived from the MER91C DNA transposon. To test this hypothesis, we first aligned the X-linked miR-506 family miRNAs from several species to a human MER91C DNA transposon.

Indeed, numerous FmiRs of almost all species analyzed aligned to the MER91C DNA transposon despite few mismatches (Fig. 2A). The phylogenetic tree further confirmed that the MER91C is the sister group of the miR-506 family miRNAs (Fig. 2B and Fig. S4), suggesting that the MER91C DNA transposon is the likely source of the older miR-506 family miRNAs. Further supporting this notion, the MER91C DNA transposons could form hairpin structures, which is a prerequisite for miRNA biogenesis (Fig. 2C and Fig. S5A). Moreover, analyses of the testis small RNA datasets from humans, marmosets, dogs, and horses revealed the peaks corresponding to these miRNAs (Fig. 2D and Fig. S5A). Finally, by overexpressing several MER91C DNA regions randomly selected from humans, dogs, and horses in HEK293T cells, we found that these DNA regions were indeed capable of producing miRNAs (Fig. 2 E and F, and Fig. S5 B and C). Co-expression of MER91C DNA regions and AGO2 significantly increased the abundance of human MER91C miRNAs, as compared with overexpression of MER91C DNA regions alone (Fig. 2 E and F), suggesting that these miRNAs could be loaded onto and protected by AGO2. Taken together, these results indicate that the X-linked miR-506 miRNAs were originally derived from the MER91C DNA transposon.

Figure 2.

X-linked miR-506 family miRNAs are predominantly expressed in spermatogenic cells and sperm

Several previous studies have shown that the X-linked miR-506 family miRNAs are predominantly expressed in the testis of multiple species (15, 16, 24, 37–39). By analyzing the publicly available small RNA sequencing (sRNA-seq) datasets from multiple species, including humans, rhesus monkeys, mice, rats, rabbits, dogs, and cows (27, 40, 41), we further confirmed that miR-506 family miRNAs were indeed highly abundant in the testis, but barely expressed in other organs (Fig. S6, and Table S3). To further determine whether these miRNAs are expressed in male germ cells in rodent testes, we conducted sRNA-seq using pachytene spermatocytes, round spermatids, and sperm purified from adult mice (Fig. S7A and Table S3). Consistent with previous data (15, 16, 24), these miRNAs were abundantly expressed in spermatogenic cells in murine testes (Fig. 3A). ∼80% of these miRNAs were significantly upregulated (FDR < 0.05) when pachytene spermatocytes developed into round spermatids, and ∼83.3% were significantly upregulated (FDR < 0.05) when developed into cauda sperm (Fig. 3A). By analyzing the publicly available sRNA-seq datasets from humans (42), marmosets (37) and horses (43), we determined the expression patterns of the miR-506 family miRNAs in the testes of these species (Fig. 3 B, C and D, Table S3). The significantly increasing abundance of the SmiRs from immature to mature testes in horses (43) supports the elevated expression in haploid male germ cells (round, elongating/elongated spermatids, and sperm) compared to meiotic male germ cells (spermatocytes) (Fig. 3D). More interestingly, the SmiRs and FmiRs appear to be differentially expressed in the testes of various species, e.g., relative levels of the SmiRs were greater than those of the FmiRs in mice (Fig. 3 A and Table S3). Both the SmiRs and FmiRs were highly expressed in horses (Fig. 3D and Table S3). Levels of the FmiRs in marmoset (37) and human testes (42) were much greater than in those of the SmiRs (Fig. 3 B and C and Table S3). Overall, the miR-506 family miRNAs are abundant in the testis and predominantly expressed in haploid male germ cells, i.e., spermatids and spermatozoa.

Figure 3.

Ablation of X-linked miR-506 family miRNAs compromises male fertility due to reduced sperm competitiveness

To define the physiological role of the miR-506 family, we sequentially deleted these miRNA genes using CRISPR-Cas9-based genome editing (Fig. 4A) (15). We first generated the KO mice lacking either the miR-883 single cluster (miR-883 sKO) or the miR-465 single cluster (miR-465 sKO) (Fig. 4A), as these two clusters are the most abundantly expressed in the mouse testes (Fig. 3 A and Fig. S6). No discernable defects were observed, and these KO males developed normally and were fertile (Fig. 4 B and C). On the miR-883 sKO background, we further deleted the miR-741 cluster, which we termed double KO (dKO) (Fig. 4 A), but no discernable abnormalities were observed in the dKO males either (Fig. S7 B-E). On the dKO background, we next deleted either the miR-465, termed triple KO (tKO), or the miR-471 and miR-470 clusters, termed quadruple KO (quadKO) (Fig. 4 A). Lastly, we ablated the miR-465 cluster on the quadKO background, named quintuple KO (quinKO) (Fig. 4 A). To reduce potential off-target effects due to multiple rounds of CRISPR-Cas9 targeting, we also generated a KO mouse line with only 4 guide RNAs (gRNAs) flanking the SmiR region, named X-linked SmiR KO (XS) (Fig. 4 A). The XS mice were genetically equivalent to the quinKO mice, and phenotypically identical to quinKOs (Fig. 4 B and C), suggesting the phenotype observed was not due to the accumulating off-target effects. To further exclude the potential off-target effect, all KO mouse strains were backcrossed with WT C57BL/6J mice for at least 5 generations before data collection. In addition, T7 endonuclease I (T7EI) assays showed no discernable off-target effects in the quinKO mice (Fig. S7F).

Figure 4.

While the litter size was still comparable between the tKO and WT control mice, the quadKO, quinKO, and XS males produced significantly smaller litters (∼5 vs. ∼8 pups/litter) (adjusted p-value < 0.05, One-Way ANOVA) (Fig. 4 B). Of interest, no significant changes were detected in litter interval, testis weight or histology in any of the four types of KOs, as compared to WT controls (Fig. 4 C, D and G). Computer-assisted sperm analyses (CASA) revealed no significant differences in sperm counts and motility parameters among the four types of KOs (Fig. 4 E and F and Fig. S7H). Overall, there appears to be an inverse correlation (R² = 0.9139, p < 0.05, F-test) between the number of miRNAs inactivated and the litter size (Fig. S7G). Interestingly, several human studies have correlated the dysregulated miR-506 family miRNAs with impaired male fertility due to maturation arrest and oligo-asthenozoospermia (Table S4)(44–48). These data suggest that the miR-506 family may play an important role in spermatogenesis and male fertility.

Most of the protein-coding genes that are exclusively or preferentially expressed in the testis with an essential role in spermatogenesis are highly conserved across species (49). Despite their male germ cell-predominant expression, the miR-506 family miRNAs appear to have evolved rapidly to diverge their sequences, suggesting that these miRNAs might control certain “non-conserved” aspects of spermatogenesis, leading to enhanced sperm competitiveness for male reproductive success. Supporting this hypothesis, previous reports have documented that females of most species throughout the animal kingdoms mate with multiple males before pregnancy, suggesting that sperm competition may serve as a selection mechanism to bias the birth of offspring sired by the males with more competitive sperm (20, 21). Studies have also shown female rodents in the wild mate with multiple males and produce litters of mixed paternity, and that pups born to the females following such polyandrous mating display greater survival rates than those produced from females following monandrous mating (50). Given that CASA detected no difference in swimming patterns between quinKO and WT sperm (Fig. S7H), we next carried out sperm competition experiments that mimic polyandrous mating in the wild. Since the miR-506 family miRNAs are X-linked, the Y sperm from the quinKO mice are genetically indistinguishable from those of WT controls. We, therefore, adopted the mTmG male mice (51) for sperm competition experiments because the embryos or offspring fathered by the mTmG males can be easily identified based on the constitutively expressed membrane-tagged tomato red (mT) fluorescence and/or PCR genotyping.

We first conducted sequential mating with two mating events ∼6-8 hours apart. Interestingly, all of the pups born were fathered by mTmG males (n=8) when the WT females were mated first with mTmG males and subsequently with the quinKO males. In contrast, when the WT females were mated first with quinKO males and subsequently with mTmG males, ∼89% of the pups born were fathered by quinKO males, and the remaining ∼11% of pups were from mTmG males (n=28) (Fig. 4 H-I). It is noteworthy that in the sequential mating experiments, the two coituses occurred ∼6-8 hours apart due to practical reasons, whereas in the wild, polyandrous mating may take place much faster. To better mimic polyandrous mating in vitro, we mixed the WT and quinKO sperm in different ratios and used the mixed sperm to perform in vitro fertilization (IVF). MII oocytes fertilized by mTmG, quinKO, or a mixture of two types of sperm at three ratios (mTmG: quinKO = 1:1, 4:1 and 1:4) all displayed comparable rates at which fertilized oocytes developed into blastocysts (Fig. S7I). Interestingly, when a 1:1 ratio (mTmG sperm: quinKO sperm) was used, ∼73% of the resulting blastocysts were derived from mTmG sperm, whereas the remaining ∼27% were from quinKO sperm (n=179) (p < 0.0001, Chi-squared test) (Fig. 4 J). When a 4:1 sperm ratio (mTmG: quinKO) was used, ∼92% of the blastocysts were from mTmG sperm and only 8% were from quinKO sperm (n=170) (p < 0.05, Chi-squared test) (Fig. 4 J). In contrast, when a 1:4 sperm ratio (mTmG: quinKO) was used, blastocysts derived from mTmG and quinKO sperm represented ∼28% and ∼72% of the total, respectively (n=135) (Fig. 4 J). We also performed co-artificial insemination (AI) using mTmG and quinKO sperm. When a 1:1 sperm ratio (mTmG: quinKO) was used, ∼62.5% of the embryos (n=96) were derived from mTmG sperm (p< 0.05, Chi-squared test) (Fig. 4 K). When a 1:4 ratio was used, ∼35.3% of the embryos (n=17) were from the mTmG mice (Fig. 4K and Fig. S7J). Together, these results indicate that the quinKO sperm are less competitive than the control mTmG sperm both in vivo and in vitro. Previous studies suggest that sperm aggregation and midpiece size might be involved in sperm competitiveness (52, 53), but no changes in these two parameters were observed in the quinKO sperm (Fig. S7 K and L). Although the blastocyst rate (out of 2-cell embryos) of quinKO was comparable to that of the mTmG mice, the 2-cell rates (out of zygotes) were significantly reduced (p< 0.05, paired t-test) in the quinKO mice (∼39%, n=67) when compared to the mTmG mice (∼88%, n=58) (Fig. S7M), implying that the quinKO sperm is indeed less efficient in fertilizing eggs and/or supporting early embryonic development, especially the first cleavage of the zygotes.

X-linked miR-506 family miRNAs mostly target the genes involved in spermatogenesis and embryonic development and compensate for each other

To identify the target genes of these X-linked miR-506 family miRNAs, we performed RNA-seq analyses using testis samples from the five types of KOs (miR-465 sKO, dKO, tKO, quadKO, and quinKO) (Fig. S8A and Table S5). Comparisons between the KO and WT testes revealed thousands of differentially expressed genes (DEGs) (fold change ≥ 2, FDR< 0.05, Fig. S8A and Table S5). The DEGs identified were then compared with the predicted miR-506 target genes using four different databases, including TargetScan (54), microrna.org (55), miRWalk (56) and mirDB (57), to predict the differentially expressed targets (DETs) of the miR-506 family miRNAs (Fig. S8A and Table S5). We obtained 2,692, 2,028, 1,973, 3,405, and 1,106 DETs from miR-465 sKO, dKO, tKO, quadKO and quinKO testes, respectively. GO terms of DETs from each KO testis revealed that the DETs were mostly involved in embryonic development, response to stimulus, centrosome cycle, epithelium morphogenesis, organelle organization, cell projection, RNA metabolic process, and DNA repair (Fig. S8B). The 431 DETs identified to be shared across all five KO testes were also enriched in similar pathways (Fig. 5 An and B). Several genes, including Crisp1, Egr1, and Trpv4, were selected for validation using qPCR, Western blots and luciferase-based reporter assays. Consistent with the RNA-seq data, qPCR showed that Crisp1, Egr1, and Trpv4 were significantly downregulated in the quinKO testes (Fig. S8C). CRISP1 is enriched in the sperm principle piece and head (Fig. S8D). Western blots also confirmed that CRISP1 is downregulated in the quinKO testis when compared to the WT testis (Fig. S8E). Luciferase assays further confirmed that Egr1 and Crisp1 are targets of the miR-506 family members (Fig. S8 F and G). Egr1 3’UTR luciferase activity was upregulated by miR-465c, while downregulated by miR-743b (Fig. S8F). miR-465a, miR-465c, miR-470, miR-741, and miR-743a upregulated Crisp1 3’UTR luciferase activity, while miR-743b exerted the opposite effect (Fig. S8G). Of interest, KO of Crisp1 in mice or inhibition of CRISP1 in human sperm appears to phenocopy the quinKO mice (58, 59). Specifically, sperm motility in the Crisp1 KO mice is comparable to that in WT mice, but their ability to penetrate the eggs was reduced in the Crisp1 KO mice (58); a similar effect was also observed in human sperm treated with anti-hCRISP1 antibody (59).

Fig. 5.

The inverse correlation between the number of miRNAs inactivated and the severity of the phenotype strongly hints that these miRNAs compensate for each other (Fig. 4B and Fig. S7G). To test this hypothesis, we performed sRNA-seq on four KO (miR-465 sKO, tKO, quadKO and quinKO) testes. The sRNA-seq data showed that these miRNAs were no longer expressed in the corresponding KOs, confirming the successful deletion of these miRNAs in these KOs (Fig. 5C and Table S6). Interestingly, in miR-465 sKO testes, miR-201, miR-547, miR-470, miR-471, miR-742, miR-871, miR-881, miR-883a, and miR-883b were all significantly upregulated (FDR < 0.05). Similarly, miR-201, miR-547, miR-470, miR-471, miR-871, and miR-883b were all significantly upregulated in the tKO testes (FDR < 0.05); miR-201 and miR-547 were all significantly upregulated in the quadKO and the quinKO testes (FDR < 0.05) (Fig. 5 C and Table S6). These results support the notion that genetic compensation exists among the X-linked miR-506 family miRNAs.

Rapid evolution of the miR-506 family is not driven by the increased complexity of 3’UTRs of the conserved targets but rather adaptive to targeting more genes

To minimize false positives, the DETs in mice were selected using the following criteria: (1) dysregulated by fold change ≥ 2 and FDR < 0.05. (2) Falling within the predicted targets. (3) Intersected with at least two different KO mouse samples. Using the 3,043 DETs identified in mice as a reference, we searched the predicted targets of the miR-506 family miRNAs in rats and humans to determine if these target genes were shared across species (Fig. 5A). While 2,098 (∼69%) target genes were shared among all three species, 2,510 (∼82%) were common to both humans and mice, and 2,202 (∼72%) were shared between mice and rats (Fig. 6 A and Table S7). To test the accuracy of the predicted targets, we selected several genes in humans and performed luciferase assays using human miR-506 family miRNAs and their corresponding target genes (Fig. S9 A and B). Among these targets, Crisp1 and Fmr1 were shared among humans, mice, and rats, and confirmed to be targeted by the miR-506 family miRNAs in mice (Fig. S8C, S8E, S8G, and ref 11 and 15)(11, 15). Luciferase assays also confirmed that human CRISP1 (hCRISP1) and FMR1(hFMR1) were targets of the miR-506 family, and miR-510 and miR-513b both could activate hCRISP1 3’UTR luciferase activity (Fig. S9A), whereas miR-509-1, miR-509-2, miR-509-3, miR-513b, miR-514a, and miR-514b could enhance hFMR1 3’UTR luciferase activity (Fig. S9B). These results confirmed the accuracy of our predicted targets in humans.

Figure 6.

We considered two likely explanations for the paradox where the majority of their target genes were shared across species despite the rapid evolution of the miR-506 family miRNAs: 1) The 3’UTR sequences in extant target genes became increasingly divergent during evolution such that the miR-506 family miRNAs had to adapt to maintain their ability to bind these 3’UTRs; or 2) that the miR-506 family miRNAs evolved rapidly in a manner that allowed them to target mRNAs encoded by additional genes involved in spermatogenesis. To distinguish the two possibilities, we first compared the extent of similarities among the 3’UTR sequences of the 2,510 shared target genes between humans and mice (Fig. 6 A and Table S7). We adopted the PhyloP scores to measure the evolutionary conservation at individual nucleotide sites in the 3’UTRs of the shared target genes. The overall conservation appeared to be greater in the regions targeted by miR-506 family miRNAs than in the non-target regions in the 3’UTRs of the shared target genes in both mice and humans (Fig. 6 B and C) with few exceptions (Fig. 6 D and E) (p < 0.05, t-test). These data suggest that the regions targeted by the X-linked miR-506 family miRNAs are under relatively stronger purifying, rather than adaptive, selection. We then tested the second hypothesis that the rapid evolution of the miR-506 family resulted in more extant mRNAs being targeted by these miRNAs. We first compared the average target numbers of each miR-506 family miRNA between humans and mice using the 2,510 shared targets between predicted targets in humans and the dysregulated targets in mice (Fig. 6 A and F). Among these shared targets, the human miR-506 family members could target ∼1,268 unique transcripts per miRNA, whereas the murine miR-506 family members could only target ∼1,068 (p < 0.05, t-test) (Fig. 6 F), indicating that the miR-506 family miRNAs target more genes in humans than in mice. Furthermore, we analyzed the number of all potential targets of the miR-506 family miRNAs predicted by the aforementioned four algorithms among humans, mice, and rats. The total number of targets for all the X-linked miR-506 family miRNAs among different species did not show significant enrichment in humans (Fig. S9C), suggesting the sheer number of target genes does not increase in humans. We then compared the number of target genes per miRNA. When comparing the number of target genes per miRNA for all the miRNAs (baseline) between humans and mice, we found that on a per miRNA basis, human miRNAs have more targets than murine miRNAs (p<0.05, t-test) (Fig. S9D), consistent with higher biological complexity in humans. This became even more obvious for the X-linked miR-506 family (p<0.05, t-test) (Fig. S9D). In humans, the X-linked miR-506 family, on a per miRNA basis, targets a significantly greater number of genes than the average of all miRNAs combined (p<0.05, t-test) (Fig. S9D). In contrast, in mice, we observed no significant difference in the number of targets per miRNA between X-linked miRNAs and all of the mouse miRNAs combined (mouse baseline) (Fig. S9D).

These results suggest that although the sheer number of target genes remains the same between humans and mice, the human X-linked miR-506 family targets a greater number of genes than the murine counterpart on a per miRNA basis. We also investigated the number of miR-506 family miRNA targeting sites within the individual target genes in both humans and mice, but no significant differences were found between humans and mice (Fig. 6 G). To determine whether increased target sites in humans were due to the expansion of the MER91C DNA transposon, we analyzed the MER91C DNA transposon-containing transcripts and associated them with our DETs. Of interest, 28 human and 3 mouse mRNAs possess 3’UTRs containing MER91C DNA sequences, and only 3 and 0 out of those 28 and 3 genes belonged to DETs in humans and mice, respectively (Fig. S9E), suggesting a minimal effect of MER91C DNA transposon expansion on the number of target sites. Taken together, these results suggest the human X-linked miR-506 family has been subjected to additional selective pressure, causing them to exert additional regulatory functions by targeting more mRNAs expressed during spermatogenesis (Fig. 6 H).

Discussion

Successful reproduction is pivotal for the perpetuation of species, and sperm are constantly facing selective pressures (60). To enhance their chance to fertilize eggs, sperm need to adapt accordingly, and miRNAs-mediated regulation of gene expression in spermatogenesis provides a rapidly adaptable mechanism toward this end. Although miRNAs were initially believed to be evolutionarily conserved, the number of non-conserved miRNAs has been steadily increasing (61). Among the non-conserved miRNAs, many are derived from TEs, suggesting that TEs may serve as a major source of miRNA sequences (61). The delayed recognition of TE-derived miRNAs, in part, results from the fact that repetitive sequences were usually excluded during the computational annotation of miRNAs. In theory, TEs can serve as a good donor of miRNA sequences for the following reasons: 1) TEs are ubiquitous and abundant in the genome and are known to contribute to the regulatory elements of the coding genes, e.g., UTRs (62). As one of the regulatory factors that target mainly the 3’UTRs, TE-derived miRNAs can regulate a larger number of mRNAs with multiple miRNA-targeting sites. 2) TEs are among the most rapidly evolving sequences in the genome (63) and thus, can continuously produce species/lineage-specific miRNA genes to diversify their regulatory effects.

Consistent with these notions, the present study provides evidence supporting that the miR-506 family miRNAs originated from the MER91C DNA transposons. Of more interest, the miR-506 family miRNAs, despite their rapid evolution, are all expressed in spermatogenic cells in the testis and sperm, supporting a lineage-specific functional diversification of TE-derived miRNAs. In fact, the miR-506 family miRNAs were among the first reported TE-derived miRNAs because of their location in a small region of the X chromosome and their confined, abundant expression in the testis (16). RNAs are much less abundant in sperm than in somatic or spermatogenic cells (∼1/100) (64). Sperm-borne small RNAs represent a small fraction of total small RNAs expressed in their precursor spermatogenic cells, including spermatocytes and spermatids (64). Therefore, when the same amount of total/small RNAs are used for quantitative analyses, sperm-borne small RNAs (e.g., miR-506 family miRNAs) would be proportionally enriched in sperm compared to other spermatogenic cells.

It has been demonstrated that under certain circumstances, genes that evolve under sexual conflicts tend to move to the X chromosome, especially when they are male-beneficial, female-deleterious, and act recessively (65, 66). The X chromosome is enriched with genes associated with male reproduction (15, 16, 67). The rapid evolution of the X-linked miR-506 family strongly suggests that these miRNA genes were under selection to expand and diversify their regulatory effects on spermatogenesis. Indeed, our data strongly suggest that targeting new mRNAs was likely the driving force for the rapid evolution of the miR-506 family of miRNAs. However, expansion and sequence divergence of the X-linked miR-506 family may simply reflect natural drifting without functional significance, similar to some of the pachytene piRNA clusters (30). We argue that the neutral drifting theory may not be true to the miR-506 family for the following reasons: 1) Despite highly divergent overall sequences of the miR-506 family, some miRNAs share the same seed regions across multiple species, suggesting that these regions may undergo strong selections. 2) The miR-506 family miRNAs, especially the FmiRs, are highly conserved in modern humans, implying a strong selection of these miRNAs. 3) Knockout the miR-506 family (either quinKO or XS) results in male subfertility, reflecting a biological function. 4) The quinKO sperm are less competitive than the WT sperm both in vivo and in vitro. 5) Several human studies have linked the dysregulation of the miR-506 family with male infertility/subfertility (44–48). Similarly, one study in Drosophila also showed that the rapidly evolving testis-restricted miRNAs underwent adaptive evolution rather than neutral drifting (68).

Since TEs are abundant in UTRs, TE-derived miRNAs can target a much greater number of mRNAs than those derived from distinct non-repetitive genomic loci (61). Indeed, thousands of the dysregulated genes detected in the miR-465 sKO, dKO, tKO, quadKO and quinKO testes are involved in multiple pathways of spermatogenesis. By analyzing the sequence divergence, we noticed that the most common sequence substitutions among all of the miR-506 miRNA sequences were U-to-C and A-to-G, which were likely mediated by ADARs (adenosine deaminases acting on RNA) that can change A to I (which is functionally equivalent to G) (69). Interestingly, ∼90% of the A-to-I editing appears to have occurred in Alu elements (belonging to the SINE family), and some of the edits occurred in miRNAs (69). Since G and U can form the so-called G-U wobble base pair (70), those U-to-C or A-to-G substitutions can, in theory, target similar sequences and exert regulatory functions (71), suggesting that the evolving miRNA sequences could target not only the original sequences but also new sites with similar sequences. Consistent with this notion, the predicted target genes of the miR-506 family in mice can also be found in rats and humans, suggesting the target genes are shared across species despite the quick divergence of the miRNA sequences across species. This is also supported by our data showing that the binding sites for the miR-506 family of miRNAs are more conserved than the surrounding, non-targeting regions in the 3’UTRs of the predicted target mRNAs. Furthermore, seed sequences among some miR-506 family miRNAs remain the same despite the high divergence of these miRNAs, and these conserved seed sequences appear to be present in the dominant mature miRNAs. Thus, the seed region of these miRNAs appears to have undergone strong selection. Supporting this notion, previous studies have shown correlations between miRNA expression and the evolution of miRNAs and target sites (72, 73). In general, miRNAs repress their target gene expression. However, numerous studies have also shown that some miRNAs, such as human miR-369-3, Let-7, and miR-373, mouse miR-34/449 and the miR-506 family, and the synthetic miRNA miRcxcr4, activate gene expression both in vitro (74, 75) and in vivo (11, 15, 76, 77). Earlier reports have shown that these miRNAs can upregulate their target gene expression, either by recruiting FXR1, targeting promoters, or sequestering RNA subcellular locations (11, 74, 75). Of interest, miRNAs with the same seed sequences may exert divergent functions. For example, the mature miR-465a, miR-465b, and miR-465c only have a few mismatches outside of the seed region, but only miR-465c exerts functional activation of Egr1. A similar effect has also been reported in the miR-465 cluster on the Alkbh1 3’UTR activity (18). Similarly, despite the same seed sequences in the miR-465 or miR-743 cluster, miR-465a and miR-465c have differential activating effects on the 3’UTR of Crisp1, and miR-743a and miR-743b exert opposite effects on the Crisp1 3’UTR, further confirming their functional divergence. Therefore, the sequences outside of the miRNA seed region may play an important role in their functions, which have also been observed in C. elegans and human HEK293 cells (18, 78, 79). To unequivocally demonstrate the physiological role of miRNAs, it would be ideal to delete not only the miRNAs but also their binding sites in their target transcripts in vivo. A few previous studies have established the miRNA: target relationship by deleting the miRNA-binding sites in target transcripts in C. elegans, Drosophila, and cell lines (80–82), but similar studies have not been reported in mice or humans. The strategy may work in mRNAs with 3’UTRs containing only one or two miRNA-binding sites, but for more complex 3’UTRs of mRNAs in mice and humans that often contain multiple binding sites for the same or different miRNAs, deletion of one miRNA binding site may not cause any discernable effects as the loss of function can easily be compensated by other miRNAs. Nevertheless, by deleting all five highly expressed clusters of the miR-506 family one by one, we were able to overcome the compensatory effects among the family members/clusters and successfully revealed the physiological role of this miRNA family. Based on small RNA-seq, some FmiRs, e.g., miR-201 and miR-547, were upregulated in the SmiRs KO mice, suggesting that this small cluster may act in concert with the other 5 clusters and thus, worth further investigation.

It is well-known that mice in the wild are promiscuous, and one female often mates with multiple males sequentially, giving rise to polyandrous litters derived from sperm from more than one sire (20, 21). Polyandrous mating establishes a situation where sperm from multiple males co-exist in the female reproductive tract, with the most competitive ones fertilizing eggs and producing offspring (22, 83).

Therefore, a male that may be fertile in the monandrous mating scheme may rarely sire offspring in a polyandrous mating scenario, rendering this male functionally “sub-fertile” or even “infertile.” Therefore, sperm competitiveness reflects the general reproductive fitness of the male (83). Although the quinKO males tend to produce smaller litters under the monandrous mating scheme, their sperm counts, sperm motility and morphology are indistinguishable from those of WT sperm. This is not surprising given that miRNAs of the miR-506 family most likely function to control certain non-essential aspects of spermatogenesis. Sperm can be subject to competition at multiple steps during fertilization, including their migration through the female reproductive tract (cervix, uterine cavity and oviduct), binding the cumulus-oocyte complexes, penetration of zona pellucida, etc. Therefore, IVF may not be ideal for evaluating sperm competition in the real world as it bypasses several key sites where sperm competition likely takes place. Artificial insemination (AI) may represent a better way to assess sperm competition than IVF, but it is probably less desirable than polyandrous mating for the following reasons: First, in the wild, sperm from two males rarely, if not never, enter the female reproductive tract simultaneously. We had tried to place two males into the cage with one female, but the two males ended up fighting, and the submissive one never mated. Second, sperm are delivered directly into the uterus or oviduct during AI (84, 85), thus bypassing the potential sites for sperm competition (e.g., cervix and uterine cavity). Although our breeding scheme also involves sperm competition, by shortening the time between the two mating events in a laboratory setting, the sequential mating method reported here may be further improved to better mimic the natural polyandrous mating in the future. Moreover, future analyses of the quinKO sperm may help identify biochemical or molecular biomarkers for sperm competitiveness.

In summary, our data suggest that the miR-506 family miRNAs are derived from the MER91C DNA transposon. These miRNAs share many of their targets and can compensate for each other’s absence, and they work jointly through regulating their target genes in spermatogenesis to ensure sperm competitiveness and male reproductive fitness.

Materials and Methods

Animal care and use

All mice used in this study were on C57BL/6J background and housed in a temperature-and humidity-controlled, specific pathogen-free facility under a light-dark cycle (12:12 light-dark) with food and water ad libitum. Animal use protocol was approved by the Institutional Animal Care and Use Committees (IACUC) of the University of Nevada, Reno and The Lundquist Institute at Harbor-UCLA, and is in accordance with the “Guide for the Care and Use of Experimental Animals” established by the National Institutes of Health (1996, revised 2011).

Generation of the knockout mice

The single, double, triple, quadruple, and quintuple miR-506 KO mice were generated as previously described (15). Briefly, Cas9 mRNA (200 ng/μl) and gRNAs flanking the miR-506 family subclusters (100 ng/μl) were mixed and injected into the cytoplasm of zygotes in the M2 medium. After injection, all embryos were cultured for 1 h in KSOM+AA medium (cat. #MR-121-D, Millipore) at 37°C under 5% CO₂ in the air before being transferred into 7–10-week-old female CD1 recipients. The miR-883 sKO or the miR-465 sKO mice were first generated. After two rounds of backcrossing with C57BL/6J mice, the miR-741 cluster was knocked out on the miR-883 sKO background to generate dKO mice. After two rounds of backcrossing of the dKO with C57BL/6J mice, the miR-465 cluster or miR-471 and miR-470 clusters were further deleted, thereby generating tKO or quadKO mice. Lastly, the miR-465 cluster was ablated on the quadKO background, leading to the production of quinKO mice. The XS mice were generated by using only four gRNAs flanking the SmiR region on the C57BL/6J background. All KO mice were backcrossed with the C57BL/6J mice for at least five generations before collecting data.

Sequential polyandrous mating

Sequential polyandrous mating was carried out based on the ovulation time point (10-13h after hCG) as previously described (84). Adult (8-12 weeks of age) C57BL/6J females were injected (i.p.) with 7 IU PMSG at 8 p.m., followed by 7 IU hCG 48h later. After hCG, the first male mouse was put into the cage of one female from 8 p.m. to 6 a.m. the next day. The first plug was marked with a marker pen. A second male mouse was then introduced into the cage of the plugged female, which was checked every 30-40 min to identify a new plug (non-marked). Females that were plugged twice were kept for producing pups for paternity analyses.

In vitro fertilization (IVF)

Adult (8-12 weeks) C57BL/6J female mice were first treated with 7 IU pregnant mare serum gonadotropin (PMSG, Cat.# HOR-272, Prospecbio) through intraperitoneal (i.p.) injection followed by i.p. injection of 7 IU hCG 48h later. Oocytes were collected from the ampulla ∼14 h after the hCG (Cat.# HOR-250, Prospecbio) treatment, and the cumulus cells surrounding oocytes were removed by treatment with bovine testicular hyaluronidase (1.5 mg/ml; Cat.# H3506, Sigma,) in M2 (Cat.# MR-015-D, Millipore) at 37°C for 2 min. The cumulus-free oocytes were washed and kept in equilibrated HTF (Cat# MR-070-D, Millipore) at a density of 20-30 oocytes per 60 µl HTF at 37°C in an incubator with air containing 5% CO₂ prior to IVF. Cauda epididymal sperm were collected in 100 µl of equilibrated HTF medium, allowing spermatozoa to capacitate for ∼30 min at 37°C in an incubator containing 5% CO2 air. After capacitation, spermatozoa (2 µl) were diluted by 10 folds and subjected to computer-assisted sperm analyses (CASA) using the Sperm Analyzer Mouse Traxx (Hamilton-Thorne). Based on the sperm concentration, an aliquot of 2.5 x 10⁸ spermatozoa was added into each HTF-oocytes drop (∼60µl) for IVF. ∼4 h later, zygotes were washed and cultured in KSOM+AA (Cat.# MR-121-D, Millipore) until the blastocyst stage at 37°C in an incubator with air containing 5% CO₂. The 2-cell embryos were counted 24-26 h after IVF, and blastocysts were counted and analyzed under a fluorescence microscope 70-72 h after IVF.

Artificial insemination (AI)

At least 2-month-old female CD1 or C57BL/6J mice were administrated with 2.5 IU of PMSG (Cat.# HOR-272, Prospecbio) at 5:30 PM 3 days before artificial insemination, followed by 2.5 IU of hCG (Cat.# HOR-250, Prospecbio) at 5:00 PM 1 day prior to AI. The next morning at 8:00 AM, ∼2-month-old mTmG and quinKO male mice were sacrificed, the cauda epididymis was dissected, and fat tissue and blood were removed before placing the cauda epididymis into 500μl or 150μl of EmbryoMax® Human Tubal Fluid (HTF) (1X) (Cat.# MR-070-D, Millipore Sigma) containing 4mg/ml BSA (Cat.# 12659-250GM, EMD Millipore Corp.) (HTF-BSA) covered with 4 ml of mineral oil (Cat.# M8410-500ML, Sigma). Three incisions were made on the cauda to allow sperm to swim out and to get capacitated for at least 30 min. 25 μl of mTmG and 25 μl quinKO sperm suspensions were mixed, and 40μl (if using 500 μl HTF-BSA) or 25μl (if using 150 μl HTF-BSA) of the mixed sperm were delivered to superovulated females using C&I Device for Mice (Cat.# 60020, Paratech) at 9:00 AM. Recipients were immediately paired with vasectomized males overnight. The next day, the plug was checked and the female mice with plugs were used for collecting embryonic day 10 (E10) embryos, and the ones without plugs were used for collecting zygotes, 2-cell embryos, morulae, or blastocysts embryos.

Mouse genotyping

Mouse tail snips were lysed in a lysis buffer (40mM NaOH, 0.2mM EDTA, pH=12) for 1h at 95°C, followed by neutralization with the same volume of neutralizing buffer (40mM Tris-HCl, pH 5.0). PCR reactions were conducted using the 2×GoTaq Green master mix (Promega, Cat. M7123). The primers used for genotyping are the same as previously described (15). For single embryo genotyping (e.g. zygotes, 2-cell embryos, 4-cell embryos, morulae, and blastocysts), each embryo was picked up by mouse pipetting and transferred into a 200μl tube, and lysed in 10 μl of lysis buffer (100 mM Tris–HCl (pH 8.0), 100 mM KCl, 0.02% gelatin, 0.45% Tween 20, 60 μg/mL yeast tRNA, and 125 μg/mL proteinase K) at 55 °C for 30min followed by inactivation at 95 °C for 10min. 2μl of the lysis was used as the template for the first round of PCR (30 cycles) in a 10μl reaction using the PrimeSTAR® HS DNA Polymerase (Cat.# R010B, Takara) or 2×GoTaq Green master mix (Cat.# M7123, Promega). Then 2μl of the first PCR was used for the second round of PCR in a 10μl reaction using 2×GoTaq Green master mix (Cat.# M7123, Promega) for 35 cycles. Primers used for embryo genotyping were included in Table S8.

openCASA

Sperm parameters were assessed using openCASA (86). After sperm capacitation in HTF containing 4%BSA at 37 °C for 30min, the video was recorded as an AVI format at 60 frames per second (FPS) for 2 seconds with a resolution of 768*576 pixels using UPlan FL N 4×/0.13 PhP Objective Lens (Olympus) and DMK 33UP1300 camera (The imaging source). Motility module in openCASA was set with the following parameters: 1.21 Microns per Pixel, the cell size of 10-200 μm^2, Progressive motility (STR>50%, VAP>50), Minimum VCL of 10 μm/s, VCL threshold of 30-200 μm/s, 60 Frame Rate (frames/S), 10 Minimum Track Length (frames), 20 μm Maximum displacement between frames, and Window Size (frames) of 4.

Analysis of conservation of miR-506 family in modern humans using the 1000 Genomes Project (1kgp)

The vcf files from the 1000 Genomes Project covering 3202 samples were downloaded from http://ftp.1000genomes.ebi.ac.uk/vol1/ftp/data_collections/1000G_2504_high_coverage/working/20220422_3202_phased_SNV_INDEL_SV/. The miRNA annotations were obtained from UCSC genome browser, and pachytene piRNA hg19 genome coordinance was obtained from (30) and converted to GRCh38 genome coordinance using LiftOver. The derived allele frequency (DAF) was retrieved from the vcf file, and mean nucleotide diversity (MND) was calculated as 2*DAF*(1-DAF) as previously described (30). Kruskal-Wallis test was used for statistical analysis, and adjusted p < 0.05 was identified as statistically significant.

Overexpression of MER91C

RNA structures for MER91C were predicted using RNAfold (87). MER91C DNA transposons from humans (hsa-miR-513a1), dogs (cfa-miR-507b), and horses (eca-miR-514a) were synthesized by IDT and inserted into pCI-Neo plasmid (Cat.# E1841, Promega) using NEBuilder® HiFi DNA Assembly Master Mix (Cat.# E2621L, NEB). 150ng of pCI-Neo (negative control) or pCI-Neo-Mer91c were transfected with or without 150ng of pcDNA3.1+_FH-AGO2-WT (Plasmid # 92006, Addgene) into HEK293T cells at ∼60% confluency in 24-well plates. 24h later, cells were harvested followed by RNA extraction using mirVana™ miRNA Isolation Kit (Cat# AM1561, Thermo Fisher Scientific), polyadenylation by E. coli Poly(A) Polymerase (Cat.# M0276L, NEB), reverse transcription using SuperScript IV Reverse Transcriptase (Cat.# 18090010, Thermo Fisher Scientific), and PCR using 2×GoTaq Green master mix (Cat.# M7123, Promega) or qPCR by PowerUp™ SYBR™ Green Master Mix (Cat.# A25742, Thermo Fisher Scientific). Primers used for RT-PCR, PCR, and qPCR were included in Table S8.

Bioinformatic analyses of transposable element (TE)

Genomic regions and GFF3 files for transcript 3’UTR annotations were downloaded from the UCSC genome browser, and GTF files for transposon annotations were downloaded from https://labshare.cshl.edu/shares/mhammelllab/www-data/TEtranscripts/TE_GTF/. For the TE containing transcripts, bedtools closest was used to extract the closest TEs to transcripts. Genomes used were GRCh38 (humans), panTro6 (chimpanzees), equCab3 (horses), GRCm39 (mice), and rn7 (rats).

Phylogenetic tree analysis of the MER91C DNA transposon and the miR-506 family miRNAs

The transposon fasta sequences from humans, dogs, horses, and guinea pigs were downloaded from the UCSC genome browser and aligned to the miR-506 family miRNAs in their corresponding species using BLAST (Basic Local Alignment Search Tool) (34). After retrieving the transposons that aligned to the miR-506 family miRNAs, the miR-506 family miRNAs and the transposons were aligned using ClustalW2 followed by phylogenetic tree building using IQ-TREE2 with default parameters (88). The final figure was generated using Geneious software.

Purification of germ cells

Pachytene spermatocytes and round spermatids were purified from adult C57BL/6J mice using the STA-PUT method. BSA gradients (2-4%) were prepared in EKRB buffer with a pH of 7.2 containing 1× Krebs-Ringer Bicarbonate Buffer (Cat.# K4002, Sigma), 1.26 g/L sodium bicarbonate (Cat.# S6761, Sigma), 1× GlutaMAX (Cat.# 35050061, Thermo Fisher Scientific), 1× Antibiotic-Antimycotic (Cat.# 15240062, Thermo Fisher Scientific), 1× MEM Non-Essential Amino Acids (Cat.# 11140050, Thermo Fisher Scientific), 1×MEM Amino Acids (Cat.# 11130051, Thermo Fisher Scientific), and 100 ng/ml cycloheximide (Cat.# 01810, Sigma). After being removed and decapsulated, testes were placed into 10 ml of EKRB buffer containing 0.5 mg/ml type IV collagenase (Cat.# C5138, Sigma) and digested at 33 °C for ∼12min to dissociate the seminiferous tubules. Once dissociated, the seminiferous tubules were washed three times using EKRB buffer to remove the interstitial cells and red blood cells followed by trypsin digestion by incubation at 33 °C for ∼12 min with occasional pipetting in 10 ml EKRB buffer containing 0.25 mg/ml trypsin (Cat.# T9935, Sigma) and 20 μg/ml DNase I (Cat.# DN25, Sigma). 1 mL of 4% BSA-EKRB was added to the 10 mL fully dispersed testicular cells to neutralize the trypsin digestion followed by centrifuge at 800g for 5min at 4°C. Testicular cells were washed two times with EKRB buffer and resuspended in 10 mL 0.5% BSA-EKRB. The cell suspension was passed through 70 μm cell strainer (Cat.# 431751, Corning) and loaded onto the STA-PUT apparatus containing 2-4% BSA-EKRB gradients for sedimentation. After 2-3 h sedimentation at 4 °C, cell fractions were collected from the bottom of the sedimentation chamber. Fractions containing the same cell types were pooled and saved for RNA sequencing.

Library construction and RNA-seq

RNA was extracted using the mirVana™ miRNA Isolation Kit (Cat.# AM1561, Thermo Fisher Scientific) following the manufacturer’s instructions. Large RNA (>200nt) and small RNA (<200nt) were isolated separately for library construction. Small RNA libraries were constructed using NEBNext® Small RNA Library Prep Set for Illumina® (Multiplex Compatible) (Cat.# E7330L, NEB) following the manufacturer’s instructions and sequenced using HiSeq 2500 system for single-end 50bp sequencing. Large RNA libraries were constructed using the KAPA Stranded RNA-Seq Kit with RiboErase (Cat.# KK8483, Roche) and the adaptor from NEBNext® Multiplex Oligos for Illumina® (Index Primers Set 1, Cat.# E7335L, NEB). The indexed large RNA libraries were sequenced using Nextseq 500 with paired-end 75bp sequencing.

Large and small RNA-Seq data analysis

For the large RNA-seq data, raw sequences were trimmed by Trimmomatic (89), followed by alignment using Hisat2 (90), and assembly using StringTie (90). Reads were summarized using featureCounts (91) and the differential gene expression was compared using DESeq2 (92). For each KO mouse sample, the genes with a fold change ≥ 2 and FDR < 0.05 were considered DEGs. The DEGs in each KO mouse were then intersected with the corresponding miRNA targets predicted by four different algorithms, including TargetScan (54), microrna.org (55), miRWalk (56), and mirDB (57). The gene is considered a putative target (differentially expressed target, DETs) as long as it intersects with the targets predicted by any method mentioned above. The DETs identified in each KO mouse sample were intersected with other KO mouse samples, and the DETs intersected at least two different KO mouse samples were selected as the “pool” of DETs in mice. Then the DETs “pool” in mice was intersected with miR-506 family predicted targets in humans and rats to determine the shared targets among humans, rats, and mice. A mRNA that any miR-506 family member is targeting is deemed as the shared target.

For the small RNA-seq data, we applied the AASRA (93) pipeline (for mice) or SPORTS1.0 (94) (for humans, monkeys, rats, and horses) to parse the raw sequencing data. The clean reads were mapped against miRbase (28). The DESeq2 (92) (for mice) or edgeR (95) (for humans, monkeys, rats, and horses) algorithm was used to compare the groupwise miRNA expression levels. The RNAs with a false discovery rate (FDR) < 5% were deemed differentially expressed. Cohen’s d was computed by the “cohensD” function within the “lsr” package.

Luciferase assay

For luciferase reporter assays, the 3′ UTR of Crisp1, Egr1, hCRISP1, and hFMR1 were amplified using C57BL/6J tail snips or HEK293 cells genomic DNA template with Q5® Hot Start High-Fidelity 2X Master Mix (Cat.# M0494L, NEB). The PCR products were inserted into psiCHECK-2 vector (Cat.# C8021, Promega) via Xho I (Cat.# R0146S, NEB) and Not I (Cat.# R3189S, NEB) restriction enzymes cutting sites downstream of the Renilla luciferase-coding sequence using either NEBuilder® HiFi DNA Assembly Master Mix (Cat.# E2621L, NEB) or T4 DNA ligase (Cat.# M0202L, NEB). For the miRNA overexpression plasmids, ∼300bp upstream and downstream of the precursor miRNA genomic region were amplified using C57BL/6J tail snips or HEK293 cells genomic DNA with Q5® Hot Start High-Fidelity 2X Master Mix (Cat.# M0494L, NEB) or PrimeSTAR® HS DNA Polymerase (Cat.# R010B, Takara), and inserted into pcDNA3.1 plasmids using NEBuilder® HiFi DNA Assembly Master Mix (Cat.# E2621L, NEB). HEK293T cells were co-transfected with 150ng pcDNA3.1-miRNA and 150 ng psiCHECK-2 containing the 3′ UTR of the target gene using Lipofectamine 3000 (Cat.# L3000015, Thermofisher Scientific) in a 24-well cell culture plate (Cat.# 3524, Corning) at ∼60% confluency. After 24 h of culture, cells were lysed and assayed with Dual Luciferase Assays (Cat.# E1910, Promega) according to the manufacturer’s instructions. Renilla luciferase signals were normalized to Firefly luciferase signals to adjust the transfection efficiency. pcDNA3.1-cel-miR-67 which has a minimal sequence identity to the miRNAs in humans, mice and rats was used as a negative control miRNA. Primers used for generating plasmids containing miRNAs or 3’UTR of the target genes were included in Table S8.

Immunofluorescence

Cauda sperm were capacitated in HTF at 37 °C for half an hour followed by spreading onto Superfrost Plus slides (Cat.# 22-037-246, Thermo Fisher Scientific). The slides were air-dried, fixed in 4% paraformaldehyde for 15 min (Cat.# J19943-K2, Thermo Fisher Scientific), then washed twice in 0.4% Photo-Flo 200 (Cat.# 1464510, Kodak) /1×PBS (5 min/wash), followed by a 5 min wash in 0.4% Photo-Flo 200/ddH₂O, and stored in –80 °C after air-dried. The slides were equilibrated to room temperature before immunofluorescence, followed by incubation in acetone for 20 min at 4°C and rehydration in 95% ethanol twice (5 min/wash), 70% ethanol twice (5 min/wash) and 1×PBS for 3 times (5 min/wash) sequentially. Heat-induced antigen retrieval was performed in citrate buffer (pH 6.0) with high power for 4 min once, and three times with low power for 4 min in microwave. Slides were cooled down to room temperature and washed in 1×PBS twice (5 min/wash). Following permeabilization with 0.25% Triton X-100 (Sigma-Aldrich, Cat.# T8787) in 1×PBS for 20 min at room temperature, the slides were washed with 1× PBS three times (5 min/wash), and incubated at 3% H₂O₂ solution to block endogenous peroxidase activity. After washing in 1× PBS twice (5 min/wash), the slides were blocked with 1× blocking solution (5% normal donkey serum, 5% fetal bovine serum, and 1% bovine serum albumin in 1×PBS) at room temperature for 1h, then incubated with the anti-CRISP-1 antibody (Cat.# AF4675-SP, R&D systems, 1: 100 in 1× blocking solution) at 4 °C overnight. After primary antibody incubation, the slides were washed in 1× PBS three times (10 min/wash), followed by incubation in the Donkey Anti-Goat IgG H&L (HRP) (Cat.# ab97110, Abcam, 1: 250 in 1× blocking solution) at room temperature for 1 h, and 3 times washes in 1× PBS (10 min/wash). Tyramide signal amplification was performed and stopped using reagents from Invitrogen™ Alexa Fluor™ 488 Tyramide SuperBoost™ Kit (Cat. # B40941, Thermo Fisher Scientific), followed by mounting and counterstaining in Antifade Mounting Medium with DAPI (Cat.# H-1800, Vector Lab). Nail polish was applied on the edge of the coverslips after 2h of mounting to prevent further evaporation and stored at 4°C before taking images. Images were taken using the Nikon ECLIPSE Ti2 Confocal microscope with the NIS-Elements Software.

Western blot

Testes from adult WT and KO mice were collected and sonicated in 2 × Laemmli buffer (Cat.# 1610737, Bio-Rad) supplemented with 2-Mercaptoethanol (Cat.# M6250, Sigma-Aldrich) and cOmplete™, Mini, EDTA-free Protease Inhibitor Cocktail (Cat.# 11836170001, Sigma-Aldrich) followed by incubating at 100°C for 10 min. The proteins were separated on 4–20% Mini-PROTEAN® TGX™ Precast Protein Gels (Cat.# 4561094, Bio-Rad) and then transferred onto Amersham™ Protran® Premium Western blotting membranes, nitrocellulose (Cat.# GE10600003, Sigma-Aldrich). The membranes were then stained with Ponceau S solution (Cat.# P7170, Sigma-Aldrich) to check the samples’ loading. After taking pictures, the membrane was de-stained with 0.1 M NaOH, and washed with water and TBS. Then the membrane was blocked with 5% skim milk in TBST (TBS containing 0.1% (v/v) Tween-20) for 1 h at room temperature and incubated with the anti-CRISP-1 antibody (Cat.# AF4675-SP, R&D systems, 1: 2000 in TBST containing 5% skim milk) and anti-GAPDH antibody (Cat.# G9545, Sigma, 1: 6000 in TBST) overnight at 4°C. After washing with TBST three times, the membrane was incubated with the Donkey Anti-Goat IgG H&L (HRP) (Cat.# ab97110, Abcam) or Goat Anti-Rabbit IgG H&L (HRP) (Cat.# ab6721, Abcam) at room temperature for 1 h. Followed by three washes with TBST, the bands were detected using the WesternBright ECL kit (Cat.# K-12045-D20, Advansta).

T7 endonuclease I (T7EI) assay

The potential off-target sites that may be induced by CRISPR-Cas9 were predicted using Alt-R Custom Cas9 crRNA Design Tool (IDT) and assessed by T7 endonuclease I (Cat.# M0302L, NEB) assay. The sequences were retrieved from the UCSC genome browser, and the primers flanking the off-target sites were designed to cover ∼600bp. Genomic DNA from WT C57BL/6J or quinKO was amplified using Q5® Hot Start High-Fidelity 2X Master Mix (Cat.# M0494L, NEB) or PrimeSTAR® HS DNA Polymerase (Cat.# R010B, Takara) with the designed off-target primers. 2μl of the unpurified PCR product was diluted in 1X NEBuffer 2 in a 9.5 μl volume and denatured at 95°C for 5min, followed by annealing at 95-85°C at a –2°C/second rate, and 85-25°C at a –0.1°C/second rate. Then 0.5 μl of T7EI was added to the 9.5 μl denatured PCR products and incubated at 37°C for 30min. The T7EI-treated PCR products were run on 1X TBE gel and stained with SYBR™ Gold Nucleic Acid Gel Stain (Cat.# S11494, Thermo Fisher Scientific) to detect the off-target effects. Primers used for T7EI were included in Table S8.

miRNA and 3’ UTR conservation analysis

The Multiz Alignment and Conservation method was used to measure miRNA sequence conservation with either the human or the mice genome as the reference (25, 96). 100 and 60 species were used for the human and mouse references, respectively. PhastCons takes the flanking sequences into consideration and doesn’t rely on fixed sliding windows; consequently, both highly conserved short sequences and moderately conserved long sequences can yield higher scores (97). PhastCons gives a value between 0 to 1, the higher the value is, the more conserved the region is. By contrast, PhyloP compares the conservation of individual nucleotides among all phylogeny clades, giving positive scores once the region is conserved and vice versa. PhyloP and PhastCons scores of all miRNAs, miR-506 family, Fmr1 CDS, Slitrk2 CDS, and the intergenic region (IGR) were retrieved from the UCSC genome browser and quantified. Kruskal-Wallis test was used for statistical analysis, and adjusted p < 0.05 was identified as statistically significant. PhyloP scheme was used to measure the evolutionary conservation level at individual nucleotide sites in the 3’UTRs of the target genes. Positive PhyloP scores suggest higher conservation and stronger purifying selection, whereas negative PhyloP scores indicate accelerated evolution and potential adaptive selection. The genomic annotations and mRNA sequences were based on the hg38 (25) (‘TxDb.Hsapiens.UCSC.hg38.knownGene’ and ‘BSgenome.Hsapiens.UCSC.hg38’) and mm10 (‘TxDb.Mmusculus.UCSC.mm10.knownGene’ and ‘BSgenome.Mmusculus.UCSC.mm10’) assemblies for human and mouse, respectively. The PhyloP scores were mapped to individual nucleotides in 3’UTRs based on the transcript coordinates (98).

Data Availability

The sRNA-seq and RNA-seq datasets have been deposited into the SRA database with accession#: PRJNA558973 and PRJNA670945. The scripts for sRNA analysis can be found on GitHub (https://github.com/biogramming/AASRA and https://github.com/junchaoshi/sports1.1).

Statistical analyses

Data are presented as mean ± SEM, and statistical differences between datasets were assessed by two samples t-test, F-test, Dunnett’s multiple comparisons test as the posthoc test following one-way ANOVA, Chi-squared test, or Kruskal-Wallis test as described in the text or figure legends. Normal distribution was assessed by quantile-quantile (QQ) plot or density plot. p < 0.05, 0.01, 0.001, and 0.0001 are considered as statistically significant and indicated with *, **, ***, and **** respectively.

Author Contributions

Z.W. and W.Y. designed the research. Z.W., Y. W., S. C., D.M., S. W., H. W., H.M., J. R. M., and H.Z., performed bench experiments. Z.W., T.Z., R.D.M.M., M.L., Y.X., J.M., and E. C. L. performed bioinformatics analysis. Z.W., and W.Y. wrote the manuscript.

Competing Interest Statement

The authors declare no conflict of interest.

Acknowledgements

We would like to thank Dr. Kevin J. Peterson, Dartmouth College, Hanover, NH, for his help with phylogenetic analyses of the miR-506 miRNA genes. This work was supported by grants from the NIH (HD071736, HD085506, HD098593, HD099924, and P30GM110767 to WY), NIH/National Center for Advancing Translational Science (NCATS) UCLA CTSI (UL1TR001881-01 to ZW and WY), and the Templeton Foundation (PID: 61174 to WY). Work in the Lai lab was supported by the NIH (R01-GM083300 and R01-HD108914) and Memorial Sloan-Kettering Institute Grant P30-CA008748.