Bdelloid rotifers use hundreds of horizontally acquired genes against fungal pathogens

Obligately asexual lineages are typically rare and short-lived. According to one hypothesis, they adapt too slowly to withstand relentlessly coevolving pathogens. Bdelloid rotifers seem to have avoided this fate, by enduring millions of years without males or sex. We investigated whether bdelloids’ unusual capacity to acquire non-metazoan genes horizontally has enhanced their resistance to pathogens. We found that horizontally transferred genes are three times more likely than native genes to be upregulated in response to a natural fungal pathogen. This enrichment was twofold stronger than that elicited by a physical stressor (desiccation), and the genes showed little overlap. Among hundreds of upregulated non-metazoan genes were RNA ligases putatively involved in resisting fungal toxins and glucanases predicted to bind to fungal cell walls, acquired from bacteria. Our results provide evidence that bdelloids mitigate a predicted challenge of long-term asexuality in part through their ability to acquire and deploy so many foreign genes.


Introduction
The maintenance of sexual reproduction is a longstanding problem in evolutionary biology [1][2][3] . In theory, sexual populations ought to be replaced by asexual competitors, which can reproduce up to twice as efficiently by investing resources in independently reproductive females 4,5 , versus males that make no contribution to growth. Despite this and other evolutionary advantages of asexuality 6 , biparental meiotic sex is nearly universal among eukaryotes 7,8 . Less than 1% of animal and plant species are obligately asexual, and almost all are of recent origin, implying that sex brings such major benefits that clonal lineages are driven extinct 9,10 .
Many genetic and ecological hypotheses have been proposed to explain the prevalence of sex 11,12 , but none has decisive support. One prominent idea is that shuffling genes through sex facilitates adaptation to the intense, ubiquitous and relentlessly changing selection imposed by coevolving parasites and pathogens [13][14][15] . This 'Red Queen Hypothesis' (RQH) has broad theoretical appeal as a mechanism to favour genetic mixing and suppress asexuality, either by itself 16,17 or in combination with other processes [18][19][20] , especially when pathogens are highly virulent 21 . It is hard, however, to find appropriate systems to test this hypothesis in nature 22,23 . Cases where sexuals and obligate asexuals coexist offer some support, but these are unusual 6 , and even studies lasting many years do not necessarily reflect stable evolutionary outcomes 24 . Given that asexuals can outcompete sexual relatives locally or temporarily 25 , and many predicted fitness costs of clonality are not easily measured 26 , a longer view is vital for understanding why they almost inevitably disappear.
One approach is to investigate taxa where males and sex seem to have been absent for extended evolutionary timescales 27 . Whatever mechanisms maintain sex in the majority of organisms ought to be absent, weakened or unusually mitigated in these groups. Understanding these cases can thus shed light on the rules that govern sex and keep clonality in check everywhere else 28 . The RQH posits a key role for coevolving parasites and pathogens. If so, 'ancient asexuals' ought to have unusual genetic or ecological mechanisms for dealing with disease, which would be of interest not only in understanding sex, but in helping to test broader models of coevolution, genetic diversity and immunity [29][30][31][32] .
Bdelloid rotifers present a useful system for addressing such questions. This class of tiny freshwater invertebrates has persisted for tens of millions of years and diversified into over 450 species 33,34 . Reproduction is only known by parthenogenetic eggs and neither males nor sperm have been reported in contemporary populations, unlike other putatively ancient asexual animals [35][36][37][38][39][40] . Genetic evidence to either confirm or refute their long-term asexuality has proved complicated to obtain, and evidence previously reported for and against this hypothesis [41][42][43][44] has been overturned by later studies [45][46][47][48][49] . Within individuals, genes occur in multiple copies that appear to share sequence homology via ongoing gene conversion 43,45 and recombination 50 , but some authors have also proposed unusual modes of genetic exchange between individuals, based on exotic mechanisms unknown in other animals 44,51,52 . While it is difficult to exclude that bdelloids engage in some form of inter-individual recombination that remains obscure 50 , they nevertheless comprise an exceptional animal clade in which conventional sex is so rare or cryptic as to elude centuries of observation 53,54 and remain opaque to direct study.
If the RQH helps explain the prevalence of sex, we might expect that bdelloids either lack virulent pathogens, or possess unusual characteristics that mitigate the impact of these natural enemies. In fact, bdelloids are attacked by over 60 species of virulent fungal (e.g. Fig.  1a) and oomycete pathogens that specialize on them 55 . These can exterminate cultured populations in a few weeks 56,57 , and significantly depress the abundance of hosts in natural habitats 58 . Controlled infection trials show substantial variation in resistance and infectivity phenotypes among host and fungal pathogen isolates ( Fig. 1b; unpublished data), suggesting the possibility of genetically specific mechanisms of immunity whose efficacy varies depending on the host-pathogen pairing, as envisaged by Red Queen models 16,59 and observed in other invertebrate systems [60][61][62][63][64][65] . Nothing is currently known, however, about the genetic basis of these resistance traits in bdelloids (or rotifers more generally), or how they might evolve and remain effective in the absence of sex.
One potential genetic mechanism that could help bdelloid rotifers diversify and sustain resistance against pathogens is the incorporation of new genes by horizontal gene transfer (HGT). Approximately 10% of genes in bdelloid genomes appear to have been acquired from bacteria, fungi, plants and other non-metazoan sources, rather than sharing recent common ancestry with orthologs in related animal groups 43,51,66 . This estimate is an order of magnitude greater than in other animals 67 , holds for all bdelloid genomes so far examined, and is robust across various approaches to identifying HGT events 46,[68][69][70] .
Comparisons among bdelloid species indicate that most HGT events are ancient but with ongoing rates of stable gene acquisition estimated to be on the order of one gain per 100,000 years 71 . At these rates, the phenomenon is too slow to be equated with the shuffling of alleles in a population by sex, but HGT has been hypothesized to introduce novel biochemical functions that help bdelloids adapt to environmental challenges. Foreign genes are indeed expressed and incorporated into metabolic pathways 66,68 , some of which are not shared by other metazoans 72 . Furthermore, these genes have been found to play a role in desiccation tolerance, nutrient exploitation and repair of DNA damage 66,68,73 . However, the intense and continual selection resulting from evolving pathogens might particularly favour the acquisition and co-option of genes involved in immune responses, such as defensive proteins or antimicrobial compounds. HGT in bacteria is commonly associated with functions and products that help combat antagonists 74,75 and isolated examples of HGT contributing to immunity are known from invertebrates [76][77][78] . The massive scale of HGT in bdelloids might therefore have brought new gene functions that could compensate in part for the theoretical challenge that pathogens pose to asexual groups, as posited by the RQH.
To test this prediction, we used RNA-seq to identify genes that are differentially expressed when bdelloid rotifers are attacked by a fungal pathogen that in the genus Rotiferophthora (Clavicipitaceae, Hypocreales), which preys specifically upon them 79 . We assessed variation in this response by comparing two species, Adineta ricciae and A. vaga, which differ in resistance by a factor of four to the same pathogen. We used functional annotation to identify putative defensive or immune-related genes and tested whether these showed evidence of horizontally acquisition from non-metazoan taxa. Finally, we compared our results with RNA-seq data obtained when bdelloids were exposed to desiccation 73 , enabling us to test whether genes of foreign origin contribute disproportionately when responding to a biotic as opposed to an abiotic stressor. We found strong support for the prediction that HGT plays a disproportionate role in defence against pathogens. Hundreds of foreign genes were differentially regulated in response to fungal attack, some producing proteins not previously implicated in immune defence in animals. Our results demonstrate the scale and range of functions that HGT adds to the genetic repertoire of bdelloids, particularly in the important context of adaptation to biotic antagonists.

