Rapid evolutionary dynamics of an expanding family of meiotic drive factors and their hpRNA suppressors

Abstract

Meiotic drivers are a class of selfish genetic elements whose existence is frequently hidden due to concomitant suppressor systems. Accordingly, we know little of their evolutionary breadth and molecular mechanisms. Here, we trace the evolution of the Dox meiotic drive system in Drosophila simulans, which affects male-female balance (sex ratio). Dox emerged via stepwise mobilization and acquisition of multiple D. melanogaster gene segments including from protamine, which mediates compaction of sperm chromatin. Moreover, we reveal novel Dox homologs and massive amplification of Dox superfamily genes on X chromosomes of its closest sisters D. mauritiana and D. sechellia. Emergence of Dox loci is tightly associated with 359-class satellite repeats that flank de novo genomic copies. In concert, we find coordinated diversification of autosomal hairpin RNA-class siRNA loci that target subsets of Dox superfamily genes. Overall, we reveal fierce genetic arms races between meiotic drive factors and siRNA suppressors associated with recent speciation.

Extended Data Fig. 1 |. BLASTn identity matrices for components of the Dox meiotic driver network.

(A) percent identity matrix for BLAST alignment for CG8664 (segment 1) shown in main Figure 1B. An example key to interpret the table is highlighted in yellow, for example the identity of MDox to CG8664 is 60.7%. (B) percent identity matrix for BLAST alignment for Prot/GD21981 (segment 2) shown in Figure 1C. The % identity shown in B is from a global alignment comparing 2.9Kb Prot/GD21981 gene to segments derived from protamine. (B’) percent identity matrix for a local alignment where there is high homology between Prot/GD21981 and segment 2 shown in Figure 1C. (C) percent identity matrix for BLAST alignment for CG15306 (segment 4) shown in Figure 1D. (D) percent identify matrix for BLAST alignment shown for Cubilin (segment 5) in Figure 1E.

Extended Data Fig. 2 |. Nucleotide and amino acid alignments for additional putative Dox ORFs (ORF13 and ORF5).

In two instances, ORF5 has accumulated premature stop codon at Nmy hpRNA left arm, and MDox. There are also two instances of frame-preserving indels at ORF5. Similarly, there are 2 instances of frame-shifting mutations, which resulted in a premature stop at UDoxA and Dox at ORF13. For ORF5 in Dsim, using a comparable ORF upstream of ProtA/B in Dmel as an outgroup, we found an excess of both non-synonymous fixed (7) and polymorphic sites (26) compared to synonymous fixed (1) and polymorphic (8) changes from MK test (MK test, P-value = 0.66). Similarly, for ORF13, we found an excess of both non-synonymous fixed (8) and polymorphic (3) sites compared to synonymous fixed (1) and polymorphic (2) changes (MK test, P-value = 0.51). While MK test does not support directional selection, frame preserving indels indicate likely functional roles for these ORFs.

Extended Data Fig. 3 |. Phylogeny of putative ORFs, ORF5 and ORF13, and the HMG-box containing ORF found at Dox.

