Teosinte Pollen Drive guides maize diversification and domestication by RNAi

Abstract

Selfish genetic elements contribute to hybrid incompatibility and bias or 'drive' their own transmission^1,2. Chromosomal drive typically functions in asymmetric female meiosis, whereas gene drive is normally post-meiotic and typically found in males. Here, using single-molecule and single-pollen genome sequencing, we describe Teosinte Pollen Drive, an instance of gene drive in hybrids between maize (Zea mays ssp. mays) and teosinte mexicana (Z. mays ssp. mexicana) that depends on RNA interference (RNAi). 22-nucleotide small RNAs from a non-coding RNA hairpin in mexicana depend on Dicer-like 2 (Dcl2) and target Teosinte Drive Responder 1 (Tdr1), which encodes a lipase required for pollen viability. Dcl2, Tdr1 and the hairpin are in tight pseudolinkage on chromosome 5, but only when transmitted through the male. Introgression of mexicana into early cultivated maize is thought to have been critical to its geographical dispersal throughout the Americas³, and a tightly linked inversion in mexicana spans a major domestication sweep in modern maize⁴. A survey of maize traditional varieties and sympatric populations of teosinte mexicana reveals correlated patterns of admixture among unlinked genes required for RNAi on at least four chromosomes that are also subject to gene drive in pollen from synthetic hybrids. Teosinte Pollen Drive probably had a major role in maize domestication and diversification, and offers an explanation for the widespread abundance of 'self' small RNAs in the germ lines of plants and animals.

Fig. 1. Single-pollen sequencing reveals selfish inheritance in TPD.

a, Anther florets (5 mm) from wild-type (WT; left) and TPD (right) plants. Scale bars, 1 mm. b, Mature pollen grains from WT (left) and TPD (right) plants. Arrowheads denote developmentally arrested pollen grains. Scale bars, 0.1 mm. c, Viable pollen grains are plump and darkly stained with iodine potassium iodide (I₂KI), whereas arrested pollen grains (arrowheads) exhibit reduced diameter and incomplete staining. Scale bars, 0.1 mm. d, Quantification of pollen abortion rates in TPD backcross (BC_11,12), WT and TPD self-fertilized (BC₈S₃) lines. Data are mean ± s.d. (n = 6–8). ****P < 0.0001 and not significant (NS; two-tailed Student’s t-test). e, Phenotypic segregation ratios in replicate reciprocal crosses. The numbers above the bar represent the sample size for each progeny population. The red dashed lines denote a perfect 2:1:1 phenotypic segregation ratio. f, Fluorescein diacetate (FDA) viability staining of tetrads from TPD plants. Pollen viability is progressively restricted to a single spore following meiosis. Panels show differential interference contrast (DIC), FDA and merged images. Scale bars, 50 µm. g, Viability scoring of TPD and WT tetrads shown in panel f. TPD spores exhibit significantly reduced viability at the tetrad stage. n = 3 biological replicates, 952 total tetrads assayed. Data are mean ± s.d. *P < 0.05 and **P < 0.01 (Welch’s t-test). h, Single-pollen grain genome sequencing. Imputed allele frequencies at mexicana markers in a population of 178 mature pollen grains collected from TPD plants. Chr. chromosome. i, Imputed mexicana marker density on chromosomes 5 and 6 for individual pollen grain genome sequences. Multiple mexicana haplotypes (blue) are selfishly inherited in viable TPD pollen grains (n = 178) but not in WT pollen grains (n = 32). Values shown are plotted using a 500-kb sliding window (h,i).

Fig. 2. A toxin–antidote system introduced from mexicana on chromosomes 5 and 6.