Results and Discussion
Clonal populations of the bdelloid rotifers Adineta ricciae and A. vaga were experimentally inoculated under controlled conditions with spores of the virulent fungal pathogen Rotiferophthora globospora 79 , as part of a larger study. In both species, 95% of animals ingested spores and contracted (Fig. 1a) within minutes of exposure; they were then monitored over 3-4 days (see Methods). The proportion of animals succumbing to fungal infection differed greatly between the species: 18% for A. ricciae and 71% for A. vaga (Fig. 1b, relative risk = 3.74, 95% CI: 2.8-5.0, z = 8.758, P <0.001). Given that these closely related species are ecologically 46,80 and morphologically similar (Fig. 1a) and were reared under standard conditions, the basis of the difference in mortality seemed to be physiological, reflecting a strong defensive response in A. ricciae, and a weaker defence in A. vaga. To investigate the genetic basis of this response, we repeated the inoculation experiment at scale, and harvested total RNA from exposed and control populations of each species in triplicate at 7h and 24h, during the early stages of fungal attack (Methods). We took advantage of previously published reference genomes for both species 43,46 to assemble, map and quantify sequenced transcripts, and compared species and timepoints for patterns of differential gene expression between inoculated and control populations (Fig. 1c, Methods). ; composite image of three A. ricciae corpses with fully developed R. globospora infections, with hyphae differentiating into irregular resting spores and long conidiophores bearing spherical infectious conidia (scale bars 100µm). b. Proportion of rotifers killed by infection 72 hours after initial exposure to R. globospora. Points indicate replicate laboratory populations with approximately 10 individuals (range 6 to 20) exposed to 1000 conidia. c. Numbers of genes showing significant upregulation (solid bars) and downregulation (hatched bars) at 7-and 24-hours post-inoculation (T7 and T24 respectively), for A. vaga (red) and A. ricciae (blue), relative to control populations inoculated with UV-inactivated spores.
We defined significant differential expression stringently, as a >4-fold difference in expression level and a False Discovery Rate-adjusted P-value < 0.001 (Fig. 1, Methods). By these criteria, 541 A. vaga genes and 962 A. ricciae genes were differentially expressed at 7 hours post-inoculation (T7). Most of these were upregulated in the pathogen treatment group relative to controls (452 in A. vaga, ~84%; 709 in A. ricciae, ~74%). At 24-hours postinoculation (T24), the number of differentially expressed genes rose to 1767 in A. vaga (with 1093 upregulated, ~62%), and 2388 in A. ricciae (1590 upregulated, ~67%) (Fig. 1). The identities of the differentially regulated genes in each species at T7 overlapped substantially with those at T24, and the direction and magnitude of differential expression at both timepoints was significantly correlated (Pearson's correlation R = 0.66-0.84; P < 2e-16 in all cases; Supplementary Fig. 1). This indicates that many of the transcriptional changes initiated at T7 were ongoing and consistent at T24, even as new genes joined the response. In both species, A. ricciae A. vaga R. globospora the numbers of differentially expressed genes showed a comparable increase between T7 and T24, but the response in A. ricciae was initially stronger and remained higher, especially among upregulated genes. The more rapid and extensive transcriptional response in A. ricciae might be a broad-scale physiological reflection of its ability to recognise the pathogen sooner and mobilise defences more extensively than A. vaga. The two species showed overlapping profiles of differentially expressed genes, as expected if some defensive or recognition pathways were conserved (Supplementary Fig. 2; Supplementary Table 3). For example, of the 1093 genes that were significantly upregulated in A. vaga at T24, 552 (50.5%) shared an ortholog that was also significant upregulated in A. ricciae, and these orthologs showed a significant correlation in the magnitude of differential expression (Pearson's correlation R = 0.62, P < 2e-16). Based on the predicted proteins from the reference genomes 46,70 , the background frequency of candidate HGT genes (HGTC) in both species is ~11%. Among the significantly upregulated genes, however, the frequency of HGTC increased to 20-30% ( Fig. 2; Table 1), which represents a 2-to 3-fold enrichment (Fisher's exact test for an association between upregulation and HGTC, P < 0.001 in all cases). The pattern was repeated at both timepoints; the proportional enrichment was stronger in the early response, 7 hours after inoculation, but the absolute numbers of upregulated HGTC genes were higher at T24, with 345 genes of putatively foreign origin upregulated in A. ricciae, and 269 in A. vaga. The proportion of HGTC was not significantly different between the two species ( Table 1). These results are robust to the threshold of fold-change used to define differential expression (1.5-fold, 2-fold, 8-fold, and 16-fold absolute differences yielded the same results as 4-fold, Supplementary Figs. 3, 4). We also observe a significant but lesser enrichment of HGTC among genes that are significantly downregulated during the experiment (Fisher's exact test, P < 0.001 in all cases; Fig. 2; Table  1; Supplementary Figs. 3,4). Horizontally acquired genes were therefore significantly overrepresented among the genes that are differentially expressed, and especially upregulated, in the infected treatments. Comparison of control populations (Supplementary Figs. 5,6) showed that HGTC genes are more likely to be expressed than native genes under control conditions, but at significantly lower levels. The enrichment of HGTC among differentially expressed gene is thus unlikely to be explained by a known false-positive bias in some RNA-seq analyses toward genes with higher expression levels 81,82 . <0.001 1 "Pathogen": analysis of DE where treatment represents animals exposed to live pathogen spores, versus controls exposed to inactivated spores, at 7h (T7) and 24h (T24) post exposure. "Desiccation": analysis of DE where treatment represents animals either entering into or recovering from desiccation, versus control animals that remained hydrated (Hecox-Lea & Mark Welch 2018). 2 "Up": genes that are upregulated in the treatment group; "Down", genes that are downregulated in the treatment group. 3 "%HGTC in DE0": the proportion of HGTC among genes that are not significantly differentially expressed, either up or down depending on DE type; "%HGTC in DE1": the proportion of HGTC among genes which are significantly differentially expressed. Significant DE is defined as genes with >4-fold absolute change in expression and FDR < 1e-3. 4 Fisher's exact test for a significant association between classifications; null hypothesis: true odds ratio not equal to 1.