When Dmel is used as an outgroup, the alignable portion of ORF13 is 130nt, and the resulting tree is unresolved. However, the phylogeny of sequences bearing ORF5 (derived from non-coding portion of protamine) resolve the relationships between extant Dox family members and their hpRNA suppressors, such that Tmy is related to the PDox family while Nmy is related to the Dox/MDox genes. A similar relationship is recapitulated in the phylogeny constructed from HMG-box ORF (derived from protein-coding portion of protamine). Numbers in the tree branches show posterior probability from MrBayes analysis (outgroup is Dmel ORF13/ORF5/HMG-ORF, substiution model: GTR, rate variation: equal, chain-length: 10000, Heated chains: 4, burn-in length: 1000, random-seed: 13804, and we used the following unconstrained branch lengths as priors: Gamma (1,01,1,1). The dotted grey line is to align the branch labeis.

Extended Data Fig. 4 |. Evolution and expression of MDox.

(A) Cubilin was acquired as part of the MDox transcript following insertion at a 359 satellite located between CG15371/GD16960 and Cubilin. RNA-seq shows MDox transcript includes a portion of Cubilin, which forms the part of segment 5 seen at both MDox and Dox. The small RNAs (sRNAs) that map to MDox do not correspond to the Cubilin section.

Extended Data Fig. 5 |. Example synteny analysis with RNA-seq, small RNA, and repeat annotation tracks.

For each Dox superfamily locus, we examined the Dmel ancestral State and their synteny in the simclade. In this example, the Dox locus is shown with reference to its flanking genes CG42797 and Ptpmeg2. The Dmel ancestral State shows the flanking genes, and the corresponding repeat tracks in Dmel, which shows a block of 359 satellite repeat. In Dsim, there is an insertion of Dox in this site flanked by GD16956/CG42797 and GD16051/Ptpmeg2. RNA-seq forward and reverse are shown in black and light green respectively at both flanking genes and the Dox locus. Small RNA mapping in the forward and reverse strands are showin in dark green and red respectively at Dox. In both Dmau and Dsech, the sequence in between GD16956/CG42797 and GD16051/Ptpmeg2 contains 359 satellite repeat similar to Dmel, but the Dox locus is absent. Y-axis values represent normalized RNA-seq and small RNAcounts in their respective tracks.

Extended Data Fig. 6 |. Homology relationships between hpRNAs Nmy and Tmy and their Dox superfamily targets.

Alignment of hpRNAarms to targets reveals Nmy is highly homologous to Dox and MDox, while Tmy is highly homologous to PDox1 and PDox2. Boxed segment in black shows stretches of nucleotide homology between Dox/MDox to Nmy. Boxed nucleotides in red shows diagnostic variants, which show higher homology of Tmy to PDoxl and PDox2. These relationships suggest origins of hpRNA suppressors. Tmy likely emerged from a PDox copy, while Nmy emerged from Dox/MDox.

Extended Data Fig. 7.. RNA-seq evidence for expression of Dox superfamily amplified genes in simulans clade species.

IGV screenshots show small RNAand RNA-seq mapping to Dox superfamily loci. The flanking 359 satellite sequences are shown in the annotation track in blue.

Extended Data Fig. 8 |. Nucleotide alignments showing relationships between hpRNAs and their targets.

hpRNA/target complementarity in Dsech (A) and Dmau (B). Highlighted background shows subset of Dox superfamily loci related to a specific hpRNA. In Dsech, a subset of UDox family loci are highly similar to Tmy at 3R, while another subset of UDox family loci are similar to the novel mini-Tmy hpRNA cluster in 2L. Similarity defining polymorphisms are shown in boxes. Similarly, in Dmau, a subset of Dox family loci are similar to Nmy, while a subset of UDox family loci are homologous to Tmy. Interestingly, one locus UDox3 shows higher homology to Tmy in its 5’ end but greater homology to Nmy in its 3’ end, indicating complex relationships and evolution of these hpRNAs and their suppressors.

Extended Data Fig. 9 |. Synteny relationships of all Dox superfamily loci.

(Top) The genomic interval on chromosome X that contains Dox superfamily loci in the three simulans clade species (Dsim, Dmau, Dsech) is syntenic and preserves gene order in Dmel, which lacks Dox genes. The presence of Dox superfamily loci (subdivided into Dox/MDox, PDox and UDox families), is designated below. (Bottom) The syntenic genomic locations of each of the Dox superfamily insertions in the four mel-complex species is shown. Note that, as a rule, Dox superfamily loci are flanked by 359 satellite repeats. Moreover, all but one of these Dox superfamily insertion regions (#12) contains a 359 satellite block in the ancestral State (Dmel). Beyond Dox superfamily genes, there is also co-insertion of certain other gene families, notably Ptpmeg2 and mkg-p.

Extended Data Fig. 10 |. Phylogenetic relationships between hpRNAs and theirtargets in Dsim, Dmau, and Dsech.

Tmy and Nmy hpRNAs in Dsim target distinct Dox superfamily members PDox and Dox/MDox, respectively. In Dmau and Dsech, a non-syntenic Tmy2 targets a distinct UDox family. However, in Dmau and Dsech, the distinction of hpRNA: target relationships at the phylogenetic level are not well resolved due to lack of appropriate outgroups, and potential gene conversion events confounding 1:1 relationships. The asterisk in Dmau hpRNA: target phylogeny is one such example of a likely gene conversion event which clouds separation of Tmy targets and Nmy targets into distinct clades. All phylogenetic trees were constructed using MrBayes plugin in Geneious software. We assessed MCMC convergence of trees with chain-lengths ranging from 100000-500000, and used a chain-length from the trace file, which represented converged tree. Values in the nodes represent posterior probability. We used the following MCMC settings [Chain length: 500000 (Dsim); 150000 (Dmau and Dsech); Subsampling frequency: 100000; Burn-in length: 1000] for the trees shown above.

Figure 1.

Structure of Dox transcript with segments acquired from various genes on the path to its origin. (A) Testis RNA-seq data shows a multi-exonic transcript from the Dox region, with several distinct segments acquired from protein-coding genes and repetitive elements. 359 corresponds to sequence with similarity to the 359 (also known as 1.688 family) satellite repeat. Segment 1 (blue) corresponds to sequence acquired from GD15682 (CG8664); embedded within this segment is 82 bp derived from DNAREP1 transposable element (turquoise). Segments 2 (green) and 3 (orange) correspond to sequences acquired from GD21981 (Protamine). Segment 3 is from the protein-coding portion of Protamine, which harbors an HMG-box domain. Inset (A1’) highlight amino acid identify/similarities between Dox family genes and both the Protamine HMG-box domain, and more distant HMG-box sequences from human (pfam definition) and Sox4 (pfam00505) as an outgroup. Segment 4 (yellow) was acquired from CG15306, and segment 5 (pink) derives from Cubilin. The key depicts Dox segment features, including their segment number, nucleotide length, and origin. (B) Overlap of various genomic regions to GD15682 (CG8664) from BLAST search. Segment 1 (blue) corresponds to 887bp from C-terminus and 3’ UTR of GD15682 (CG8664). BLAST hits of various lengths to different genomic features on ChrX and Chr3R are shown as light blue bars with lengths of nucleotide homology indicated. 82 bp of Segment 1, which corresponds to DNAREP1 transposable element, retrieves 126 BLAST hits in the Dsim PacBio genome (C) Genomic matches to the ancestral protamine (GD21981) gene, include regions with similarity to its upstream/5’ UTR/intronic regions (green), and others bearing the HMG box domain (orange). (D) Segment 4 from Dox was acquired from CG15306. CG15306 is no longer extant in Dsim, but relics from the insertion can be identified from BLAST search at Dox superfamily genes and their hpRNA suppressors. (E) Segment 5 from Dox was acquired from C terminus of Cubilin (pink). This segment is found only at Dox and MDox. Note that Cubilin matches to the antisense strands of Dox and MDox.

Figure 2.

Stepwise origins of Dox from a Protamine-like ancestor. Upper right, key for the labeling of gene names and structures. Note that many regions correspond to extant genomic loci in Dmel (purple) and Dsim (yellow), but mobilizations occurred in a simulans-clade ancestor; they are not meant to indicate mobilizations occur between contemporary Dsim and Dmel. In some cases, the inferred events are no longer present in contemporary species (dotted blue boxes). (A) Dmel Protamine is tandemly duplicated (ProtA/B), while Dsim has a single copy (GD21981); Dmel is a derived state. Protamine segments acquired by Dox genes are green (segment 2) and orange (HMG box). (B) Juxtaposition of segments 1–2-HMG occurred upon insertion of Protamine between CG8664 (turquoise box in CG8664 3’ UTR corresponds to DNAREP1 fragment) and forked genes, which we term the “original-Dox” (ODox). (C) ODox is inferred as an ancestral intermediate, since the contemporary Dsim GD15682 locus lacks the HMG segment and only contains fused segments 1–2. (D) Another Dsim locus exhibits segments 1–2 without the HMG box, derived from insertion into CG5004 (forming Dsim GD17329). (E) The Dox lineage with HMG acquired segment 4 by ancestral insertion of ODox into CG15306. We refer to HMG-bearing insertion into the Dsim ancestor as ODox2, again to reflect that it is not retained in present day Dsim. (F) The contemporary Dsim genome contains an unannotated gene referred to as Ur-Dox, which lacks the HMG box. However, ancestry to “ODox2” is reflected in the fact that the syntenic regions in Dsech and Dmau contain HMG box-containing Dox superfamily genes (see also Fig. 3). (G) Inferred insertion of ODox2, which juxtaposes segments 1–2-HMG-4 into a 359 segment, in the intron of GD24701 (CG43740) yielded ParaDox (PDox1). (H) A nearly identical, dispersed copy (PDox2) is present in Dsim. (I) Mobilization of PDox between Dsim homologs of CG15317 and Cubilin generated MDox. (J) Dox was generated by mobilization of MDox between Dsim ancestors of Ptpmeg2 and CG42797. (K) Summary of mobilizations from an ancestral autosomal protamine copy through multiple regions of the X chromosome, ultimately yielding the contemporary Dsim Dox gene.

Figure 3.

Evolution and diversification of the Dox superfamily in simulans-clade species. (A) Chromosomal view of expansion of Dox superfamily and non-Dox family expansions in ~ 1Mb genomic window on the X. Blue tiles show flanking genes as genomic regions to orient expanded copies of Dox superfamily members. (B) Classification of Dox superfamily loci into three subfamilies (PDox, Dox, and UDox) was based on amino acid similarity. The highlighted window within the protein alignment indicates the conserved HMG-box domain shared between Protamine and Dox family genes. (C) Phylogenetic tree showing similarity relationships among Dox superfamily loci. Highlighted boxes in the tree show clustering of Dox superfamily members into three subfamilies (PDox, Dox, and UDox) based on sequence similarity. Numbers in the tree nodes indicate posterior probability obtained from MrBayes analysis. Dmel ProtA/B and Dsim Prot/GD21981 were used as outgroups.

Figure 4.

Examples of modes of Dox superfamily expansions in the simulans-clade. (A) Insertion of Dox superfamily loci is associated with 359 satellite repeats. Synteny analysis in the mel-complex revealed 359 sequences at insertion sites are conserved evidenced by their syntenic presence in Dmel. (B) Within simulans-clade examples of independent insertions of Dox superfamily members were found indicating their active spread within species. Novel, species-specific insertions/expansions are also associated with 359 satellite repeat at insertion sites, and the synteny of 359 repeat is preserved in other species that lack an insertion. (C) Spread of Dox superfamily loci is also linked to co-amplification of two non-Dox family genes on the X chromosome. Ptpmeg2 and mkg-p gene amplification harbors signatures similar to Dox superfamily expansion, where these non-Dox family genes preferentially inserted at 359 satellite regions. Ptpmeg2 and mkg-p co-amplification events show synteny at some instances (insertion between CG32694 and CG33557), but also harbor independent insertions similar to Dox superfamily genes. A detailed synteny analyses of expansions of Dox superfamily with amplification of Ptpmeg2 and Mkg-p are shown in Extended Data Fig. 9.

Figure 5.

hpRNA:target evolution in the Dox superfamily network. (A) Chromosome map of recently-emerged autosomal hpRNAs targeting X-linked Dox superfamily in the simulans-clade. (B) Nmy hpRNA is syntenic only in Dsim and Dmau, but not Dsech. Flanking sequences reveal a gap at the corresponding Dsech region, while Dsim and Dmau Nmy share flanking genes CG14369 and CG31337. (C) Tmy hpRNA is nearby Nmy on Chr3R and is unique to Dsim. Alignments show presence of Tmy only in Dsim, and gaps in Dmau and Dsech. However, the flanking genes CG4525 and CG5623 are preserved in all three species. (D) Tmy2 hpRNA is syntenic between Dmau and Dsech. Tmy2 emerged via insertion of a Dox superfamily member between Gr98d and Klp98A, followed by duplication to generate an hpRNA in the ancestral species. Dsim Gr89b and Klp98A exhibit Dmel-like ancestral state, but these genes are disrupted in Dmau and Dsech due to hpRNA birth at this locus. (E) The mini-Tmy-like hpRNA complex (mTmy-C) in Dsech. Gene models show Dmel ancestral state and location of the emergence of duplicated Tmy-like hpRNA cluster. The ancestor to the hpRNA cluster disrupted CG13131 and in the contemporary state, this hpRNA is flanked by Ndf and CG13127 genes. Local duplications also affect the flanking gene Trp1. Secondary structure for one hpRNA unit of the mTmy-C cluster. (F) Example of hpRNA:target relationships, with resulting small RNAs with antisense complementarity to targets. Number of small RNAs observed in w[XD1] testis dataset is shown in parentheses. For panels B, C, D, and E, normalized sRNA and RNA-seq tracks depict the structure and expression from hpRNA. Forward strand sRNA reads are shown in dark green and reverse strand sRNAs in red. Similarly, forward strand RNA-seq reads shown in black, and reverse strand RNA-seq reads in light green. Y-axis values indicate normalized read counts for sRNA and RNA-seq.

Figure 6.

Genomic expansion of protamine genes in D. melanogaster. As noted, the autosomal protamine locus is in a derived state in Dmel, as it is locally duplicated, unlike simulans-clade and outgroup Drosophila species (Figure 2A). (A) Genome browser tracks of small RNA data (yellow) and RNA-seq data (purple) from control (ctrl) and ago3 heterozygous (over TM6) testis, as well as from piRNA pathway mutant testis (ago3 and aubergine/aub). Two regions of the assembled Y chromosome (central red box) are shown as enlargements. Top, expansion of Mst77 genes. Protamine belongs to the MST-HMG box family, for which the autosomal member Mst77F was previously observed to have broadly expanded on the Y chromosome (Mst77Y cluster). Bottom, adjacent to the Mst77Y cluster, in the h17 cytoband, is another cluster bearing repeated portions of multiple protein-coding genes, including protamine. The annotated genes in the Mst77Y cluster are associated with testis RNA-seq evidence, but not small RNA data. By contrast, the h17 cluster is associated with abundant small RNA data, but not RNA-seq data. Most of these reads seem to be piRNAs, since they are depleted in piRNA mutant testes. (B) Alignment of MST-HMG box members indicates that the h17 copies on the Y are more closely related to Protamine than to Mst77 members or other MST-HMG box members. (C) Small RNAs from the h17 cluster are mostly piRNA-sized (~23–28 nt in these data) and are depleted in testis mutated for the piRNA factor aubergine (aub). The minor population of h17 cluster small RNAs remaining in aub mutant testis appear siRNA-sized (21 nt).