a, Representative tassels from fertile, semi-sterile and sterile plants in a maternally segregating population. Scale bar, 1 cm. b, I₂KI viability staining of pollen from the same genotypes as in panel a. Scale bar, 0.1 mm. c, Measurement of days to anthesis in fertile, semi-sterile and sterile phenotypic classes. n is given at the bottom of the plots. ****P < 0.0001 (two-tailed Mann–Whitney test). d, Genotypic segregation ratios in reciprocal crosses. The numbers at the top of the bars represent the sample size for each progeny population. The red dashed lines denote a perfect 1:1:1:1 genotypic segregation ratio. Normal segregation is only observed in maternal progeny. e, Bulk segregant analysis of fertile and sterile progeny pools indicates that Tpd1 (red arrowhead) is necessary and sufficient for dominant male sterility (toxin), whereas Tpd2 (blue arrowhead) is associated with fertility (antidote). FDR ≤ 0.01 (Benjamini–Hochberg method). f, Dot plots of chromosomes 5 and 6 showing multiple alignment between the TPD and W22 reference genomes. The blue lines and shaded regions correspond to five fully scaffolded intervals of mexicana introgression (indicated by arrowheads). As in panel e, the red and blue arrowheads mark the Tpd1 and Tpd2 intervals, respectively. The small purple arrowheads indicate breakpoints of an approximately 13-Mb paracentric inversion present within the Tpd1 haplotype on chromosome 5L. g, Schematics summarizing the Tpd1 and Tpd2 intervals, as well as associated markers. The 13-Mb inversion is indicated as a reverse arrow.

Fig. 3. Dcl2 from teosinte is a linked antidote for toxic 22-nt siRNA.

a, Genome-wide mexicana SNP density in bulk-sequenced Tpd1;Tpd2 (BC₈S₃) plants. A subset of mexicana introgression intervals (in addition to Tpd1 and Tpd2) are selectively maintained and include RNAi factors. A mexicana-derived allele of Dcl2 (dcl2^T) with a high rate of nonsynonymous substitution is maintained in linkage to Tpd1. dsRBD, double-stranded RNA-binding domain. b, Rates of recombination between dcl2^T and Tpd1 in replicate reciprocal crosses. dcl2^T exhibits tight pseudolinkage with Tpd1 when propagated as male (0 cM), but not as female (18.7 ± 1.6 cM). The numbers above the bars represent the sample size for each progeny population. c, Measurements of pollen viability in Tpd1/tpd1 and tpd1 plants containing combinations of Dcl2, dcl2^T and dcl2-mu1. Addition of the dcl2-mu1 hypomorphic allele is sufficient for suppression of Tpd1-mediated pollen abortion. Data points correspond to measurements from individual plants (n = 6–10). **P < 0.01 and ***P < 0.001 (two-tailed Mann–Whitney test). d, Volcano plots of 21-nt (n = 9), 22-nt (n = 212) and 24-nt (n = 6) small RNA (sRNA) clusters that are differentially expressed in WT and TPD pollen. The accumulation of ectopic 22-nt siRNAs occurs specifically in TPD pollen. log₂ fold change ≥ 2, FDR ≤ 0.01. e, Representative ears from replicate crosses containing WT Dcl2 (W22 × Tpd1/tpd1) or dcl2-mu1 (W22 × dcl2-mu1 Tpd1/dcl2-mu1 tpd1) in linkage to Tpd1. Pollen parents homozygous for dcl2-mu1 restore the seed set. Scale bar, 4 cm.

Fig. 4. 22-nt siRNAs from a mexicana-derived hairpin (Tpd1) target Tdr1, an essential pollen gene.

a, 22-nt siRNA levels at the Tdr1 locus in leaf and pollen tissue from WT and TPD genotypes. Ectopic 22-nt siRNAs accumulate in TPD pollen specifically. CPM, counts per million. b, iPARE-seq depicting the accumulation of 3′-OH cleavage products at the Tdr1 locus. Tick marks indicate predicted target sites for hp-siRNAs derived from the Tpd1 hairpin. Sites with (red) and without (grey) iPARE read support are shown. c, 22-nt hp-siRNA accumulation at the Tpd1 hairpin. The hairpin locus is disrupted by transposable element insertions in the W22 genome. Data shown are normalized CPM (panels a–c). d, 22-nt hp-siRNA abundance at the Tpd1 hairpin locus in WT and TPD pollen. n = 3 replicates per condition. ****P < 0.0001 (Mann–Whitney test). e, Average size distribution of reads mapping to the Tpd1 hairpin. f, Small RNA target site prediction at the Tdr1 locus using psRNATarget. Counts indicate unique hp-siRNAs from Tpd1 that target each cleavage site. g, Homology between the guide strand (black) and the target strand (orange) is shown for the four most abundant hp-siRNAs. The tenth (red) and eleventh nucleotides in the guide strand flank the site of AGO-mediated cleavage. Tpd1-hp-siRNAb is predicted to suppress translation. h, CRISPR–Cas9 targeting of the Tdr1 locus. Edits corresponding to tdr1-1 and tdr1-2 (blue) are shown. 1F, 1R, PCR primers; sgRNA, single guide RNA. i, Developmentally synchronized tassels from WT and tdr1-mutant T0 plants. tdr1 mutants exhibit severely delayed anthesis. Scale bars, 3 cm. j, Mature 5-mm anthers from WT and tdr1-mutant T0 plants. Scale bars, 1 mm. k, I₂KI viability staining of pollen from WT and tdr1-mutant T0 plants. Scale bars, 0.1 mm.