Fig. 2. Up-and downregulation of native and foreign genes in response to fungal infection in two species of bdelloid rotifer.
Data are shown for a. A. vaga T7; b. A. vaga T24; c. A. ricciae T7 and d. A. ricciae T24. Each point on the 'volcano' plot represents a gene plotted by log2 fold-change in expression level on the X-axis, and significance level (expressed as -log10 FDR) on the Y-axis. Positive X-axis values indicate genes that were upregulated in the treatment groups (i.e., those exposed to the live pathogen) relative to the control groups, whereas negative values indicate genes that were downregulated in the treatment group relative to the control group. Genes with significant changes in expression level (absolute fold-change > 4 and FDR < 1e-3, dashed lines) are shown in colours (red: A. vaga; blue: A. ricciae); non-significant changes in expression are indicated with grey (note that the majority of genes are non-significant at these thresholds; see legends for details). P-values in blue show the probability of non-association between HGTC and the corresponding up-or downregulated subset (Fisher exact tests, see Table 1 for details).
We performed further analyses to investigate whether enrichment of HGTC genes is specific to the response to fungal infection, or a more widespread phenomenon among differentially expressed bdelloid genes. For example, HGTC genes might be more likely to be differentially expressed in general if they are overrepresented among effectors at the peripheral ends of gene regulatory networks or if their expression is less tightly regulated than that of typical host genes 83,84 , or systematically higher, leading to read-count bias 81 . Finally, HGTC genes might be overrepresented among general stress-response genes, rather than associated specifically with the response to a pathogen; for instance, evidence from earlier analyses has suggested a role for such genes in DNA repair, detoxification and desiccation tolerance in bdelloids 66,68,73 .
We therefore repeated our analyses on published RNA-seq data from a recent study 73 that investigated gene expression in the same cultured strain of A. vaga responding to desiccation. Most bdelloid species are remarkably tolerant to drying out, entering a state of physiological dormancy known as anhydrobiosis 85,86 . Hecox-Lea and Mark Welch (2018) compared levels of gene expression in hydrated animals to those that were either entering desiccation or recovering from a period of dormancy induced by desiccation. Thus, we were able to test whether HGTC genes are also enriched in the response to desiccation, using the same criteria applied above for differential expression during infection. Overall, the proportion of HGTC genes involved in the desiccation response was substantially lower than that seen in the pathogen experiment (~9-19%, depending on the subset, compared to ~22-30% previously; Table 1). There was a significant enrichment of HGTC when animals were recovering from desiccation (Fisher's exact test, P < 0.001; Table 1; Supplementary Fig. 7), but not when animals were entering desiccation. Thus, HGT candidates appear to be disproportionately represented among genes responding strongly to both kinds of stresses, but the effect is ~2-fold stronger during the response to pathogens compared with desiccation, consistent with a more prominent role for HGTC in immune defence.
Next, we compared the identities of genes that were differentially expressed in the fungal infection and desiccation experiments, to test how far these responses are specific to each challenge, rather than reflecting general stress response pathways in A. vaga common to both treatments. In general, the proportion of upregulated genes shared between any two sample sets was low (mean = ~26%; Supplementary Table 4), with the exceptions of T7 versus T24 in the infection experiment (~90% of T7 genes also were upregulated at T24), and entering versus recovering from desiccation (~63% of 'entering' genes also were upregulated in 'recovering'). HGTC genes showed a similar pattern. For example, of the 256 HGT candidates upregulated during recovery from desiccation, only 31 (~12%, T7) and 61 (~24%, T24) were also upregulated during pathogen response (Supplementary Table 4). Even fewer genes were shared between downregulated subsets (none and 24 for T7 and T24, respectively). The majority of differentially expressed HGTC in A. vaga therefore are unique to either the desiccation or infection experiment.
Having identified hundreds of differentially expressed HGTC genes that appear to be associated specifically with infection, we investigated putative functions by annotating the translated transcriptomes of A. vaga and A. ricciae based on the presence of Pfam and InterPro domains (see Methods), and found several proteins of particular interest. In A. vaga, at T7 and T24, the two foreign genes with the strongest evidence of upregulation are homologous copies of each other, encoding putative RNA ligase (Pfam accession PF13563) and endo/exonuclease domains (PF03372). A further 16 HGTC genes with putative RNA ligase activity (including PF09414 and PF09511) were significantly upregulated across either timepoint. A similar pattern was seen in A. ricciae, and gene ontology (GO) analysis confirmed enrichment of RNA ligase/repair activity for both species. Phylogenetic analysis of rotifer-encoded RNA ligase domains supported the classification of these genes as HGTC, showing high support for the clustering of rotifer copies with non-metazoan lineages (including bacteria, fungi, viruses, and non-metazoan eukaryotes; Fig. 3a-c). The hypothesis that these sequences arise from contamination is refuted by the presence of corresponding paired copies on genomic scaffolds for A. vaga and A. ricciae whose divergence matches the rest of the genomes, by the fact that these genes are quite distinct from the nearest sequenced non-metazoan copies, and by the presence of introns in genes encoding proteins whose closest homology is otherwise to bacteria. Many upregulated gene copies were closely related to each other and showed expected patterns of homologous, homoeologous or orthologous relationships given the ancestral tetraploidy in bdelloid species 45,87 .
While the precise function of these genes is unknown, a role for RNA ligase activity in the response of bdelloids to fungal infection has intriguing links to work in other systems. Two close relatives of Rotiferopthora in the fungal order Hypocreales attack their insect hosts using lethal ribotoxins as virulence factors [88][89][90][91] . These secreted ribonucleases cleave the critical and highly conserved sarcin-ricin loop in the large ribosomal subunit of host cells, which inhibits protein synthesis and causes death 89 . If Rotiferopthora hyphae secrete similar virulence factors against Adineta, rapid upregulation of proteins with RNA ligase activity might help the host counter this element of the attack by protecting and repairing ribosomes 92 . According to a recent review, RNA toxins, RNA ligases and related RNA repair systems are "extensively disseminated by lateral transfer between distant prokaryotic and microbial eukaryotic lineages consistent with intense inter-organismal conflict" 92 . Might bdelloid rotifers have borrowed machinery from this microbial RNA warfare to defend against their own pathogens? Future work will be needed to explore this speculative hypothesis, by testing whether R. globospora expresses such a toxin and determining its effects on rotifer hosts.
Another potential source of immunity to fungi is production of carbohydrate-active enzymes (CAZymes) to target components of fungal cell walls such as chitin and glucans, either acting as recognition receptors or degrading them directly [93][94][95] . We surveyed 35 glycosyl hydrolase and carbohydrate binding gene families for significant up-or down-regulation, but found no large-scale upregulation of putative CAZymes in response to fungal infection in bdelloids, including genes with putative chitinase and chitin-binding properties. However, a small number of HGTC genes with putative glucanase activity were upregulated in one or both species. These included six genes in both species containing GH16 domains (active on β-1,4 or β-1,3 glycosidic bonds in various glucans and galactans) with closest matches to the bacterial phylum Bacteroidetes, and two genes restricted to A. vaga with glycosyl hydrolase family 64 (GH64, a β-1,3-glucanase) domains that cluster with homologs from firmicute bacteria (Fig.  3d, e). Glucanases are secreted by plants and bacteria as antifungal compounds 96,97 , and related domains are implicated in insect antifungal defences against Metarhizium, a close relative of Rotiferophthora 94 . Future work could explore whether these acquired genes play a similar role in bdelloids.
We searched the upregulated native and foreign genes for candidates with known roles in immunity in other invertebrate models 60,98,99 , such as Toll-like receptor genes (TLRs), caspases, thioester proteins, scavenger receptors, Gram-negative bacterial proteins (GNBPs) and peptidoglycan recognition proteins (PGRPs). We found little evidence for major upregulation of these gene families during exposure to the pathogen, though many of these pathways are detectable in bdelloid genomes. For example, for TLRs, there are 120 and 133 proteins with significant matches to the Toll/interleukin-1 receptor-like (TIR) domains (PF01582 and PF13676) in A. vaga and A. ricciae, respectively. However, we found only three putative A. vaga TLR genes to be significantly DE on exposure to the pathogen, and these were downregulated at T24; nine were found in A. ricciae, all but one downregulated. These observations do not exclude a role in immunity, since TLRs act as signalling components in other invertebrates 99,100 , and may be expressed constitutively rather than upregulated in response to pathogens. For caspases, we detected 101 and 70 proteins with a significant match to the Peptidase_C14 (PF00656) domain in A. vaga and A. ricciae respectively. Many of these were marked as HGTC from bacteria (i.e., putative metacaspases) but in fact are highly divergent from any known caspase homologs found in UniProt (Fig. 3f), and their origin and evolution remain obscure. Similar to TLRs, only a relatively small number of putative caspases were found to be significantly differentially expressed on exposure to the pathogen, and the majority were downregulated. One possibility is that some genes with a role in immunity were downregulated in our experiment because they respond to non-fungal pathogens, such as bacteria. In nematodes, induction of antifungal defences is correlated with selective repression of antibacterial immune response genes 101 , perhaps to focus resources on an immediate threat or balance biochemical trade-offs in defence mechanisms 102 . Functional inferences about these genes are complicated by the substantial phylogenetic distance of rotifers from wellcharacterised models of invertebrate immunity 60,103 , as well as putative horizontal origins in some cases.  (PF13563)

Conclusions
We find that genes acquired by bdelloid rotifers from non-metazoan sources play a disproportionate role in a specific transcriptional response to infection by a fungal pathogen. The speed and scale of this response was greater in the more resistant of the two host species we compared, though expression differences were subtle relative to a fourfold difference in susceptibility. The enrichment of foreign genes was twice as pronounced in the response to pathogens as compared with desiccation, suggesting that their expression patterns are specific to biotic defence rather than a general reaction to stress or physical damage. Since virtually nothing is currently known about the genetics of immunity in rotifers, and tools for assessing the functions of candidate genes in bdelloids are lacking, the precise roles of the differentially regulated genes remain to be determined. Nonetheless, we hypothesise that RNA ligase and glucanase functions acquired from bacteria and other sources have been co-opted for immune defence against fungi in particular. While some putative horizontally acquired genes might in fact represent novel gene families with a deep or obscure metazoan heritage, the horizontal origin of the genes we highlight is well supported by phylogenetic evidence grouping them closely with non-metazoan orthologs.
Overall, our findings are consistent with the hypothesis that horizontal gene transfer contributes to bdelloid defence against pathogens. This could be interpreted to support the view that diseases pose a particular challenge for asexual lineages 13,14 , with special measures required for longstanding parthenogenetic lineages to keep up 9,27,28 . But to what extent would this process help bdelloids overcome the presumed deficiencies of long-term asexuality in combatting disease? Many of the genes that were upregulated in infected populations were acquired prior to the common ancestor of A. vaga and A. ricciae, and some appear to be shared even more deeply among bdelloid rotifers. Such genes could provide an enhanced biochemical repertoire for combatting pathogens overall, potentially including functions unavailable to other animals, but they would not create variation in bdelloid populations at the speed or scale theoretically required to sustain rapid coevolution 15,104-108 . Nonetheless, a handful of upregulated foreign genes were restricted to just one Adineta species and therefore interpreted as more recently acquired. It is possible that new genes do arrive at a fast enough rate to diversify immune function at the scale of lineages and species, and thereby provide some variation to counter long-term extinction, especially if ecological factors 56,58 help susceptible rotifer clones avoid coevolving pathogens in the shorter term.

Rotifer and pathogen isolates
Animals belonging to the species Adineta ricciae 80 and A. vaga 109,110 were isolated respectively in 1998 from mud in Australia 80 and ca. 1984 from moss in Italy 111 . These have been propagated clonally in continuous long-term culture across several laboratories 43,48,112,113 and annotated genome assemblies are available for both lines 43,46 . Our cultures were maintained in 60mm plastic Petri dishes in sterilised distilled water, fed with Escherichia coli (OP50) and Saccharomyces cerevisiae (S288c) and subcultured approximately once per month. They were stored at 20°C in an illuminated incubator (LMS, Kent, UK) with a 12:12 hour light:dark cycle.
The fungal pathogen Rotiferopthora globospora 79 was found attacking co-occurring rotifers of the genus Adineta in soil in northern New York 57 . A pure culture on potato dextrose agar (PDA) was obtained in December 2008 using methods described elsewhere 114 , and deposited in April 2009 with the USDA Agricultural Research Service Collection of Entomopathogenic Fungal Cultures (ARSEF) for long-term cryogenic storage under the accession code ARSEF 8995. Frozen mycelium was retrieved and revived from this collection in 2013, and since maintained in serial subculture on PDA at 20°C with a 12:12 hour light:dark cycle and with transfers every 4 months. This pathogen strain has no recent history of laboratory co-passage with either of the two rotifer hosts tested here and all three were isolated on different continents. However, R. globospora appears to be globally distributed-it was originally described from New Zealand 79 and has also been recorded in Japan (CGW, pers. obs.), in both cases attacking Adineta. Therefore, either or both host species could plausibly have encountered this pathogen in nature, especially given the high dispersal capacity and global distribution of bdelloid rotifers 58,[115][116][117] .
Infection by R. globospora is initiated when a host ingests infectious spores (conidia), which are spherical and approximately 3.5µm in diameter 79 . These lodge in the mouth or oesophagus, at which point the animal stops feeding and contracts within a few minutes (Fig.  1a). Over the next 6-8 hours, the conidium produces a thin germ tube that penetrates the gut wall and swells into assimilative hyphae, which begin to invade and digest the surrounding host tissue after about 12-24 hours, as the infection becomes established. By 36-48 hours, fungal hyphae have filled most of the body cavity and the host is dead. Between 48-72 hours, hyphae begin to emerge through the host integument, and will eventually differentiate to form two types of spores: a handful of thick-walled resting spores, and hundreds of fresh infectious conidia (Fig. 1a).
Because Rotiferophthora can only complete its life cycle by killing its host, these fungi are inherently highly virulent pathogens, and have been described as "devastating" to populations of Adineta 79 . R. globospora (ARSEF 8995) has been shown to exterminate laboratory populations of a sympatric Adineta c.f. vaga clone within 28 days of initial exposure to a low density of conidia 57 . However, even if rotifer individuals are inevitably killed once the pathogen becomes established, the possibility remains that an individual host can resist the initial attack and prevent an ingested spore from establishing a successful infection, or delay its progression. Even partial resistance at an early stage could dramatically slow the spread of an epidemic 118 , giving clonal relatives time to escape the infested habitat 56,58 before the local population is exterminated. Early-acting mechanisms of resistance to Rotiferophthora ought therefore to be favoured in bdelloids, corresponding to the innate immune pathways identified in other invertebrates facing virulent fungal pathogens 101,119-121 .

Infection assays
To quantify resistance to R. globospora under standard conditions, rotifers were transferred by pipette from stock populations to a 2mL droplet of sterile distilled water, where eggs, corpses and food from the cultures were washed away. Adult individuals were transferred to 96-well plates (Thermo-Fisher), with approximately 11 animals per well (mean: 11.0, SD: 3.5) in 60µL of sterilised, distilled water. Rotifers were counted using a compound microscope (Nikon Eclipse E400), noting whether each animal was active (feeding or locomoting, Fig. 1a), contracted (with head withdrawn) or dead. Animals that died during the transfer (<1.5% of the total) were physically removed where possible, or recorded so they could be excluded from later counts.
To obtain pure suspensions of conidia, approximately 5x5mm of freshly subcultured sporulating mycelium on PDA was moved to 500µL of sterile distilled water in a 1.5mL Eppendorf tube. After vortexing, 400µL of conidial suspension was transferred to a fresh sterile tube and conidial density was measured using a haemocytometer and lactophenol cotton blue stain. An inactivated inoculum was prepared simultaneously as a control, by exposing an aliquot of conidia in water to a germicidal ultraviolet lamp (25W, 253.7nm) at a distance of 10cm for 45 minutes, with regular vortexing to ensure all spores were irradiated. This treatment was successful, as none of the rotifers in control wells treated with irradiated spores became infected.
Wells were inoculated with 8µL of freshly prepared conidial suspension at a density of 125 spores µL -1 . Negative control wells received 8µL of distilled water or inactivated spore suspension. The final density of spores in each well was high (ca. 15 conidia µL -1 ) to ensure every animal was exposed to the pathogen in a synchronised pulse. This appeared to work, because 94% of animals in experimental wells were contracted 8 hours after exposure, versus a baseline of 3% of animals contracted in control wells. Plates were then stored in incubators (LMS 300NP) at 20 o C, with a 12:12 hour light:dark cycle.
After 48 and 72 hours, rotifers were counted again, and classed as active, contracted, killed by infection or otherwise dead. A rotifer was considered to be killed by infection if at least one hypha had emerged through the integument from the interior. This criterion was unambiguous, in contrast with the difficulty of determining whether a fungal infection had established inside a contracted animal. By 72h, the proportion of infected animals in experimental wells appeared to have stabilised. In a subset of wells recounted at 96h, only a small fraction of animals (~6.5%) had newly developed visible infections. The relative risk of infection for A. vaga versus A. ricciae at 96h (RR 2.94, 95% CI 1.93-4.46; Z = 5.06, P < 0.0001) had not narrowed significantly since 72h (ratio of RR 0.78, 95% CI 0.47-1.31, Z = -0.926, P = 0.354), but the survivors had reproduced, making further tracking of the originally exposed cohort difficult. We therefore took the 72h timepoint as the standard measure of infection mortality. Individuals in negative control wells never became infected, whether they received water or sterilized spores, and the background death rate was ~2%.
Inoculation trials with both A. ricciae and A. vaga were replicated on multiple occasions using at least three separate source dishes for each species prepared at different times (in some cases different years), and separately prepared cultures and inocula of fungi, with each set of trials replicated in at least three wells with accompanying negative controls in appropriately blocked designs. The total number of animals exposed across all trials was 216 for A. ricciae and 189 for A. vaga (Fig. 1b). Rates of mortality from R. globospora infection were consistent across these trials, as was the dramatic difference in susceptibility between the species (Fig.  1b). One set of well trials was run simultaneously and in close linkage with the RNA-seq experiment described below, using leftover rotifers from the same source populations, and inoculating wells with the same live and irradiated pathogen suspensions at the same time. This enabled us to infer the timings and presumed final outcomes of infections in the RNA tubes, even though those animals could not be directly counted and were sacrificed for RNA by 24h.
The results of the RNA-linked well experiments were similar to earlier trials: mean infection mortality at 72h was 18% for A. ricciae and 79% for A. vaga.

RNA-seq experimental design
Rotifers for the RNA-seq experiment were reared in eight replicate Petri dishes per species, with ca. 50 founders per dish, all from the same clonal laboratory line of each species. These were fed only with E. coli (OP50, 5 x 10 8 cells per dish) in distilled water. S. cerevisiae was omitted to avoid introducing further eukaryotic transcripts, and because gene expression by rotifers metabolising a fungal food source might complicate inferences about transcription in response to a fungal pathogen. The RNA-linked well assays described above yielded similar results to earlier trials, so the omission of S. cerevisiae as a food did not affect infection outcomes. When the A. ricciae dishes were initially established, each was inoculated with 50µL of 5µm-filtered water from the A. vaga source population, and vice-versa, to homogenise any co-cultured bacterial communities whose composition might differ between the cultured lines of the two species, and limit this as a potential source of gene expression differences. Dishes of the two species were stored in evenly interspersed blocks while the populations were growing.
Rotifers were counted and harvested after 4 weeks, when the mean population size was 2965 per dish for A. ricciae and 1688 for A. vaga, which reproduces more slowly 85 . Each dish was gently washed, using several water changes to float and pour off eggs, bacterial cells and corpses, while active animals adhered to the plastic. Cleaned rotifers were detached from the plastic by dropping distilled water onto them from a 10mL syringe at 40cm height, followed by repeated aspiration and forcible expulsion of medium using a P1000 pipetter and a 1mL tip. Physical detachment was preferred over a salt or cold shock approach 43 , to reduce the chance of inducing transcriptomic, immunological or behavioural changes. Suspended rotifers were poured into a 50mL centrifuge tube and pelleted at 3000 x g for 5 minutes using a swing-bucket cradle. All but 1.5mL of supernatant was removed, then a vortex shaker was used to resuspend the rotifer pellets and 1mL was transferred immediately by pipette to a 1.5mL Eppendorf tube. This process was repeated to yield 16 replicate 1.5mL tubes for each species. Each tube was centrifuged at 17000 x g for 3 minutes, and all but 100µL of water was removed from the rotifer pellet. The pelleted rotifers were left overnight to recover and redistribute themselves around the submerged interior of the tube. We estimate that the efficiency of rotifer recovery via this method is about 70%, based on numbers of animals left over in plates or tubes, giving approximately 1000 animals per tube for A. ricciae and 600 for A. vaga.
The 16 tubes for each species were randomly allocated to receive either live or irradiated pathogen spores, and to have RNA extracted either 7 or 24 hours later, with each combination replicated four times. Irradiated spore suspensions, as described above, were used as a control treatment to account for physical, chemical, or nutritional effects of ingesting fungal cells, so that all else was equal except for pathogenic activity. RNA sampling times were chosen to correspond to the early stages of infection 114 -by 7h, germ tubes from ingested spores have attempted to infiltrate the host, but assimilative hyphae would not yet have become established. By 24h, assimilative hyphae would be present if the germ tube had been successful, but these would not yet have colonised the host extensively or killed it. Each tube was inoculated with 20µL of live or irradiated spore suspension at a density of 500 spores µL -1 , for a total of 10,000 spores and a final density of 80 conidia µL -1 , to ensure every animal was exposed as synchronously as possible. Tubes were incubated upright at 20 o C in a blocked layout until RNA extraction.

RNA extraction and sequencing
Total RNA was extracted from each tube at the appropriate timepoint using an RNeasy Mini kit (Qiagen #74104), following the manufacturer's protocol for animal tissues. To lyse and homogenise the rotifers, tubes were centrifuged at 17,000 x g for 3 minutes and 100µL of excess water was removed, leaving rotifers pelleted in approximately 20µL of water. After adding 50µL of Buffer RLT, the pellet was immediately disrupted and homogenised for 1 minute by pulsed application of a pellet pestle attached to a cordless motor (Kimble). A further 250µL of Buffer RLT was added to stabilise the lysate and rinse residue from the pestle, before proceeding with the manufacturer's protocol, including an on-column DNase I digestion step. Lysis and stabilisation were completed for all tubes within 30 minutes of the target timepoint, in a balanced order with respect to treatment group and species.
RNA was eluted in 32µL of RNase-free water and 1.5µL aliquots were analysed using a Nanodrop 2000 (ThermoFisher). Spectrophotometric measurements were used to select the three replicates with the highest RNA concentrations from each treatment group for further analysis and sequencing, so that downstream steps had three biological replicates. These 24 tubes were frozen at -80 o C and shipped on dry ice to the University of Edinburgh (Edinburgh, UK), where further quality control was undertaken by Edinburgh Genomics, including RNA quantitation with a Qubit 2.0 fluorometer in duplicate for each sample, and RNA ScreenTape analysis with an Agilent 2200 TapeStation system to assess RNA integrity. All 24 samples proceeded to cDNA library preparation, using the TruSeq stranded mRNA kit (Illumina) to enrich for polyadenylated transcripts. The indexed libraries were sequenced in multiplex on an Illumina NovaSeq 6000 at Edinburgh Genomics, using an SP flow cell to generate 50-base paired end read libraries with 200 bp inserts.
A total of 102.9 Gb of raw sequencing data were generated over the 24 RNA libraries. After data filtering and quality-control, 78.5 million reads were retained per library (94.2 Gb total data, Supplementary Table 1). Over 99% of filtered reads mapped to the A. vaga and A. ricciae reference genomes (Supplementary Table 2). Filtered data were de novo assembled using the Trinity software v2.8.9 122,123 , resulting in 118,860 and 63,749 transcripts representing 41,527 and 35,079 'genes' for A. ricciae and A. vaga, respectively (Supplementary Table 2). Transcriptomes showed G + C proportions that were consistent with their respective genomes (A. ricciae 37.6%; A. vaga 33.0%) and high 'completeness' scores as measured by BUSCO analysis (>97% of core eukaryote genes completely recovered in both cases). Transcript-togenome mapping rates were >99% in both species, with 87.1% (A. ricciae) and 81.8% (A. vaga) of introns correctly called based on previously published gene models 46 .

Data filtering and quality control
Raw RNA-seq sequencing reads were quality and adapter trimmed using BBTools 'bbduk' v38.73 with the parameters 'ktrim=r k=23 mink=11 hdist=1 tpe tbo', and error-corrected using BBTools 'tadpole' with the parameters 'mode=correct tossjunk tossuncorrectable' (https://sourceforge.net/projects/bbmap/). Unwanted reads derived from ribosomal RNA (rRNA) were removed by mapping the data to the SILVA rRNA database 124 using BBTools 'bbmap' with the parameters 'local=t outu=filtered_R#.fq.gz' (i.e. retaining only unmapped pairs of reads). Sequences derived from fungal contaminants were also removed using the same approach, mapping to sequenced genomes of fungi in the family Clavicipitaceae (NCBI:txid34397). The proportion of reads mapping to the A. ricciae or A. vaga target genome was assessed using STAR v2.7.3a 125 with the parameter '--twoPassMode Basic'. Final data quality was assessed visually using FastQC 126 and MultiQC 127 . All raw sequencing data have been deposited in the relevant International Nucleotide Sequence Database Collaboration (INSDC) database with the Study ID PRJEB39927 (see Supplementary Table 1 for run  accessions).

Differential expression analysis
Transcript quantification was performed using Salmon 'quant' v0.14.1 128 , using the gene annotations of Nowell et al. 46 as the target transcriptomes. Short transcripts <150 bases were removed prior to analysis, and genomic scaffolds were appended to each transcriptome data as 'decoys' prior to quantification as recommended in the Salmon documentation. Relationships between biological replicates within and between samples were checked visually using utility scripts from the Trinity software 122,123 . Statistical analysis of the resulting count matrix was performed with DESeq2 129 , which uses negative binomial generalized linear models to test for differential expression. P-values were adjusted for multiple testing using the Benjamini-Hochberg method 130 to control the false discovery rate (FDR). Stringent thresholds of FDR < 1e-3 and log2 fold-change > 2 (i.e. 4-fold difference in expression) were used to define an initial set of differentially expressed genes for downstream analysis.
To assess the quality of the A. vaga and A. ricciae gene predictions 46 , a transcriptome was assembled de novo from the filtered RNA-seq data using Trinity v2.8.9 122 with default parameters, and mapped to the genomes of A. ricciae 46 and A. vaga 43 using Minimap2 131 v2.17 with the parameters '-ax splice -C5'. The degree of correspondence between the predicted gene models and the assembled transcriptome was then assessed using 'paftools.js junceval', which provides a count of the number of correctly predicted introns based on the exact agreement between the location of the intron-exon junction in the annotation file and gaps in the transcriptome alignments. Gene completeness of assembled transcriptomes was also measured by Benchmarking Universal Single-Copy Orthologs (BUSCO) analysis 132 , using the core eukaryote and metazoa gene sets (n = 303 and 978 query genes, respectively).

Horizontal gene transfer detection
Putative HGT candidate genes (HGTC) were determined from the recent analysis of Jaron et al. 70 . Briefly, proteins were aligned to the UniProt90 sequence database 133 using Diamond 'blastp' v0.9.21 140 with the parameters '--sensitive -k 500 -e 1e-10'. For each protein, a HGT 'index' (hU) was computed based on the alignment bitscores to the best-matching hits from the Metazoa (BIN) and non-Metazoa (BOUT) with the formula: hU = BOUT -BIN 68 . A 'consensus hit support' (CHS) score was also calculated as the proportion of secondary hits in agreement with the result based on the hU 46,141 . Proteins were classified as HGTC if hU > 30, CHS > 90% and were physically linked (i.e., found on the same genomic scaffold) to a gene of unambiguous metazoan origin (hU < 0). Phylogenetic trees of selected HGTC and their putative co-orthologs from the UniRef90 database were constructed using IQ-TREE v1.6.12 with the parameters 'alrt 1000 -bb 5000 -m TEST' [142][143][144] .

Gene orthology
Orthologous relationships between A. ricciae and A. vaga genes was determined using OrthoFinder v2.3.12 145 with default parameters. Proteomes from nine other rotifer species 46,47,146 were included in the analysis to aid orthology inference. Explanation of tests: "Timepoints", correlation in differential expression (DE) in the same gene across timepoints T7 and T24 ( Supplementary Fig. 1); "Within genomes", correlation in DE between inferred gene copies within genomes (including putative homologs, homoeologs and other gene copies) ( Supplementary Fig. 2a-d); "Between genomes"; correlation in DE between inferred orthologs across species ( Supplementary Fig. 2e-f). 2 Timepoints: "T7", 7h post-infection; "T24", 24h post-infection. 3 DE type: "Up", genes that are upregulated in the treatment group (i.e., positive DE); "Down", genes that are downregulated in the treatment group (negative DE). Note that all genes are included in correlations regardless of significance. Values in bold on diagonals indicate the total number of genes in category (e.g., 136 HGTC genes upregulated at T7). Percentages in each cell then indicate the proportion of shared genes relative to the row total (e.g., 89.7% of 136 T7 upregulated HGTC genes are shared with T24 upregulated set). Green and red shading is to highlight the relative size of the shared fraction for up-and downregulated subsets, respectively.