Introduction {#s1}
============

Mendel's law of equal segregation stipulates that paternal and maternal alleles of a gene should have an equal chance of being transmitted to progenies. This law guarantees a fair competition between different alleles and enables beneficial ones to prevail during natural selection. Meiotic drivers, a type of selfish genetic element, break Mendel's law by skewing transmission ratios to their advantage, and thus can spread in a population even when having a deleterious effect on organismal fitness ([@bib33]; [@bib66]). The term 'meiotic drive' was initially coined to describe segregation distortion resulting from preferential inclusion in the gamate during asymmetric female meiosis but has now been used more broadly to include biased transmission caused by postmeiotic mechanisms. In fact, some of the best-known meiotic drivers, such as *Segregation Distorter* in *Drosophila* and the *t* haplotype in mouse ([@bib34]; [@bib30]), act postmeiotically to disable male gametes (sperms) that do not inherit them. This type of meiotic driver, called gamete killer, exists in animals, plants, and fungi ([@bib7]). Fungal gamete killers, or spore killers, have been found in several filamentous ascomycetes, most notably *Neurospora* and *Podospora anserina* ([@bib11]; [@bib16]; [@bib18]; [@bib62]). It is unclear how widespread spore killers are among fungal species.

The fission yeast *Schizosaccharomyces pombe* is a prominent model organism for molecular and cell biology and has been increasingly used to study natural variation and genome evolution ([@bib6]; [@bib22]; [@bib24]; [@bib53]). *S. pombe* natural isolates, which are nearly all haploids and have pair-wise nucleotide differences of less than 1% ([@bib24]; [@bib53]), can readily mate with each other to form hybrid diploids, but the viability of spores derived from inter-isolate crosses is often below 5% and in many instances under 1% ([@bib17]; [@bib27]; [@bib42]). This is in stark contrast to *Saccharomyces cerevisiae*, whose natural isolates have a similar level of nucleotide diversity but much better spore viability when inter-crossed ([@bib21]). One explanation for the within-species hybrid sterility in *S. pombe* is chromosomal rearrangement ([@bib2]; [@bib6]; [@bib69]). However, because one rearrangement reduces spore viability at most by half ([@bib2]; [@bib21]), other factors are likely in play. It was shown recently that when the *S. pombe* laboratory strain, which was isolated from French grape juice ([@bib22]; [@bib49]), was crossed to a strain isolated from fermented tea (initially called *Schizosaccharomyces kambucha* \[[@bib57]\], later renamed *S. pombe var. kambucha* \[[@bib53]\]), at least three spore killers contributed to hybrid sterility ([@bib4]; [@bib69]), but the identities of the killer genes were unknown.

In this study, through investigating the causes of intraspecific hybrid sterility of fission yeast, we uncovered the molecular identities of two active spore killers, which both belong to the *wtf* gene family. We analyzed their killing behaviors in both native and non-native genomic contexts. Furthermore, we performed comparative genomic analysis of *wtf* genes and revealed interesting patterns of divergence among them.

Results {#s2}
=======

We have previously studied CBS5557, an *S. pombe* strain isolated from grapes in Spain ([@bib22]; [@bib52]). Like most *S. pombe* isolates, CBS5557 differs from the laboratory strain, whose genome is the *S. pombe* reference genome, in a 2.23 Mb pericentric inversion on chromosome I ([@bib6]). After introducing the same inversion into the laboratory strain ([@bib22]), the viability of spores derived from the cross between the laboratory strain and CBS5557 increased from 14% to 22% ([Figure 1---figure supplement 1](#fig1s1){ref-type="fig"} and [Supplementary file 1](#SD4-data){ref-type="supplementary-material"}). To uncover the causes of the remaining sterility, we performed next-generation sequencing (NGS)-assisted bulk segregant analysis ([@bib22]) and found that among the viable progenies from the cross, SNP alleles of the laboratory strain (reference alleles) on chromosome III were markedly under-represented ([Figure 1A](#fig1){ref-type="fig"}). Two broad transmission distortion peaks, one on the left arm (hereafter referred to as the left peak) and one on the right arm (hereafter referred to as the right peak), were observed, suggesting that there are at least two spore killer genes on chromosome III of CBS5557.10.7554/eLife.26057.003Figure 1.The identification of *cw9* and *cw27* genes as spore killers.(**A**) NGS-assisted bulk segregant analysis of viable progenies from a cross between DY9974, a heterothallic derivative of CBS5557, and DY8531, a laboratory strain harboring the 2.23 Mb chromosome I inversion ([@bib22]). Reference allele frequencies at SNP positions are plotted as dots. Red trend lines are based on a rolling median calculation. Two vertical dashed lines denote the positions of *cw9* and *cw27* genes on chromosome III. Black triangles mark the positions of centromeres. (**B and C**) Bulk segregant analysis of viable progenies from crosses between a laboratory strain and the backcrossed-1 strain, and between a laboratory strain and the backcrossed-2 strain, respectively. Blue and orange colored segments in the diagrams and on the X-axes of the plots denote chromosome III regions with reference genome sequence and those with CBS5557 genome sequence, respectively, in the two backcrossed strains. (**D and E**) Schematics of the regions surrounding *cw9* and *cw27* in the CBS5557 genome and the corresponding regions in the reference genome. Gene structures of *cw9* and *cw27* were predicted using the AUGUSTUS web server ([@bib58]). Solo LTRs were annotated based on BLAST analysis.**DOI:** [http://dx.doi.org/10.7554/eLife.26057.003](10.7554/eLife.26057.003)10.7554/eLife.26057.004Figure 1---figure supplement 1.Hybrid sterility between CBS5557 and the laboratory strain can be partially rescued by eliminating the difference in the 2.23 Mb chromosome I inversion.(**A**) Representative tetrads from a cross between LD775, an *h^-^* laboratory strain, and DY9974, an *h^+^* derivative of CBS5557. (**B**) Representative tetrads from a cross between DY8531, an *h^-^* laboratory strain harboring the 2.23 Mb chromosome I inversion ([@bib22]), and DY9974. (**C**) Quantitation of spore viability assessed by tetrad analysis. p-value was calculated using Fisher's exact test. Numerical data are provided in [Supplementary file 1](#SD4-data){ref-type="supplementary-material"}.**DOI:** [http://dx.doi.org/10.7554/eLife.26057.004](10.7554/eLife.26057.004)10.7554/eLife.26057.005Figure 1---figure supplement 2.Alignment of the 5' portions of the *cw9* gene of CBS5557 and the *wtf9* gene of the reference genome.The sequences were aligned using MAFFT via Jalview ([@bib25]; [@bib65]). Identical bases are indicated by blue background. Red and green boxes denote the first and second predicted protein coding exons, respectively. Predicted start codons are highlighted with red letters. Gene structure of *cw9* was predicted using the AUGUSTUS web server ([@bib58]).**DOI:** [http://dx.doi.org/10.7554/eLife.26057.005](10.7554/eLife.26057.005)10.7554/eLife.26057.006Figure 1---figure supplement 3.Alignment of the two directly oriented LTRs flanking *cw27* in the CBS5557 genome (top and bottom sequences, respectively) and the LTR at the corresponding location in the reference genome (middle sequence).The reference genome LTR was referred to as PCRW_055 or C3_1895301 in [@bib5]. According to the classification proposed in Bowen et al., these three LTR sequences belong to the ζ clade. The sequences were aligned using MAFFT via Jalview. Identical bases are indicated by grey background. 7-nt-long homology at the breakpoint is highlighted in red. Matched sequences upstream and downstream of the breakpoint are colored in green.**DOI:** [http://dx.doi.org/10.7554/eLife.26057.006](10.7554/eLife.26057.006)

To narrow down the genomic regions containing the spore killer genes, we performed multiple rounds of backcross and obtained strains with only a portion of chromosome III originating from CBS5557 and the rest of the genome coming from the laboratory strain (see \'Materials and methods\'). We selected two backcrossed strains that retained the ability to skew inheritance of markers inserted in the left and right peak regions, respectively. These two strains, referred to as backcrossed-1 and backcrossed-2, respectively, were crossed to the laboratory strain and viable progenies were subjected to NGS-assisted bulk segregant analysis, which revealed the locations of the CBS5557 genomic segments present in the two backcrossed strains and confirmed the allele transmission bias (left-side plots in [Figure 1B and C](#fig1){ref-type="fig"}).

Through inspecting the Illumina and PacBio sequencing data of the CBS5557 genome, we discovered that, within the CBS5557 genomic segments in the two backcrossed strains, sequences most divergent from the reference genome are those containing *wtf* family genes ([Figure 1D and E](#fig1){ref-type="fig"}). *wtf* (for [w]{.ul}ith [Tf]{.ul} LTRs) genes were so named because they are often in close proximity to solo long terminal repeats (LTRs) derived from Tf retrotransposons ([@bib5]; [@bib67]). *wtf* genes are predicted to contain multiple introns and encode multi-transmembrane proteins and have been shown to be highly up-regulated transcriptionally during meiosis and sporulation by large-scale analyses ([@bib5]). No homologs of *wtf* genes have been found in any other species, including other *Schizosaccharomyces* species ([@bib53]). We deleted from the backcrossed-1 strain a *wtf* gene locating at the position of *wtf9* in the reference genome but exhibiting a low sequence identity (42%) to *wtf9* in the 5' portion ([Figure 1D](#fig1){ref-type="fig"} and [Figure 1---figure supplement 2](#fig1s2){ref-type="fig"}) and found that the deletion notably reduced the segregation distortion (right-side plot in [Figure 1B](#fig1){ref-type="fig"}), suggesting that this *wtf* gene, which we named *cw9* (the nomenclature of *wtf* genes in the CBS5557 genome will be explained later in this paper), is a spore killer. The residual allele frequency bias is probably due to additional spore killer gene(s) in the CBS5557 genomic segments present in the backcrossed-1 strain. Within the CBS5557 genomic segment in the backcrossed-2 strain, we found a *wtf* gene absent in the reference genome, possibly owing to a deletion event mediated by homologous recombination between LTRs ([Figure 1E](#fig1){ref-type="fig"} and [Figure 1---figure supplement 3](#fig1s3){ref-type="fig"}). Removing this *wtf* gene, which we named *cw27*, from the backcrossed-2 strain completely abolished segregation distortion (right-side plot in [Figure 1C](#fig1){ref-type="fig"}).

To determine how *cw9* and *cw27* behave in their native strain background, we generated *h^+^*/*h^-^* CBS5557 diploid strains containing homozygous or heterozygous deletion of *cw9* or *cw27*, and performed tetrad analysis to evaluate spore viability ([Figure 2A](#fig2){ref-type="fig"} and [Supplementary file 1](#SD4-data){ref-type="supplementary-material"}). Homozygous deletion of *cw9* or *cw27* had no effect on spore viability ([Figure 2A](#fig2){ref-type="fig"} and [Figure 2---figure supplement 1](#fig2s1){ref-type="fig"}). In contrast, heterozygous deletion of either *cw9* or *cw27* caused a significant reduction of spore viability, which resulted mainly from the death of spores carrying the deletion ([Figure 2B](#fig2){ref-type="fig"}). When *cw9* and *cw27* were both in the heterozygous state, a more severe loss of spore viability was observed, with the double deletion mutant spores suffering a particularly high level of death ([Figure 2A and B](#fig2){ref-type="fig"}, and [Supplementary file 1](#SD4-data){ref-type="supplementary-material"}). These data indicate that *cw9* and *cw27* can act as spore killers in the self-cross of CBS5557 when either of them is in the heterozygous state, and inheriting one killer can prevent its killing effect but not that of the other killer.10.7554/eLife.26057.007Figure 2.*cw9* and *cw27* act as spore killers when in the heterozygous state in the CBS5557 background.(**A**) Heterozygous deletion but not homozygous deletion of *cw9* or *cw27* in CBS5557 *h^+^/h^-^* diploid caused spore viability loss. Spore viability was measured using tetrad analysis. Representative tetrads are shown in [Figure 2---figure supplement 1](#fig2s1){ref-type="fig"} and in panel B. p-values were calculated using Fisher's exact test. Numerical data are provided in [Supplementary file 1](#SD4-data){ref-type="supplementary-material"}. (**B**) Among the spores derived from CBS5557 diploids with heterozygous deletion, loss of viability mainly occurred to spores with deletion. Asterisks indicate significant deviation from 50% (p=1.94E-7 and 1.42E-11 for *cw9∆* spores from *cw9∆*/*cw9*^+^ diploid and *cw27∆* spores from *cw27∆*/*cw27*^+^ diploid, respectively, exact binomial test). Numerical data are provided in [Supplementary file 1](#SD4-data){ref-type="supplementary-material"}.**DOI:** [http://dx.doi.org/10.7554/eLife.26057.007](10.7554/eLife.26057.007)10.7554/eLife.26057.008Figure 2---figure supplement 1.Representative tetrads from a wild type CBS5557 *h^+^/h^-^* diploid strain (DY21782), a CBS5557 *h^+^/h^-^* diploid strain with homozygous *cw9* deletion (DY21838), and a CBS5557 *h^+^/h^-^* diploid strain with homozygous *cw27* deletion (DY21842).**DOI:** [http://dx.doi.org/10.7554/eLife.26057.008](10.7554/eLife.26057.008)

To determine whether *cw9* and *cw27* can cause spore killing in a non-native genomic context, we constructed integrating plasmids carrying *cw9* or *cw27*, integrated the plasmids at various genomic locations in the laboratory strain, and subjected the resulting strains to crosses. In crosses where only one of the parental strains harbored *cw9* or *cw27*, spores inheriting the killer had normal viability, whereas spores not inheriting the killer suffered significant loss of viability, regardless of the genomic location of the killer ([Figure 3A](#fig3){ref-type="fig"} and [Figure 3---figure supplement 1](#fig3s1){ref-type="fig"}). We then performed crosses in which both parental strains harbored a killer at the *leu1* locus ([Figure 3B](#fig3){ref-type="fig"} and [Figure 3---figure supplement 2](#fig3s2){ref-type="fig"}). When both parental strains contained the same killer, spore viability was normal. In contrast, when one parental strain had *cw9* and the other parental strain had *cw27*, spores inheriting either *cw9* or *cw27* experienced a dramatic loss of viability. These results demonstrate that *cw9* and *cw27* can act independently of genomic context and confirm that they do not exhibit mutual resistance.10.7554/eLife.26057.009Figure 3.*cw9* and *cw27* act as spore killers when inserted into the genome of the laboratory strain.(**A**) In crosses of laboratory strains, when only one parental strain had an insertion of *cw9* or *cw27* at *his3*, *lys1*, or *leu1* locus, spores without the insertion suffered viability loss. Representative tetrads are shown in [Figure 3---figure supplement 1](#fig3s1){ref-type="fig"}. p-values were calculated using Fisher's exact test. Numerical data are provided in [Supplementary file 1](#SD4-data){ref-type="supplementary-material"}. (**B**) In crosses of laboratory strains, when both parental strains had an insertion of *cw9* or *cw27* at the *leu1* locus, spore viability was normal when parents had the same killer, but was severely low when parents had different killers. The two parental alleles were distinguished by *leu1*-linked antibiotic resistance markers. Representative tetrads are shown in [Figure 3---figure supplement 2](#fig3s2){ref-type="fig"}. Numerical data are provided in [Supplementary file 1](#SD4-data){ref-type="supplementary-material"}. (**C**) Electron microscopy analysis of laboratory-background *h^+^/h^-^* diploid cells undergoing synchronous meiosis and sporulation.**DOI:** [http://dx.doi.org/10.7554/eLife.26057.009](10.7554/eLife.26057.009)10.7554/eLife.26057.010Figure 3---figure supplement 1.Representative tetrads from laboratory-background crosses in which only one of the parental haploid strains had a vector or a killer-containing plasmid integrated at the *his3*, *lys1*, or *leu1* locus.**DOI:** [http://dx.doi.org/10.7554/eLife.26057.010](10.7554/eLife.26057.010)10.7554/eLife.26057.011Figure 3---figure supplement 2.Representative tetrads from three laboratory-background *h^+^/h^-^* diploid strains homozygous for plasmid integration at the *leu1* locus and a laboratory-background *h^+^/h^-^* diploid strain heterozygous for plasmid integration at the *leu1* locus, with one allele containing *cw9* and the other allele containing *cw27.kanMX* and *hphMX* markers were inserted between the coordinates 1963090 and 1963096 on chromosome II so that they are tightly linked to *leu1*.**DOI:** [http://dx.doi.org/10.7554/eLife.26057.011](10.7554/eLife.26057.011)

To observe the ultrastructural details of the spore killing process, we performed electron microscopy analysis on laboratory-background *h^+^*/*h^-^* diploid cells undergoing synchronous meiosis and sporulation induced by a shift to a nitrogen-free synthetic sporulation medium (SSL-N medium) ([Figure 3C](#fig3){ref-type="fig"}). Two diploid strains were used: a strain heterozygous for a *cw27*-containing plasmid integrated at the *his3* locus and a control strain heterozygous for an empty vector integrated at the same locus. 12 hr after the shift, *cw27*-containing cells formed four-spored asci indistinguishable from those formed by the control cells. At this time point, organelles inside spores were clearly visible. 18 hr after the shift, the cytoplasm of spores derived from the control cells became electron-dense and organelles were no longer discernable, owing to a little-understood spore maturation process ([@bib68]). Interestingly, at the same time point, asci derived from the *cw27*-containing cells often exhibited a 2:2 pattern of spore morphology: two spores in an ascus looked like the matured spores in the control sample, whereas the other two spores exhibited cytoplasm with a lower electron density and visible organelles, an electron-lucent spore wall, and frequently a non-spherical shape. These data suggest that killer-affected spores fail to undergo a proper maturation process.

*cw9* and *cw27* are both flanked by solo LTRs, and upstream of their predicted start codons they share a nearly identical 288-bp-long sequence (hereafter referred to as the conserved_up sequence) ([Figure 4---figure supplement 1](#fig4s1){ref-type="fig"}), which is also conserved upstream of the coding sequences of many other *wtf* genes ([@bib5]). To determine whether these and other sequence features are important for spore killing, we performed truncation analysis ([Figure 4A](#fig4){ref-type="fig"} and [Figure 4---figure supplement 2](#fig4s2){ref-type="fig"}). The truncated versions were named Ta, Tb, Tc, Td, and Te. Their effects on spore viability in a heterozygous cross of laboratory strains were examined. The Ta versions, which lack the LTRs, caused spore killing phenotypes similar to those caused by longer constructs containing the flanking LTRs, indicating that the LTRs are not important. Interestingly, the Tb versions missing the sequence between the upstream LTR and the conserved_up sequence caused spores containing the killer to suffer a moderate viability loss, suggesting that these truncated killer genes partially lost the ability to protect killer-containing spores from being harmed. Removing from the Ta versions, the sequence downstream of the predicted stop codon (Tc versions) did not have a notable effect. Removing from the Tc versions, the last 10 amino-acid-coding codons (Td versions) greatly diminished spore killing, whereas removing from the Tc versions the sequences upstream of the start codons (Te versions) resulted in more indiscriminate killing than that caused by the Tb versions. These results suggest that intact C termini are important for the killing activity, and the upstream sequences are important for the protecting activity.10.7554/eLife.26057.012Figure 4.Sequence requirement for the killing and the protecting activities of *cw9* and *cw27*.(**A**) Truncation analysis to assess the involvement of the 5' and 3' sequences of *cw9* and *cw27* in spore killing. In the diagrams on top, blue arrows represent LTRs and yellow bars represent the conserved_up sequence. Representative tetrads are shown in [Figure 4---figure supplement 1](#fig4s1){ref-type="fig"}. p-values were calculated using Fisher's exact test. Numerical data are provided in [Supplementary file 1](#SD4-data){ref-type="supplementary-material"}. (**B**) The Td versions of *cw9* and *cw27* are able to effectively protect against killing despite their lack of killing activity. Representative tetrads are shown in [Figure 4---figure supplement 2](#fig4s2){ref-type="fig"}. p-values were calculated using Fisher's exact test. Numerical data are provided in [Supplementary file 1](#SD4-data){ref-type="supplementary-material"}.**DOI:** [http://dx.doi.org/10.7554/eLife.26057.012](10.7554/eLife.26057.012)10.7554/eLife.26057.013Figure 4---figure supplement 1.Alignment of the 5' portions of *cw9* and *cw27*.The starting positions of these two sequences correspond to those of the Ta truncations. The sequences were aligned using MAFFT via Jalview. Identical bases are indicated by blue background. Yellow and red boxes denote the conserved_up sequence and the first predicted protein coding exon, respectively. Predicted start codons are highlighted with red letters.**DOI:** [http://dx.doi.org/10.7554/eLife.26057.013](10.7554/eLife.26057.013)10.7554/eLife.26057.014Figure 4---figure supplement 2.Representative tetrads from laboratory-background *h^+^/h^-^* diploid strains heterozygous for plasmid integration at the *leu1* locus.**DOI:** [http://dx.doi.org/10.7554/eLife.26057.014](10.7554/eLife.26057.014)10.7554/eLife.26057.015Figure 4---figure supplement 3.Representative tetrads from laboratory-background *h^+^/h^-^* diploid strains heterozygous for a Tc-version killer-containing plasmid integrated at the *leu1* locus on chromosome II, and also heterozygous for a vector or a Td-version killer-containing plasmid integrated at the *ars1* replication origin region upstream of the *hus5* gene on chromosome I.**DOI:** [http://dx.doi.org/10.7554/eLife.26057.015](10.7554/eLife.26057.015)

To determine whether the Td versions, which are largely devoid of the killing activity, can still protect, we performed crosses in which a strain with a Td-version-containing plasmid integrated at the *ars1* locus was crossed to a strain with a Tc-version-containing plasmid integrated at the *leu1* locus, and found that the Td versions of both *cw9* and *cw27* can effectively protect spores harboring them from being killed by the respective Tc-version killers ([Figure 4B](#fig4){ref-type="fig"} and [Figure 4---figure supplement 3](#fig4s3){ref-type="fig"}). Thus, perturbing the C termini resulted in separation-of-function mutants that retain the protecting activity but can no longer kill. Such kind of mutants, if present in natural populations, would behave similarly to the resistant alleles of other previously studied gamete killers, and may also represent an intermediate state toward the extinction of a killer ([@bib7]).

In the reference genome, there are 25 *wtf* genes at 20 locations (15 singletons and five tandem pairs) ([@bib5]) ([Figure 5A and B](#fig5){ref-type="fig"}). Using PacBio sequencing, we found that in the CBS5557 genome, these 20 locations are occupied by 29 *wtf* genes (12 singletons, seven pairs, and one triplet), with four locations containing an extra *wtf* gene compared to the reference genome ([Figure 5A and B](#fig5){ref-type="fig"}). We systematically named these 29 genes according to their syntenic relationship with the 25 *wtf* genes in the reference genome, with each name consisting of the prefix *cw*, a number (from 1 to 25), and a suffix for genes at locations where extra genes exist in the CBS5557 genome ([Figure 5A and B](#fig5){ref-type="fig"}). In addition, there are three CBS5557 singleton *wtf* genes without syntenic counterparts in the reference genome and we named them *cw26*, *cw27*, and *cw28* according to their order in the genome ([Figure 5A](#fig5){ref-type="fig"}). Like *cw27*, *cw26* and *cw28* also appeared to have been lost in the laboratory strain through LTR-mediated recombination ([Figure 5---figure supplement 1](#fig5s1){ref-type="fig"}), suggesting that this is a common mode of *wtf* gene turnover.10.7554/eLife.26057.016Figure 5.*wtf* genes vary both in numbers and sequences between the reference genome and the CBS5557 genome.(**A**) Genomic locations of *wtf* genes in the reference genome and in the CBS5557 genome. The 25 *wtf* genes in the reference genome have been named according to their order in the genome ([@bib5]). Their locations are depicted as 20 blue vertical bars, including five thick bars denoting five tandem pairs. Our PacBio sequencing analysis revealed that in the CBS5557 genome, there are 32 *wtf* genes, whose locations are depicted as 23 vertical bars, including three thin red bars at locations where no *wtf* genes exist in the reference genome and four thick red bars at locations where compared to the reference genome one extra *wtf* gene is found. The *wtf* genes in the CBS5557 genome are named with the prefix *cw* and a number from the name of the syntenic gene in the reference genome. Two genes of a tandem pair corresponding to a singleton in the reference genome are distinguished using the suffixes a and b. Among the three genes of the triplet, two are named *cw11* and *cw12* based on their homology to *wtf11* and *wtf12*, respectively, and the gene situated between *cw11* and *cw12* is named *cw11x*. Genes at new locations are named *cw26*, *cw27*, and *cw28*. Chromosome lengths are not drawn to scale. (**B**) Diagrams depicting the eight genomic locations with more than one *wtf* gene in at least one of the two genomes. Genes are shown as arrows. Introns are not shown. *wtf17* is depicted not according to its annotation at PomBase, with its 5' boundary revised based on sequence alignment. (**C**) Maximum likelihood phylogenetic tree of 57 *wtf* genes of the reference genome and the CBS5557 genome. DNA sequences including the conserved_up regions, predicted coding sequences, and associated introns were aligned using the L-INS-i iterative refinement algorithm of MAFFT ([@bib26]) ([Figure 5---source data 1](#SD1-data){ref-type="supplementary-material"}). Maximum likelihood analysis was performed using IQ-TREE ([@bib43]). The tree was rooted by midpoint rooting ([@bib20]). Red dots on nodes indicate IQ-tree-calculated ultrafast bootstrap (UFBoot) support values \< 95%. Colored rectangles highlight phylogenetic neighbors. Two genes are considered phylogenetic neighbors if they are separated by a single internal node with a support value \>= 95%. Green rectangles indicate the 11 pairs of neighbors each composed of a reference *wtf* gene and a syntenic CBS5557 *wtf* gene. Yellow rectangles indicate the five pairs each composed of two *wtf* genes locating at different genomic positions. Magenta brackets denote the 12 genes that share a 150 bp sequence within the predicted intron 1 (see [Figure 5---figure supplement 3](#fig5s3){ref-type="fig"}). These 12 genes are divided into two subtypes, Intron-1-ATG genes and Exon-2-ATG genes, based on the locations of the first ATG codons downstream of the 150-bp-long sequence. The brown dashed box indicates genes with six exons. Scale bar, 0.1 nucleotide substitutions per nucleotide site.**DOI:** [http://dx.doi.org/10.7554/eLife.26057.016](10.7554/eLife.26057.016)10.7554/eLife.26057.017Figure 5---source data 1.MAFFT-aligned DNA sequences of 57 *wtf* genes of the reference genome and the CBS5557 genome.The sequences include the conserved_up regions, predicted coding sequences, and associated introns.**DOI:** [http://dx.doi.org/10.7554/eLife.26057.017](10.7554/eLife.26057.017)10.7554/eLife.26057.018Figure 5---source data 2.Gene structure predictions for 57 *wtf* genes of the reference genome and the CBS5557 genome.**DOI:** [http://dx.doi.org/10.7554/eLife.26057.018](10.7554/eLife.26057.018)10.7554/eLife.26057.019Figure 5---source data 3.MAFFT-aligned amino acid sequences of the predicted protein products of 57 *wtf* genes of the reference genome and the CBS5557 genome.**DOI:** [http://dx.doi.org/10.7554/eLife.26057.019](10.7554/eLife.26057.019)10.7554/eLife.26057.020Figure 5---figure supplement 1.Schematics of how the three CBS5557-only singleton *wtf* genes may have been lost in the laboratory strain through LTR-mediated recombination.Shown at the bottom of each diagram is the breakpoint junction sequence with 100% identity between the two directly oriented LTRs flanking the CBS5557-only singleton *wtf* gene and its position in the reference genome.**DOI:** [http://dx.doi.org/10.7554/eLife.26057.020](10.7554/eLife.26057.020)10.7554/eLife.26057.021Figure 5---figure supplement 2.Schematics depicting the high levels of diversity between five singleton CBS5557 *wtf* genes and their counterparts in the reference genome.Genes are shown as maroon arrows. Introns are not shown. Sequence alignment was performed using MAFFT via Jalview. Pair-wise identity was calculated by dividing the numbers of identical bases by the alignment length.**DOI:** [http://dx.doi.org/10.7554/eLife.26057.021](10.7554/eLife.26057.021)10.7554/eLife.26057.022Figure 5---figure supplement 3.*cw9*, *cw27*, and 10 other *wtf* genes in the reference genome and the CBS5557 genome share a 150 bp conserved sequence in the predicted intron 1.Red boxes denote the 150 bp sequence. The transcription start site of a short *wtf23* isoform identified by the Iso-Seq method ([@bib28]) is situated closely downstream of this 150 bp sequence, and at an optimal distance from a canonical TATA box (TATAAA) ([@bib32]). The first ATG codons downstream of the 150 bp sequence are located within the predicted intron 1 of the top five genes, which we call Intron-1-ATG genes, and within the predicted exon 2 of the bottom 7 genes, which we call Exon-2-ATG genes. These ATG codons are all in frame with the predicted coding sequences. Proteins initiated from these ATG codons correspond to the poison isoforms in the model proposed by [@bib46]). Sequence coordinates and alignment are based on the MAFFT-aligned DNA sequences in [Figure 5---source data 1](#SD1-data){ref-type="supplementary-material"}.**DOI:** [http://dx.doi.org/10.7554/eLife.26057.022](10.7554/eLife.26057.022)10.7554/eLife.26057.023Figure 5---figure supplement 4.The predicted protein products of 6-exon-containing *wtf* genes in the reference genome and the CBS5557 genome.Alignment is based on the MAFFT-aligned protein sequences in [Figure 5---source data 3](#SD3-data){ref-type="supplementary-material"}. Protein products of *wtf8*, *wtf22*, *wtf24*, and *cw16*, which belong to the 6-exon group but appear to have suffered pseudogenizing mutations, are not shown. Green brackets denote the eight transmembrane helices predicted by PolyPhobius (<http://phobius.sbc.su.se/poly.html>) ([@bib29]). Helices 3 and 8 are indicated with light green brackets, because the former is absent in some of the aligned sequences, and the latter has low residue-wise posterior probabilities in the PolyPhobius prediction.**DOI:** [http://dx.doi.org/10.7554/eLife.26057.023](10.7554/eLife.26057.023)10.7554/eLife.26057.024Figure 5---figure supplement 5.The predicted protein products of 5-exon-containing *wtf* genes in the reference genome and the CBS5557 genome.Alignment is based on the MAFFT-aligned protein sequences in [Figure 5---source data 3](#SD3-data){ref-type="supplementary-material"}. Protein products of the eight most divergent genes (*wtf7*, *wtf11*, *wtf14*, *wtf15*, *cw7*, *cw11*, *cw14*, and *cw15*) and eight genes that appear to have suffered pseudogenizing mutations (*wtf1*, *wtf2*, *wtf3*, *wtf6*, *wtf12*, *wtf17*, *cw3*, and *cw12*) are not shown. Green brackets denote the five transmembrane helices predicted by PolyPhobius. Helix three is indicated with a light green bracket, because it is absent in some of the aligned sequences.**DOI:** [http://dx.doi.org/10.7554/eLife.26057.024](10.7554/eLife.26057.024)

Among the 12 CBS5557 singleton *wtf* genes with syntenic counterparts in the reference genome, *cw9* and five others exhibit exceptionally high DNA sequence diversity (\<85% identity) from their counterparts in the reference genome ([Figure 1D](#fig1){ref-type="fig"}, [Figure 1---figure supplement 2](#fig1s2){ref-type="fig"}, and [Figure 5---figure supplement 2](#fig5s2){ref-type="fig"}). Given that the genome-wide nucleotide difference between CBS5557 and the laboratory strain is only 3.1 SNPs/kb ([@bib22]), *wtf* genes must have diverged at a much faster pace than average genes.

To examine the phylogenetic relationship among all 57 *wtf* genes of the two genomes, we performed a maximum likelihood analysis ([Figure 5C](#fig5){ref-type="fig"}). In the resulting phylogenetic tree, we used colored rectangles to highlight 16 pairs of genes, with each pair consisting of phylogenetic neighbors, defined as two genes separated by a single internal node that has an ultrafast bootstrap support value \>= 95% ([Figure 5C](#fig5){ref-type="fig"}). Interestingly, only 11 of these 16 pairs conform to typical orthologous relationship where two members of a pair are syntenic (highlighted by green rectangles in [Figure 5C](#fig5){ref-type="fig"}), whereas the five other pairs each consist of two genes that are non-syntenic (highlighted by yellow rectangles in [Figure 5C](#fig5){ref-type="fig"}). This pattern suggests that ectopic gene conversion (also called non-allelic or interlocus gene conversion) may contribute to the inter-isolate difference of *wtf* genes ([@bib8]; [@bib50]). Interestingly, genes belonging to the same gene cluster (pair or triplet), despite being subject more strongly to the homogenizing effect of gene conversion due to their physical proximity, are without exception very distinct from each other. One possible explanation is that diversifying selection may have prevented sequence homogenization within gene clusters.

The substantial difference in the numbers and sequences of *wtf* gene between the laboratory strain and CBS5557 suggests that rapid sequence divergence and independent loss/gain of *wtf* genes have resulted in different *S. pombe* natural isolates harboring distinct sets of active spore killers, which become reproductive barriers during inter-isolate crosses and contribute to the extensive and severe hybrid sterility within this species.

Discussion {#s3}
==========

Only a small number of gamete killers have been characterized at the molecular level. Two well-studied gamete killers in animals, *Segregation Distorter* in *Drosophila* and the *t* haplotype in mouse, both contain at least two key loci, a killer locus and a target locus ([@bib34]; [@bib7]; [@bib30]; [@bib33]; [@bib3]). In fungi, spore killers of *Neurospora* also appear to comprise of two genes, *rsk* (*resistant to Spore killer*) and a yet-to-be-identified killer gene, which are proposed to act together in a fashion similar to the toxin-antitoxin (TA) systems of bacterial plasmids, in which a stable toxin and a labile antitoxin together ensure that bacterial cells cannot lose the plasmid expressing them ([@bib18]; [@bib39]; [@bib55]). The involvement of more than one gene renders these gamete killer systems vulnerable to recombination that may separate the component genes, and as a consequence, they have all evolved recombination-suppressing features such as chromosomal inversions ([@bib7]; [@bib19]). On the other hand, two types of single-gene spore killers, \[Het-s\] prion and the *Spok* genes, have been found in the filamentous fungus *Podospora anserina* ([@bib11]; [@bib16]). The killing activity of the \[Het-s\] prion requires the presence of the HET-S protein, whereas the *Spok* genes can act autonomously. Fission yeast *wtf* killers share two prominent features with the *Spok* killers: being composed of a single predicted coding sequence and acting in a self-sufficient manner. These features arguably make *wtf* and *spok* killers the most compact and efficient forms of gamete killers.

Despite the similarities they share, *wtf* and *spok* killers are very different entities. Firstly, *wtf* genes only exist in a single fission yeast species *S. pombe*, indicating that they probably have arisen after the divergence of *S. pombe* from its closest relatives *S. octosporus* and *S. cryophilus* around 120 million years ago ([@bib53]), and have not jumped across the species barriers. In contrast, *Spok* family genes have a patchy but widespread distribution among Pezizomycotina species ([@bib16]), whose last common ancestor dates to 400 million years ago ([@bib51]), suggesting that *Spok* killers have originated much earlier than *wtf* killers or have been especially prone to transfer horizontally across species. Secondly, *spok* family genes encode soluble proteins with a C-terminal kinase-like domain (InterPro:IPR011009), whereas *wtf* genes are predicted to encode multi-transmembrane proteins, suggesting that the detailed molecular mechanisms of these two types of killers are likely to be quite distinct.

The single-gene nature of *wtf* killers begs the question: how does a single gene encode both the killing and the protecting activities? Our findings on the mutations that can disrupt one but not the other activity suggest the possibility that more than one type of expression products can be generated from a *wtf* killer gene. Consistent with this idea, Nuckolls et al. showed that the *wtf4* gene from the *S. pombe var. kambucha* strain (*Sk wtf4*) acts as a spore killer by expressing two overlapping transcripts, with the shorter transcript encoding a spore-killing poison and the longer transcript encoding an antidote that protects the spores from the killing effect of the poison ([@bib46]). This dual-transcript model can satisfactorily explain why the Tb and Te forms of *cw9* and *cw27* lost protecting activities but not killing activities ([Figure 4A](#fig4){ref-type="fig"}). Presumably, the transcription of the longer isoforms, which encode antidote proteins initiated from a start codon in the predicted exon 1, partially relies on the sequence between the upstream LTR and the conserved_up sequence and absolutely requires the conserved_up sequence; in contrast, the transcription of the shorter isoforms, which encode poison proteins initiated from a downstream alternative start codon, is not affected by removing the upstream sequences.

We hypothesized that the transcription of the poison isoforms may depend on a special promoter residing within the predicted intron 1 of *cw9* and *cw27*. Indeed, through inspecting the DNA sequence alignment ([Figure 5---source data 1](#SD1-data){ref-type="supplementary-material"}), we found that *cw9*, *cw27*, and 10 other *wtf* genes in the reference genome and the CBS5557 genome share a 150 bp conserved sequence in the predicted intron 1, and this sequence is closely upstream of the transcription start site of a shorter-than-predicted *wtf23* transcript isoform uncovered by the Iso-Seq long-read RNA sequencing method ([Figure 5---figure supplement 3](#fig5s3){ref-type="fig"}) ([@bib28]). These 12 genes can be classified into two subtypes, which we call Intron-1-ATG genes (*wtf4*, *wtf13*, *cw4*, *cw18*, and *cw23*) and Exon-2-ATG genes (*wtf19*, *wtf23*, *cw8a*, *cw9*, *cw21*, *cw22*, and *cw27*), based on the locations of the first ATG codons downstream of the 150 bp sequence ([Figure 5C](#fig5){ref-type="fig"} and [Figure 5---figure supplement 3](#fig5s3){ref-type="fig"}). These ATG codons are all in frame with the predicted coding sequences. The within-intron-1 ATG codons in the 5 Intron-1-ATG genes are only 7 bp away from the predicted exon 2, and if used as start codons, would result in proteins beginning with three intron-1-coded amino acids (Met, Leu, and Ser). The two active spore killer genes found in the *S. pombe var. kambucha* strain by Nuckolls et al., *Sk wtf4* and *Sk wtf28*, also possess this 150 bp sequence and can be classified as an Intron-1-ATG gene and an Exon-2-ATG gene, respectively. We propose that this 150 bp sequence is a hallmark of the 'poison-and-antidote' genes that can, at least potentially, generate both poison and antidote isoforms according to the dual-transcript model proposed by Nuckolls et al.

The other 45 *wtf* genes in the reference genome and the CBS5557 genome lack this 150 bp sequence. In addition, except for the eight most divergent genes (*wtf7*, *cw7*, *wtf11*, *cw11*, *wtf14*, *cw14*, *wtf15*, and *cw15*), they do not have an in-frame ATG codon in intron one or in the 5'-section of exon 2, suggesting that most of these genes are not able to generate a shorter isoform according to the dual-transcript model, and are thus 'antidote-only' genes. Alternatively, it is possible that they may generate two different protein products using other mechanisms, such as alternative splicing ([@bib28]), overlapping proteins translated in different frames ([@bib12]), and post-translational processing and modification.

Five *wtf* genes in the reference genome, *wtf1*, *wtf2*, *wtf3*, *wtf22*, and *wtf24*, are annotated by PomBase as pseudogenes ([@bib38]). Our analysis of the gene structures of the 57 *wtf* genes in the reference genome and the CBS5557 genome showed that in addition to these five genes, four other *wtf* genes in the reference genome (*wtf6*, *wtf8*, *wtf12*, and *wtf17*), and three *wtf* genes in the CBS5557 genome (*cw3*, *cw12*, and *cw16*), appear to have suffered pseudogenizing mutations (see Materials and methods and [Figure 5---source data 2](#SD2-data){ref-type="supplementary-material"}), which include loss of start codons (*wtf3* and *cw3*), premature stop codons (*wtf1*, *wtf12*, *cw12*, and *cw16*), frameshifting indels (*wtf6*, *wtf8*, *wtf22*, and *wtf24*), and segmental deletions (*wtf2* and *wtf17*). These 'dead' *wtf* genes presumably are not able to generate either antidote or poison proteins.

Most of the *wtf* genes in the reference genome and the CBS5557 genome are predicted to contain either 5 or 6 exons. Interestingly, the 6-exon genes, including *cw9*, *cw27*, and all but one of the other ten 150-bp-sequence-containing genes, clustered together in the DNA sequence-based phylogenetic tree ([Figure 5C](#fig5){ref-type="fig"}), and their protein products are predicted to be longer and contain two or three more transmembrane helices than the protein products of 5-exon genes ([Figure 5---source data 3](#SD3-data){ref-type="supplementary-material"}, [Figure 5---figure supplement 4](#fig5s4){ref-type="fig"}, and [Figure 5---figure supplement 5](#fig5s5){ref-type="fig"}). For each of the 12 150-bp-sequence-containing genes, the antidote protein product is approximately 50 amino acids longer than the poison protein product at the N-terminal soluble region. The 10 C-terminal amino acids missing from the protein products of the Td forms of *cw9* and *cw27* are downstream of the last transmembrane helix, and are strongly conserved among the predicted protein products of 6-exon genes ([Figure 5---figure supplement 4](#fig5s4){ref-type="fig"}), indicating that Td-like mutations do not occur naturally to a significant extent.

How do transmenbrane proteins act as killers? In the prokaryotic toxin-antitoxin (TA) systems, many toxins are transmembrane proteins and kill cells by damaging cellular membranes ([@bib56]). The best-known examples include the type I TA toxins Hok and TisB, and the type V TA toxin GhoT ([@bib9]; [@bib15]; [@bib63]). We speculate that Wtf proteins may also possess abilities to disrupt the integrity of certain cellular membranes.

*wtf* genes were named after their association with the solo LTRs of the Tf retrotransposons ([@bib67]). It was initially hypothesized that retrotransposition may aid the expansion of the *wtf* gene family ([@bib67]). However, the presence of introns in *wtf* genes and the lack of sequence similarity between *wtf*-flanking LTRs indicate that *wtf* genes have not been retrotransposed ([@bib5]). Bowen et al. have put forward an explanation for the enrichment of LTRs near *wtf* genes: the regions flanking *wtf* genes may have been favored insertion targets of Tf retrotransposons, perhaps owing to the high transcription levels of *wtf* genes during meiosis and sporulation ([@bib5]). Our findings that LTRs flanking *cw9* and *cw27* are dispensable for spore killing support the idea that the *wtf*-LTR association does not benefit *wtf* genes but may rather reflect an exploitation of *wtf* genes by Tf retrotransposons for their own expression and/or inheritance benefit. Furthermore, we show that three singleton *wtf* genes (*cw26*, *cw27*, and *cw28*) appear to have been lost in the laboratory strain due to LTR-mediated recombination, suggesting that the *wtf*-LTR association promotes *wtf* gene loss and is detrimental to the long-term persistence of *wtf* genes. In addition, it may not be coincidental that both *cw9* and *cw27* require an approximately 90-bp-long sequence between the upstream LTR and the conserved_up sequence for the full protecting activity, raising the possibility that some *wtf* genes that lack such an intervening sequence (e.g. *wtf4* and *cw4*) may have suffered functional loss due to a Tf insertion too close to the conserved_up sequence or may have become functionally dependent on the upstream LTR.

In terms of the number of genes per genome, *wtf* genes are probably the most evolutionarily successful gamete killers known to date. During a relatively short evolutionary time span, *wtf* gene family has emerged and become the largest species-specific gene family in *S. pombe* ([@bib31]). The compactness (minimal sizes of fully functional *cw9* and *cw27* are both around 1.9 kb) and the self-sufficiency of the *wtf* killers likely have been advantageous to its family size expansion. Furthermore, given the dramatic divergence of *wtf* genes between the reference genome and the CBS5557 genome, and the mutual-killing behaviors of *cw9* and *cw27*, it is possible that two identical copies of a *wtf* gene generated by a gene duplication event can undergo relatively rapid sequence diversification to become mutually non-resistant, so that their selective advantage is not diminished by resistance. Ectopic gene conversion is known to play a causal role in the sequence diversification of animal major histocompatibility complex (MHC) genes ([@bib47]; [@bib61]), plant disease-resistance genes (R genes) ([@bib41]), and *Plasmodium falciparum var* genes ([@bib10]). Our phylogenetic analysis suggests that ectopic gene conversion may also contribute to the historical expansion and the ongoing dynamic change of the *wtf* gene family.

It is intriguing that most of the *wtf* genes locate on chromosome III, the smallest of the three fission yeast chromosomes. To explain this phenomenon, Bowen et al. proposed that the laboratory strain may have inherited only chromosome III but not the other two chromosomes from an ancestor strain in which *wtf* genes were abundant on all three chromosomes ([@bib5]). Our finding here of the same kind of chromosomal distribution bias in CBS5557, whose nucleotide diversity from the laboratory strain is higher on chromosome III than on the other two chromosomes ([@bib22]), makes this explanation less likely. The only type of viable aneuploidy in fission yeast is chromosome III disomy ([@bib45]; [@bib44]), which occurs at a frequency of \>1% during normal meiosis ([@bib40]). We speculate that the tolerance of chromosome III disomy may in some way be advantageous to the expansion of the *wtf* gene family on this chromosome.

Spore killers and other meiotic drivers, despite having only been found in a small number of taxa, are believed to exist in a wide range of species, where they are an important force influencing genome organization, gametogenesis, and speciation ([@bib7]; [@bib33]; [@bib66]). The molecular basis of their action and the evolutionary mechanisms behind their birth, spread, and extinction are scarcely known. The unveiling of the molecular identities of spore killers in *S. pombe* by this study and that of Nuckolls et al. ([@bib46]), promises to bring the power of a highly tractable model organism to the pursuit of the many mysteries of these evolutionary wonders.

Materials and methods {#s4}
=====================

Strains {#s4-1}
-------

Fission yeast strains used in this study are listed in [Supplementary file 2](#SD5-data){ref-type="supplementary-material"} in Microsoft Word format. Genetic methods and composition of media are as described ([@bib14]). Gene deletion was constructed by PCR-based gene targeting. Integrating plasmids were based on pDUAL and related vectors ([@bib35], [@bib36]). A pDUAL-based plasmid was linearized with NotI digestion and integrated at the *leu1* locus, or linearized with MluI digestion and integrated at the *ars1* replication origin region upstream of the *hus5* gene.

To convert homothallic strains to heterothallic and to mark the mating type locus, we created the *mat1Δ17* mutation ([@bib1]), by replacing a 140 bp sequence between the H1 homology box at the *mat1* locus and the nearby SspI restriction site with an antibiotic resistance marker using PCR-based gene targeting. *h^+^/h^-^* diploid strains were constructed by mating on SPAS medium *h^+^* and *h^-^* strains with different antibiotic resistance markers inserted at the *mat1* locus, and a few hours later restreaking the mating mixture onto YES medium containing two types of antibiotics.

To construct strains with only a portion of chromosome III originating from CBS5557 and the rest of the genome coming from the laboratory strain, we first introduced *rec12∆* mutation, which blocks meiotic recombination, into the laboratory strain and CBS5557, respectively. Through crossing the two resultant strains with each other, we obtained a strain whose chromosome I and II are from the laboratory strain and whose chromosome III is from CBS5557. This strain was then crossed to a *rec12^+^* laboratory-background strain carrying an antibiotic resistance marker at a locus on the left arm of chromosome III (*mug123:hphMX*) and another antibiotic resistance marker at a locus on the right arm of chromosome III (*rps2802:natMX*) (first backcross). Progeny strains with only one of the loci marked were kept. They were then crossed to laboratory-background strains with the other locus marked by *kanMX* (second backcross). By scoring the marker segregation ratios of the second backcross, we selected those first-backcross progenies that conferred strong segregation distortion at the *kanMX*-marked locus. Their second-backcross progenies that harbored the *kanMX* marker were tested by crossing, and those that conferred segregation distortion at the *kanMX*-marked locus were selected. One of the second-backcross progenies that conferred segregation distortion at the *mug123* locus is the backcrossed-1 strain used in the bulk segregant analysis shown in [Figure 1B](#fig1){ref-type="fig"}. Two chromosome III genomic regions from CBS5557, one between coordinates 247743 and 568805, and the other between coordinates 814348 and 1711635, are present in the backcrossed-1 strain. Another second-backcross progeny that conferred segregation distortion at the *rps2802* locus underwent additional rounds of backcross to yield the backcrossed-2 strain used in the bulk segregant analysis shown in [Figure 1C](#fig1){ref-type="fig"}. The CBS5557 chromosome III genomic region present in the backcrossed-2 strain is between coordinates 1839246 and 1916857.

Bulk segregant analysis {#s4-2}
-----------------------

For Illumina sequencing-assisted bulk segregant analysis, \>150 viable progeny colonies derived from random spore analysis were pooled together. Genomic DNA extraction, sequencing libraries construction, and Illumina sequencing were performed as described ([@bib22]). Sequencing data were deposited at NCBI SRA under the following accession numbers: viable progenies from CBS5557 × laboratory strain (DY9974 × DY8531), SRR5131575; viable progenies from backcrossed-1 × laboratory strain (cross-derived diploid strain was named DY26097), SRR5131578; viable progenies from backcrossed-1 *cw9∆* × laboratory strain (cross-derived diploid strain was named DY26100), SRR5131579; viable progenies from backcrossed-2 strain × laboratory strain (cross-derived diploid strain was named DY26095), SRR5131577; viable progenies from backcrossed-2 *cw27∆* strain × laboratory strain (cross-derived diploid strain was named DY26092), SRR5131576. Sequencing reads were mapped to the reference genome using BWA-MEM (version 0.7.15-r1140, RRID:[SCR_010910](https://scicrunch.org/resolver/SCR_010910)). SAM files were converted to BAM files and duplicates were removed using SAMtools (version 0.1.19--44428 cd, RRID:[SCR_002105](https://scicrunch.org/resolver/SCR_002105)). The average read depths were 14.9× for DY9974 × DY8531, 10.4× for DY26097, 18.4× for DY26100, 16.2× for DY26095, and 23.1× for DY26092. The software bam-readcount version 0.7.4 was run with the option -q 60 -b 35 (-q 60 -b 27 for DY26097, which had lower depth coverage and lower sequencing quality) to obtain the read counts of four different bases at each of the 38783 SNP positions identified previously ([@bib22]). Reference allele frequency at each SNP position was calculated by dividing the read count of the reference base by the sum of the read counts of the reference base and the variant base. Only SNP positions with the sum of the read counts of the reference base and the variant base \>= 10 (\> = 7 for DY26097) were used in the plots. A trend line was drawn in each plot by calculating a rolling median of reference allele frequencies at 45 consecutive SNP positions using the rollapply function of R's zoo package (version 1.7--13).

PacBio sequencing {#s4-3}
-----------------

For PacBio sequencing of CBS5557 genome, we prepared genomic DNA from DY9971, an *h^-^* derivative of CBS5557, by first grinding cells in liquid nitrogen and then extracting DNA using the Maxi Column Fungal DNAOUT kit (Tiandz, Beijing, China). Sequencing library construction and sequencing on the PacBio RS II platform using P6C4 chemistry were performed by BGI (Shenzhen, China). After filtering out low-quality reads and removing adaptor sequences, we obtained 91,333 subreads no shorter than 1000 bp (mean length 7370 bp). These reads were used for de novo assembly by the SMRT Analysis software. The assembly was polished using Pilon version 1.21 (RRID:[SCR_014731](https://scicrunch.org/resolver/SCR_014731)) together with our previously published CBS5557 Illumina sequencing data ([@bib22]; [@bib64]). PacBio sequencing data were deposited at NCBI SRA under the accession number SRR5133273.

Electron microscopy analysis {#s4-4}
----------------------------

For the induction of synchronous meiosis and sporulation, *h^+^/h^-^* diploid strains were pre-grown in liquid YES medium to log phase, cultured in SSL+N synthetic liquid medium for 10 hr, and then shifted to nitrogen-free SSL-N medium ([@bib13]). Four-spored asci represented about half of the cells 10 hr after shifting to SSL-N ([@bib48]). Electron microscopy was performed as described ([@bib60]).

Spore viability analysis {#s4-5}
------------------------

Spore viability was assessed by tetrad analysis using a TDM50 tetrad dissection microscope (Micro Video Instruments, Avon, USA). At least 40 tetrads (160 spores) were analyzed for each cross. Based on empirical evidence, such a sample size is large enough for reliably assessing spore viability. Numerical data of the tetrad analysis in Excel format are provided as [Supplementary file 1](#SD4-data){ref-type="supplementary-material"}. For statistical analysis of the spore viability data, Fisher's exact test and exact binomial test of goodness-of-fit were performed using Excel spreadsheets downloaded from <http://www.biostathandbook.com/fishers.html> and <http://www.biostathandbook.com/exactgof.html>, respectively ([@bib37]).

Sequence and phylogenetic analysis {#s4-6}
----------------------------------

Except for the phylogenetic analysis, DNA sequences were aligned using MAFFT via Jalview ([@bib25]; [@bib65]). Pair-wise identity was calculated by dividing the numbers of identical bases by the alignment length using the Sequence Manipulation Suite web server (<http://www.bioinformatics.org/sms2/ident_sim.html>) ([@bib59]).

For the phylogenetic analysis shown in [Figure 5C](#fig5){ref-type="fig"}, an alignment of DNA sequences containing the conserved_up regions, the coding sequences, and the associated introns was generated using the L-INS-i algorithm of MAFFT version 7.310 for Mac OS X (RRID:[SCR_011811](https://scicrunch.org/resolver/SCR_011811)) with the command line option \--localpair \--maxiterate 16 \--reorder ([@bib26]), and manually adjusted to correct an obvious misalignment of the *wtf2* sequence. The final alignment of 57 *wtf* genes of the reference and CBS5557 genome in FASTA format is provided as [Figure 5---source data 1](#SD1-data){ref-type="supplementary-material"}. Maximum likelihood analysis was performed using IQ-TREE version 1.5.3 for Mac OS X with the command line option -m TEST -alrt 1000 -bb 1000 ([@bib43]). The tree was rooted by midpoint rooting using FigTree version 1.4.2 (<http://tree.bio.ed.ac.uk/software/figtree/>) and was visualized using Phylo.io ([@bib54]).

Gene structures of the 57 *wtf* genes were predicted using the AUGUSTUS web server (version 3.2.3, RRID:[SCR_008417](https://scicrunch.org/resolver/SCR_008417)) with the parameter settings \--genemodel=exactlyone \--sample=100 \--keep_viterbi=true \--alternatives-from-sampling=true \--minexonintronprob=0.08 \--minmeanexonintronprob=0.3 \--maxtracks=20 ([@bib58]). When multiple predications were made for a gene, we manually inspected the predictions and selected the one that best conforms to conserved patterns of intron positions. For the majority of genes, the selected prediction is the top-scoring prediction, and for 10 genes (*wtf6*, *wtf7*, *wtf14*, *wtf24*, *cw6a*, *cw7*, *cw12*, *cw14*, *cw16*, *cw25*), the selected prediction is the one with the second best score. The final predictions for all 57 genes in plain text format are provided as [Figure 5---source data 2](#SD2-data){ref-type="supplementary-material"}, which can be loaded into Jalview as an annotation file for the DNA sequence alignment in [Figure 5---source data 1](#SD1-data){ref-type="supplementary-material"}. The failure to predict an exon flanked by conserved intron-exon junction sequences is an indication that pseudogenizing mutation(s) exist. Among the predicted protein sequences obtained based on these gene structure predictions, the sequences of Wtf4, Wtf6, Wtf7, Wtf12, Wtf13, Wtf15, Wtf17, and Wtf18 are different from those annotated at PomBase. A MAFFT (L-INS-i algorithm)-generated alignment of the 57 predicted protein sequences in FASTA format is provided as [Figure 5---source data 3](#SD3-data){ref-type="supplementary-material"}.

Data availability {#s4-7}
-----------------

Illumina and PacBio sequencing data reported in this paper are available at the Sequence Read Archive (SRA) under the accession number SRP095878 (BioProject PRJNA358837) ([@bib23]). The DNA sequences and annotations of the 32 *wtf* genes of the CBS5557 genome have been deposited at GenBank under the accession numbers KY926712-KY926743.

Funding Information
===================

This paper was supported by the following grant:

-   http://dx.doi.org/10.13039/501100002855Ministry of Science and Technology of the People\'s Republic of China National Basic Research Program of China 2014CB849901 to Li-Lin Du.

We thank Sarah Zanders for communicating results prior to publication. We are grateful to Meng-Qiu Dong for critical reading of the manuscript. This work was supported by the National Basic Research Program of China (973 Program, 2014CB849901).

Additional information {#s5}
======================

WH: A patent application related to this work has been filed.

FS: A patent application related to this work has been filed.

L-LD: A patent application related to this work has been filed.

The other authors declare that no competing interests exist.

WH, Conceptualization, Investigation, Writing---original draft.

Z-DJ, Investigation.

FS, Formal analysis, Writing---original draft.

J-XZ, Investigation.

W-ZH, Formal analysis.

L-LD, Conceptualization, Funding acquisition, Writing---original draft, Writing---review and editing.

Additional files {#s6}
================

10.7554/eLife.26057.025

###### Numerical data of the tetrad analysis.

**DOI:** [http://dx.doi.org/10.7554/eLife.26057.025](10.7554/eLife.26057.025)

10.7554/eLife.26057.026

###### Fission yeast strains used in this study.

**DOI:** [http://dx.doi.org/10.7554/eLife.26057.026](10.7554/eLife.26057.026)

Major datasets {#s7}
--------------

The following dataset was generated:

Hu W,Suo F,Du LL,2017,Schizosaccharomyces pombe Genome sequencing,<https://www.ncbi.nlm.nih.gov/bioproject/PRJNA358837>,Publicly available at the NCBI BioProject (accession no: PRJNA358837)

10.7554/eLife.26057.029

Decision letter

Rokas

Antonis

Reviewing editor

Vanderbilt University

,

United States

In the interests of transparency, eLife includes the editorial decision letter and accompanying author responses. A lightly edited version of the letter sent to the authors after peer review is shown, indicating the most substantive concerns; minor comments are not usually included.

Thank you for submitting your article \"A large gene family in fission yeast encodes spore killers that subvert Mendel\'s law\" for consideration by *eLife*. Your article has been favorably evaluated by Diethard Tautz (Senior Editor) and three reviewers, one of whom is a member of our Board of Reviewing Editors. The following individuals involved in review of your submission have agreed to reveal their identity: Nick Rhind (Reviewer \#2); Michael Lichten (Reviewer \#3).

The reviewers have discussed the reviews with one another and the Reviewing Editor has drafted this decision to help you prepare a revised submission.

Summary:

Hu et al. have identified, in the genome of *S. pombe* strain CBS5557 two members of the *wtf* repeated gene family that serve as meiotic drive elements when crossed to the reference strain. Using deletions of the two loci, they show that the protective function can be separated from the killing function, and show that the C-terminus of the putative open-reading frame is important for killing but not for protection. In addition, they show that the two genes actually kill their sister spores, using both viability and EM-based assays. They then analyze the ensemble of *wtf/cw* loci in CBS5557 and show that these elements appear to frequently \"move\" by ectopic gene conversion. Combined with Nuckolls et al., this paper makes an intriguing and interesting story.

Essential revisions:

1\) The switch in nomenclature, from *k1wtf* and *k2wtf* to *cw9* and *cw27*, is needlessly confusing. Why not simply call these loci *cw9* and *cw27* from the start, and just say that the reasons for doing so will be presented below? Along similar lines, the *wtf* nomenclature depicted in [Figure 5](#fig5){ref-type="fig"} seems unnecessarily complicated. In situations in which similarity and synteny make homology unambiguous, ambiguous nomenclature is not useful. For instance, *cw\[3-4\]a* and *cw\[3-4\]b* should be named *cw3* and *cw4*, respectively, likewise the homologs of *wtf11, 12, 14, 15, 17* and *18* are unambiguous.

2\) [Figure 1B/C](#fig1){ref-type="fig"} is confusing, especially since it hides the fact that there actually are no SNPs segregating in the regions where the red line is at 1. Why not simply display the rolling median for regions where there actually are SNPs? The first and last 40 SNPs will not be displayed, but there appear to be so many that this will not materially affect the presentation. The terms \"left-peak\" and \"right-peak\" strain are also confusing; we suggest adding a panel that shows a chr 3 map for the two strains with the reference and CBS5557-introgressed sequences in different colors, and including that color code on the X-axis of the plots in panel B.

3\) It\'s not clear why *k1wft*∆ does not restore reference allele recovery to 50% in the left-peak strain x laboratory strain cross ([Figure 1B](#fig1){ref-type="fig"}). Is it possible that *cw\[11-12\]c* also has drive activity?

4\) [Figure 5](#fig5){ref-type="fig"}, panel C would be considerably more useful if it contained a horizontal scale, to roughly indicate the extent of divergence between different alleles. Criteria for considering two alleles to be paired should also be indicated; for example, it\'s not clear why *cw8a* and *cw24* are considered to be a homeologous pair, but *wtf24* and *cw6b* are not. Furthermore, the phylogeny should be shown as an unrooted tree -- the authors reconstructed the phylogeny as unrooted and they should depict it that way too.

5\) For readers interested in further numerical analysis, it would be useful if the actual numerical data underlying all of the bar graphs were included in a supplementary file. For example, we were wondering if the loss of spore viability in *k1*∆ *k2*∆ double mutant spores ([Figure 2](#fig2){ref-type="fig"}) is greater than what one would expect on an additive basis of spore viability loss in the two single mutant spore types-it looks as if it is, but with the current presentation it\'s impossible to check.

6\) It would be useful to use the translation start site-criterion outlined in Nuckolls et al. for whether the *cw* loci are \"dead\", \"antidote-only\", or \"antidote-toxin\". Similarly, since the C-termini of *cw9* and *cw27* appear to be important for toxin activity, is this region more conserved among the wtf loci than other parts of the protein? In \"antidote-only\" or in \"antidote-toxin\" loci?

7\) It is not clear why the variation of allele frequency depicted in [Figure 1](#fig1){ref-type="fig"} is so large. If over 150 segregants were sequenced, how can some loci in regions that segregate 50:50 (e.g. on the left end of Chromosome I) be represented by only one allele?

8\) The authors should extend the discussion of their findings and take into account the findings of Nuckolls et al. In particular:

What do the other *wtf* genes do? Could they be acting as rescuers but not as killers? Could they have assumed other unrelated functions?

[Figure 5A](#fig5){ref-type="fig"} shows an unusual concentration of *wtf* genes in chr III. Why is that the case? (Could it have to do with the fact that it is the only chromosome where disomy is tolerated?)

On the second page of the Discussion, the authors speculate as to how one locus could produce both a toxin and an antidote. It should be possible to eliminate, or at least deprioritize, one or more of the proposed mechanisms by inspecting the sequence of the loci and perhaps publicly available RNA-seq data. It would be of great interest to have some speculation on how a putative trans-membrane protein could function in a toxin-antidote system.

10.7554/eLife.26057.030

Author response

Essential revisions:

1\) The switch in nomenclature, from k1wtf and k2wtf to cw9 and cw27, is needlessly confusing. Why not simply call these loci cw9 and cw27 from the start, and just say that the reasons for doing so will be presented below? Along similar lines, the wtf nomenclature depicted in [Figure 5](#fig5){ref-type="fig"} seems unnecessarily complicated. In situations in which similarity and synteny make homology unambiguous, ambiguous nomenclature is not useful. For instance, cw\[3-4\]a and cw\[3-4\]b should be named cw3 and cw4, respectively, likewise the homologs of wtf11, 12, 14, 15, 17 and 18 are unambiguous.

We have followed the suggestions of the reviewers and changed the nomenclature of the *wtf* genes in the CBS5557 genome.

2\) [Figure 1B/C](#fig1){ref-type="fig"} is confusing, especially since it hides the fact that there actually are no SNPs segregating in the regions where the red line is at 1. Why not simply display the rolling median for regions where there actually are SNPs? The first and last 40 SNPs will not be displayed, but there appear to be so many that this will not materially affect the presentation. The terms \"left-peak\" and \"right-peak\" strain are also confusing; we suggest adding a panel that shows a chr 3 map for the two strains with the reference and CBS5557-introgressed sequences in different colors, and including that color code on the X-axis of the plots in panel B.

[Figure 1B and C](#fig1){ref-type="fig"} have been revised according to the suggestions of the reviewers. We have renamed the two backcrossed strains backcrossed-1 and backcrossed-2, respectively.

3\) It\'s not clear why k1wft∆ does not restore reference allele recovery to 50% in the left-peak strain x laboratory strain cross ([Figure 1B](#fig1){ref-type="fig"}). Is it possible that cw\[11-12\]c also has drive activity?

We agree with the reviewers that the residual allele frequency bias is probably due to additional active spore killer gene(s) in the CBS5557 genomic regions still present in the backcrossed-1 strain. We have added a sentence about this possible explanation in the Results section. There are 12 *cw* genes remaining in the backcrossed-1 strain. Future work will be needed to determine which one(s) are active killer genes.

4\) [Figure 5](#fig5){ref-type="fig"}, panel C would be considerably more useful if it contained a horizontal scale, to roughly indicate the extent of divergence between different alleles. Criteria for considering two alleles to be paired should also be indicated; for example, it\'s not clear why cw8a and cw24 are considered to be a homeologous pair, but wtf24 and cw6b are not. Furthermore, the phylogeny should be shown as an unrooted tree -- the authors reconstructed the phylogeny as unrooted and they should depict it that way too.

During the revision, we have improved the de novo assembly of CBS5557 genome by adding a polishing step using the software Pilon. With polished sequences of the *cw* genes, we have generated a new phylogenetic tree, which has a topology similar to, but not exactly the same as, the tree in the original manuscript. This new tree is shown in [Figure 5C](#fig5){ref-type="fig"} of the revised manuscript. We now include a scale bar in [Figure 5C](#fig5){ref-type="fig"}, as suggested by the reviewers. In the Results section and in the legend of [Figure 5](#fig5){ref-type="fig"}, we have now more precisely defined the criteria for highlighting a pair of genes, which we now term phylogenetic neighbors. Two genes are considered phylogenetic neighbors if they are separated on the phylogenetic tree by a single internal node with an IQ-tree-calculated ultrafast bootstrap support value \>= 95%. We agree with the reviewers that it is inappropriate to draw unrooted tree in a rooted format. We have now rooted the tree using midpoint rooting, which is a widely used method for rooting a phylogenetic tree when no suitable outgroup is available (see the reference Hess & De Moraes Russo, 2007).

5\) For readers interested in further numerical analysis, it would be useful if the actual numerical data underlying all of the bar graphs were included in a supplementary file. For example, we were wondering if the loss of spore viability in k1∆ k2∆ double mutant spores ([Figure 2](#fig2){ref-type="fig"}) is greater than what one would expect on an additive basis of spore viability loss in the two single mutant spore types-it looks as if it is, but with the current presentation it\'s impossible to check.

Numerical data are now provided in [Supplementary file 1](#SD4-data){ref-type="supplementary-material"} of the revised manuscript. The double mutant spores indeed appeared to have suffered a greater loss of viability than expected from a simple additive effect, but the exact reason is not clear at this moment.

6\) It would be useful to use the translation start site-criterion outlined in Nuckolls et al. for whether the cw loci are \"dead\", \"antidote-only\", or \"antidote-toxin\". Similarly, since the C-termini of cw9 and cw27 appear to be important for toxin activity, is this region more conserved among the wtf loci than other parts of the protein? In \"antidote-only\" or in \"antidote-toxin\" loci?

We have more carefully examined the sequences of the 57 *wtf* genes in the reference genome and the CBS5557 genome, and uncovered a within-intron-1 150-bp sequence shared by all currently known active *wtf* killers (*cw9* and *cw27* found in our study and *Sk wtf4* and *Sk wtf28* found by Nuckolls et al.). We added a new figure ([Figure 5---figure supplement 3](#fig5s3){ref-type="fig"}) about this sequence. We propose that this sequence may serve as a promoter driving the transcription of the poison isoforms and the presence of this sequence is a useful criterion for identifying \"poison-and-antidote\" genes. In the Discussion section of the revised manuscript, we added three paragraphs to discuss \"poison-and-antidote\" genes that share this 150-bp sequence, \"antidote-only\" genes lacking this 150-bp sequence, and the \"dead\" genes that have suffered pseudogenizing mutations, respectively.

Based on the predicted gene structures of the *wtf* genes, most of them have either 5 or 6 exons. All currently known active *wtf* killers are 6-exon genes. We have added two new figures ([Figure 5---figure supplement 4](#fig5s4){ref-type="fig"} and [Figure 5---figure supplement 5](#fig5s5){ref-type="fig"}) depicting the sequence alignments of the protein products of the 6-exon genes and the 5-exon genes, respectively. From the alignment, it can be seen that the C-terminal sequences of Cw9 and Cw27 proteins are strongly conserved among the protein products of 6-exon genes, regardless of whether they are \"poison-and-antidote\" genes or \"antidote-only\" genes.

7\) It is not clear why the variation of allele frequency depicted in [Figure 1](#fig1){ref-type="fig"} is so large. If over 150 segregants were sequenced, how can some loci in regions that segregate 50:50 (e.g. on the left end of Chromosome I) be represented by only one allele?

For cost reasons, we usually limit the Illumina sequencing of the bulk segregant analysis (BSA) samples to an average read depth of between 10× and 20×. Therefore, the read depth at any SNP loci is much lower than the number of segregants, and as a consequence, the allele frequency data shown in [Figure 1](#fig1){ref-type="fig"} are much noisier than expected from the number of segregants. We have now listed in the Materials and methods section the average read depths of the five BSA sequencing datasets.

*8) The authors should extend the discussion of their findings and take into account the findings of Nuckolls et al. In particular:*

What do the other wtf genes do? Could they be acting as rescuers but not as killers? Could they have assumed other unrelated functions?

We have added several paragraphs in the Discussion section to describe the classification the 57 *wtf* genes in the reference genome and the CBS5557 genome into \"poison-and-antidote\" genes, \"antidote-only\" genes, and \"dead\" genes.

[Figure 5A](#fig5){ref-type="fig"} shows an unusual concentration of wtf genes in chr III. Why is that the case? (Could it have to do with the fact that it is the only chromosome where disomy is tolerated?)

We agree with the reviewers that tolerance of chromosome III disomy may be a reason why *wtf* genes concentrate on chromosome III. We have added a paragraph in the Discussion section to discuss the possible reasons behind this intriguing phenomenon.

On the second page of the Discussion, the authors speculate as to how one locus could produce both a toxin and an antidote. It should be possible to eliminate, or at least deprioritize, one or more of the proposed mechanisms by inspecting the sequence of the loci and perhaps publicly available RNA-seq data. It would be of great interest to have some speculation on how a putative trans-membrane protein could function in a toxin-antidote system.

In the newly added paragraphs of the Discussion section, by more carefully inspecting the sequences, we have classified the 57 *wtf* genes in the reference genome and the CBS5557 genome into \"poison-and-antidote\" genes, which conform to the dual-transcript model of Nuckolls et al., \"antidote-only\" genes, and \"dead\" genes. We speculate that it remains possible that the \"antidote-only\" genes may be able to generate poison isoforms through mechanisms other than alternative transcription start sites. We have also added a paragraph in the Discussion section to speculate on how transmembrane proteins can act as spore killers.