Fig. 5. Tpd1 hp-siRNA target site deletion in tdr1 has spread to modern maize from teosinte.

a, Sequence complement of the Tpd1-hp-siRNAa and Tpd1-hp-siRNAb target sites in Tdr1, indicating a 27-bp in-frame deletion found in modern maize, maize traditional varieties and in teosinte that removes the Tpd1-hp-siRNAa seed sequence, and a SNP on the eleventh nucleotide of Tpd1-hp-siRNAb that is predicted to reduce binding. b, Pie charts indicating the frequency of the deletion in 1,483 resequenced genomes from maize and teosinte, aligned with the B73 reference genome (GATK 3.0). The deletion allele (blue) arose in teosinte and quickly spread through maize traditional varieties in Central and South America, before fixation in modern stiff stalk, but not in tropical maize inbred lines. High frequencies of heterozygosity in mexicana and parviglumis are consistent with recent or ongoing pollen drive.

Extended Data Fig. 1. Teosinte Pollen Drive and genetic mapping of Tpd1 and Tpd2.

a, Crossing scheme of the TPD phenotype. When back-crossed as male, all the progeny of semi-sterile TPD plants display the semi-sterile pollen phenotype instead of the expected 1:1 fertile:semi-sterile ratio. Graphics in panel a were created using BioRender (https://biorender.com). b, Representative ears from Tpd1 Bt1/bt1; Tpd2 Y1/y1 reciprocal crosses with bt1; y1 testers, demonstrating severe segregation distortion (“drive” of Bt1), but only through the male. bt1 (brittle1, collapsed kernels); y1 (yellow1, white kernels). c, Summary of molecular and morphological mapping of the Tpd1 interval. Molecular mapping was performed using Tpd/++ x W22 segregating progeny, whereas morphological mapping was performed by crossing Tpd1 Bt1/tpd1 bt1 plants to bt1 testers. d, Molecular mapping of the Tpd2 interval. SNP markers are shown in blue with recombination frequencies in red.

Extended Data Fig. 2. mexicana intervals introgressed into maize carry RNAi genes.

a, Whole genome plots of homozygous mexicana SNP density present within Tpd1; Tpd2 lines. The upper plot corresponds to data from bulked seedlings after 8 backcrosses and 3 self-pollinations (BC₈S₃) whereas the lower plot is from BC₅S₂ plants. SNP density is consolidated in 250 kb genomic bins. Physical locations for morphological markers Bt1 and Y1, as well as mexicana derived RNAi genes, are labelled in red. 7/13 introgression intervals overlap in both independently maintained homozygous lines. b, Allele frequency at mexicana markers in 96 pollen grains from four different TPD plants subjected to single pollen grain sequencing. Regions highlighted in red were over-represented in viable pollen grains. c, Quantification of pollen viability in Tpd1 + /−; Tpd2 + /−; rgd1/+ and Tpd1 + /−; Tpd2 + /−; Rgd1 pollen demonstrating gametophytic suppression via germline segregation of the rgd1 null allele. n ≥ 9 plants per genotype, ≥ 200 pollen grains per plant. **** p < 0.0001 (Welch’s t-test).

Extended Data Fig. 3. DCL2-dependent 22nt siRNAs from hairpins are prevalent in maize pollen.

a, Distribution and relative abundances of small RNA size classes in WT pollen libraries. Bars indicate mean ± SD. n = 3 biological replicates. b, Comparison of relative abundances for 21-nt, 22-nt, and 24-nt sRNA size classes in WT leaf and pollen samples. Both 22-nt and 24-nt sRNAs show significant increases in pollen. Bars indicate mean ± SD. n = 3 biological replicates. **** p < 0.0001 (Welch’s t-test). c, Comparison of 22-nt sRNA levels in Dcl2, dcl2^T, and dcl2-mu1 pollen at 804 pollen-specific loci. Values shown are log₂ transformed counts per million (CPM) averaged across replicates. n = 3 replicates per genotype. **** p < 0.0001 (ANOVA test). d, Summary of relative contributions for 22-nt sRNA producing loci in WT pollen. Hairpin/inverted repeat (IR) hp-siRNAs represent the largest fraction of 22-nt species. e, Heatmap showing 22-nt hp-siRNA levels at hpRNA loci in leaf and pollen. f, Heatmap showing 22-nt siRNA levels at protein-coding genes in leaf and pollen. g, Browser shots showing 22-nt hp-siRNA accumulation at a hpRNA locus on chromosome 1 (left) and 22 nt siRNA silencing at a representative protein-coding gene. Scale is CPM.

Extended Data Fig. 4. Validation of highly abundant pollen hairpin precursors.

a, hp-siRNAs are expected to show strand bias. Measurement of strand score (min[plus, minus]/max[plus, minus]) at 28 putative hairpin precursors and randomly selected siRNA clusters from wild-type (W22) maize. A value of zero indicates complete strand bias, whereas a value of 1 indicates unbiased accumulation from both strands. n = 28. **** p < 0.0001 (Welch’s t-test). b, log2 read count at hairpin precursors indicate 22 nt size bias. n = 28. **** p < 0.0001 (ANOVA). c, Example of a 73 nt stretch from a 4,480 bp hairpin precursor demonstrating near-complete reverse complementarity. d, Mountain plots measuring thermodynamic stability of the Tpd1 hairpin from mexicana and another randomly selected hairpin structure.

Extended Data Fig. 5. Origins and targets of 22 nt small RNAs in TPD Pollen.

a, RNAi genes (Sgs3/Rgd1, Rdr6, Ago1e and Dcl2) associated with 22 nt biogenesis and function are upregulated in TPD pollen. Expression is shown in TMM normalized counts. Bars show mean ± SD. n = 5 replicates per condition. **** p < 0.0001, *** p < 0.001, ** p < 0.01 (FDR). b, Relative abundances of TPD-dependent 22 nt siRNAs mapping to annotated elements. Pie chart inset shows proportions of 22 nt siRNAs targeting genes in CPM. c, Browser shot showing transcriptional activation at PIF/Harbinger elements in TPD pollen as well as 22 nt siRNA accumulation. d, Quantification of mRNA expression at 258 PIF/Harbinger superfamily elements in WT and TPD pollen. **** p < 0.0001 (Mann-Whitney test). e, 22 nt siRNA levels at 42 PIF/Harbinger elements in WT and TPD pollen. **** p < 0.0001 (Mann-Whitney test).

Extended Data Fig. 6. TPD-dependent silencing of a GDSL lipase disrupts lipid metabolism.

a, Browser shots showing ectopic accumulation of 22-nt siRNAs at protein-coding genes in TPD pollen. Scale in counts per million (CPM). b, RNA-seq tissue expression of 22-nt siRNA targets specific to TPD pollen, data from ref. . Bars show mean ± SD. n = 3 replicates per tissue. c, RNA-seq expression of 22-nt siRNA targets in WT and TPD pollen. Bars show mean ± SD. n = 5 replicates per condition. d, Western blot comparing TDR1 protein levels in WT and TPD pollen, anthers, and leaf. Protein levels were normalized using Heat Shock Protein 90 (HSP90). e, p-nitrophenyl butyrate esterase activity assay in 5 mm anthers and pollen from WT, TPD, and Tpd1 + /− plants. f, g, GO term biological processes up-regulated in f, TPD and g, WT pollen (FDR ≤ 0.001). Upregulated genes in TPD pollen were associated with RNA metabolism, ribosome assembly, and cytoplasmic translation as well as G2 mitotic arrest. This could reflect translational repression via 22-nt siRNAs. Interestingly, a subset of genes associated with endoplasmic reticulum (ER)-nucleus signalling was also up-regulated, while genes associated with glycerol metabolism, the primary backbone for TAG synthesis, were downregulated. In pollen, the accumulation of TAGs in lipid droplets (LDs) is critical for proper membrane expansion and pollen tube growth.

Extended Data Fig. 7. Tpd1 and Dcl2 are expressed pre-meiotically, whereas Tdr1 is expressed in microspore and pollen.

a, RT-qPCR of the Tdr1 transcript throughout anther development and in mature pollen. Bars show mean ± SD. n = 2 replicates per condition. b, RT-qPCR of the Tpd1 transcript during anther development and in mature pollen in WT and Tpd. Bars show mean ± SD. n = 2 replicates per condition. c, d, Single-cell expression at different stages of meiosis of c, Tdr1 and d, Dcl2, using single cell RNAseq data. Early and late expression of Dcl2 coincides with Tpd1 and Tdr1, respectively.

Extended Data Fig. 8. Tpd2 suppresses 22nt secondary small RNAs.

a, Browser shots showing ectopic accumulation of 22nt siRNAs at Tdr1 (left) and Tpd1 hairpin (right) in fertile (tpd1; tpd2, grey), drive (Dcl2^T Tpd1/Dcl2 tpd1; Tpd2/tpd2, blue) and sterile (Dcl2 Tpd1/Dcl2 tpd1; tpd2, red) pollen from maternal segregants. Scale in counts per million (CPM). Tpd2 and Dcl2^T reduce small RNAs from Tdr1 (left) but not from the Tpd1 hairpin (right), consistent with a cell autonomous role in secondary small RNA biogenesis and silencing. b, Table summarizing genotypes transmitted when TPD is backcrossed to W22 as female (left column) or male (right column). Only the combination of dcl2^T Tpd1 (linked) and Tpd2 is transmitted through pollen. Recombinants between dcl2^T and Tpd1 occur at the expected frequency but are not transmitted through pollen (Fig. 3b), presumably because of higher siRNA production during and after meiosis. c, Rdm1 (Zm00004b029511) is one of six genes in the Tpd2 interval expressed in pollen, and is overexpressed in TPD pollen.

Extended Data Fig. 9. Signatures of Teosinte Pollen Drive in modern maize, maize traditional varieties and sympatric mexicana populations.

a, Frequency of mexicana-derived alleles were calculated for 1 Mb intervals associated with TPD on chromosomes 1,2,3,4,5,6 and 10. Correlations are shown between population means from each of 14 maize traditional varieties (left) and sympatric mexicana populations (right). Intervals on chromosomes 5, 6 and 10 include Dcl2 (5.19), Tdr1 (5.40), Tpd1 (5.79), Rgd1/Sgs3 (6.3) and candidate genes Ago1a (6.3), Tpd2/Rdm1 (6.98), Ago1b and Ago2b (10.134). Correlations were observed for most of the intervals in maize traditional varieties, except for Tdr1 (green arrow), but only for intervals including Tpd1, Rgd1 and Tpd2 in mexicana. Spearman correlation coefficients are displayed as a heatmap. b, mexicana-derived ancestry in each of 14 maize traditional varieties (above) and sympatric mexicana populations (below) in Dcl2, Tdr1, Tpd1 and Tpd2 intervals. The Tdr1 interval (green) is monomorphic in most of the maize traditional varieties, but shows extreme dimorphism in 7 out of 14 sympatric mexicana populations.

Extended Data Fig. 10. Mechanistic model of Teosinte Pollen Drive.

a, The TPD system is defined by mexicana introgression intervals on chromosomes 5 and 6. Tpd1 encodes a pre-meiotically expressed mexicana-specific hairpin that produces abundant 22nt hp-siRNAs. b, These hp-siRNAs trigger secondary siRNAs amplification by RDR6 and SGS3/RGD1 at the Tdr1 gene when it starts being transcribed at the late tetrad stage, which in turn target Tdr1 for translational repression (red ribosomes). In surviving microspores (dark yellow background) Tpd2 and dcl2^T repress secondary siRNAs processing, restoring translation and fertility (green ribosome). c, Only pollen grains of the genotype dcl2^T Tpd1; Tpd2 are viable, and all other competing gametes are eliminated. Other RNAi genes (Sgs3/Rgd1, Ago1, Ago2, Ago5) can act as partial suppressors by affecting levels of siRNAs.

Extended Data Fig. 11. Evolutionary model of Teosinte Pollen Drive.

After the antidotes arise in an ancestral teosinte population, the Tpd1 toxin arises and gains a transmission advantage when linked to the antidote genes. In extant populations of Z. mexicana and Z. mays, some antidotes are fixed, while others are polymorphic or lost. The demographic model was based on ref.  and the conceptual framework of selfish evolution was adapted from ref. . Graphics were created using BioRender (https://biorender.com).