
==== Front
Proc Natl Acad Sci U S A
Proc Natl Acad Sci U S A
PNAS
Proceedings of the National Academy of Sciences of the United States of America
0027-8424
1091-6490
National Academy of Sciences

38466855
202307780
10.1073/pnas.2307780121
datasetDatasetresearch-articleResearch ArticleevolutionEvolution418
Biological Sciences
Evolution
Intergenomic signatures of coevolution between Tasmanian devils and an infectious cancer
Gallinson Dylan G. dgallinson@usf.edu
a b 1
Kozakiewicz Christopher P. c d
Rautsaw Rhett M. a c
Beer Marc A. c https://orcid.org/0000-0003-1374-0705

Ruiz-Aravena Manuel e f https://orcid.org/0000-0001-8463-7858

Comte Sebastien e g https://orcid.org/0000-0001-7984-8159

Hamilton David G. e https://orcid.org/0000-0001-5883-0136

Kerlin Douglas H. h
McCallum Hamish I. h https://orcid.org/0000-0002-3493-0412

Hamede Rodrigo e i
Jones Menna E. e
Storfer Andrew c https://orcid.org/0000-0002-9414-8588

McMinds Ryan a b
Margres Mark J. margres@usf.edu
a 1 https://orcid.org/0000-0002-6153-6701

aDepartment of Integrative Biology, University of South Florida, Tampa, FL 33620
bCollege of Public Health, University of South Florida, Tampa, FL 33620
cSchool of Biological Sciences, Washington State University, Pullman, WA 99163
dW.K. Kellogg Biological Station, Department of Integrative Biology, Michigan State University, Hickory Corners, MI 49060
eSchool of Natural Sciences, University of Tasmania, Hobart, TAS 7001, Australia
fDepartment of Public and Ecosystem Health, Cornell University, Ithaca, NY 14853
gNew South Wales Department of Primary Industries, Vertebrate Pest Research Unit, Orange, NSW 2800, Australia
hCentre for Planetary Health and Food Security, Griffith University, Nathan, QLD 4111, Australia
iCANECEV Centre de Recherches Ecologiques et Evolutives sur le Cancer, Montpellier 34394, France
1To whom correspondence may be addressed. Email: dgallinson@usf.edu or margres@usf.edu.
Edited by Stephen Goff, Columbia University Irving Medical Center, New York, NY; received June 7, 2023; accepted January 17, 2024

11 3 2024
19 3 2024
11 9 2024
121 12 e230778012107 6 2023
17 1 2024
Copyright © 2024 the Author(s). Published by PNAS.
2024
https://creativecommons.org/licenses/by-nc-nd/4.0/ This article is distributed under Creative Commons Attribution-NonCommercial-NoDerivatives License 4.0 (CC BY-NC-ND).

Significance

Since the discovery of devil facial tumor disease (DFTD) nearly 30 y ago, Tasmanian devil populations have declined precipitously. Despite the importance of host–pathogen coevolutionary interactions in epidemiological disease progression, coevolution has not yet been investigated in this system. We show that interactions between devil and DFTD genomes are more important in explaining variance in how quickly susceptible devils become infected than either genome alone, providing evidence of (or minimally the potential for) devil–DFTD coevolution. Our results demonstrate the feasibility, and importance, of considering intergenomic interactions when investigating disease dynamics in natural systems. Overall, we provide insight into devil–DFTD (co)evolutionary dynamics and can thus inform both epidemiological models and devil management strategies.

Coevolution is common and frequently governs host–pathogen interaction outcomes. Phenotypes underlying these interactions often manifest as the combined products of the genomes of interacting species, yet traditional quantitative trait mapping approaches ignore these intergenomic interactions. Devil facial tumor disease (DFTD), an infectious cancer afflicting Tasmanian devils (Sarcophilus harrisii), has decimated devil populations due to universal host susceptibility and a fatality rate approaching 100%. Here, we used a recently developed joint genome-wide association study (i.e., co-GWAS) approach, 15 y of mark-recapture data, and 960 genomes to identify intergenomic signatures of coevolution between devils and DFTD. Using a traditional GWA approach, we found that both devil and DFTD genomes explained a substantial proportion of variance in how quickly susceptible devils became infected, although genomic architectures differed across devils and DFTD; the devil genome had fewer loci of large effect whereas the DFTD genome had a more polygenic architecture. Using a co-GWA approach, devil–DFTD intergenomic interactions explained ~3× more variation in how quickly susceptible devils became infected than either genome alone, and the top genotype-by-genotype interactions were significantly enriched for cancer genes and signatures of selection. A devil regulatory mutation was associated with differential expression of a candidate cancer gene and showed putative allele matching effects with two DFTD coding sequence variants. Our results highlight the need to account for intergenomic interactions when investigating host–pathogen (co)evolution and emphasize the importance of such interactions when considering devil management strategies.

coevolution
co-GWAS
genomics
host–pathogen
joint phenotype
National Science Foundation (NSF) 100000001 DEB 2027446 Hamish Ian McCallumRodrigo HamedeMenna E. JonesAndrew StorferMark J Margres Harvard University (Harvard College) 100007229 Sarah and Daniel Hrdy Visiting Fellowship in Conservation Biology Mark J Margres Australian Research Council (ARC) Future Fellowship FT100100250 Rodrigo HamedeMenna E. Jones Australian Research Council (ARC) Large Grants A00000162 Rodrigo HamedeMenna E. Jones Australian Research Council (ARC) DE 170101116 Rodrigo HamedeMenna E. Jones Australian Research Council (ARC) LP 170101105 Rodrigo HamedeMenna E. Jones Agence Nationale de la Recherche (ANR) 501100001665 ANR Blanc Project TRANSCAN Rodrigo Hamede Centre National de la Recherche Scientifique (CNRS) 501100004794 International Associated Laboratory Grant Rodrigo Hamede University of South Florida (USF) 100008900 USF Ryan McMindsMark J Margres National Science Foundation (NSF) 100000001 DEB 2324456 Hamish Ian McCallumRodrigo HamedeMenna E. JonesAndrew StorferMark J Margres
==== Body
pmcTransmissible or infectious diseases greatly impact biodiversity (1), and characterizing the genetics underlying disease-related traits (i.e., genomic architecture) in both host and pathogen is essential for developing effective treatments (2), anticipating disease spread (3, 4), and predicting long-term epidemiological trajectories by understanding how the host and pathogen may evolve (5, 6). Genome-wide association studies (GWAS) are critical for determining the genetics underlying disease-related traits, associating specific variants with phenotypes such as host susceptibility (7), pathogen infectivity (8), and pathogen virulence (9). Typical GWA models, however, almost exclusively consider a single genome (10–12), ignoring joint or extended phenotypes (e.g., refs. 13–15) that are the products of interactions between the genomes of organisms which may be coevolving. Coevolution is particularly important for host–pathogen relationships, as rapid reciprocal evolution frequently governs the epidemiological progression of a disease (16–18). Although new co-GWAS techniques incorporating both host and pathogen genomes have been developed for characterizing the genetic basis of disease-related joint phenotypes (18–22), such approaches have yet to be applied to natural systems because of sampling limitations and the lack of genetic resources in most nonmodel species (23).

The Tasmanian devil (Sarcophilus harrisii; Fig. 1A) has been afflicted by devil facial tumor disease (DFTD), a species-specific transmissible cancer which likely emerged circa 1985 (24, 25) and was discovered in 1996 (26). DFTD is nearly 100% fatal, is found throughout most devil populations in Tasmania, and has caused total population declines of ~80% (27). DFTD, which typically generates large tumors on the face and mouth of devils, is transmitted as a clonal allograft when an uninfected devil bites the tumor mass of an infected devil (28). Devil biting behavior is common during social interactions, making DFTD transmission predominantly frequency-dependent because such interactions typically do not depend on the density of devils (29, 30). The combination of universal devil susceptibility to DFTD, high devil mortality, and frequency-dependent transmission led early models to predict devil extinction (29), yet long-diseased populations persist. Indeed, newer models now predict devil–DFTD coexistence at a higher probability than devil extinction (31), suggesting coevolutionary interactions will be key drivers of disease dynamics and thus important for long-term disease management.

Fig. 1. (A) Photo of a Tasmanian devil (photo by R.H.). (B) Pie charts show sampling locations, including the proportion of individual devils (blue), DFTD (red), and devil–DFTD pairs (purple) collected at each site. Pie chart size indicates relative sample size at each site. Primary sampling sites are labeled as FR: Freycinet, WPP: West Pencil Pine, TK: Takone, BR: Black River, and AR: Arthur River. (C) Violin plots representing the posterior distributions of the total variance in force of infection explained by all SNPs (PVE) and the total variance in force of infection explained only by large-effect SNPs (PGE) for devils (n = 313) and DFTD (n = 331) obtained through Bayesian Sparse Linear Mixed Model fitting. Dots are median values and lines are Bayesian 95% credible intervals.

Devils and DFTD represent an ideal natural system for the application of a co-GWAS. First, extensive mark-recapture data exist for devils over a >15-y period, providing a wealth of ecological information and phenotypic data (32–36). Second, strong linkage disequilibrium (LD) in the devil genome [~200 kb (34) and small devil population sizes facilitate variant detection with less sampling of both the population and individual genomes (23). Third, DFTD is species-specific (37), eliminating the possibility of a diffuse coevolutionary interaction caused by multiple hosts. Fourth, the high virulence of DFTD (27) likely imposes a strong selective pressure on devils, potentially leading to a strong reciprocal selective pressure on DFTD as devils respond to this high virulence. Fifth, the high mutation rate in DFTD (25) produces sufficient de novo variation to be detectable through GWAS, despite asexual reproduction in the cancer (Materials and Methods and SI Appendix). Last, evidence of evolution has been found in both devils (34, 36, 38, 39) and DFTD (24, 25, 40) independently, strongly supporting the possibility of coevolutionary interactions.

Here, we investigated genomic signatures of devil–DFTD coevolution underlying a force of infection proxy, a measure of how quickly, in days, susceptible devils became infected with DFTD. Our force of infection proxy is analogous to the inverse of force of infection (Materials and Methods) and is hereafter referred to as FOI. For devils, FOI measures resistance to infection, whereas for DFTD, FOI measures infectivity. Hence, FOI represents an extended or joint phenotype (41) as selective pressures on host resistance can incur reciprocal selective pressures on pathogen infectivity (or vice versa). Because FOI is an extended phenotype, its genetic basis is expected to be at least partially determined by intergenomic interactions, whereby intergenomic interactions will likely explain a greater proportion of FOI with stronger reciprocal selective pressures. We used a genome-wide capture panel to sequence 456 devils and 504 tumors (339 of which were matched devil-tumor pairs) collected throughout Tasmania (Fig. 1B). First, we used a traditional GWA approach to determine the variance in FOI that could be explained independently by devil and DFTD genomes. Next, we used Analysis with a Two-Organism Mixed Model (ATOMM), a two-way mixed-effects model developed to simultaneously analyze host and asexual pathogen genomes (20), to look at the contribution of devil–DFTD intergenomic interactions to explaining FOI variance. Finally, we characterized the top genotype-by-genotype interactions to test whether the underlying genetics were biologically meaningful (i.e., cancer-related and/or exhibited signatures of selection) and thus supported the relevance of the identified intergenomic interactions underlying variation in FOI. By including devil and DFTD samples from multiple trapping sites across Tasmania, similar to previous work in devils (34), we were able to identify range-wide signatures of coevolution through concordant allele frequency changes across multiple host populations (42–44) and/or distinct DFTD lineages (24, 40).

Results and Discussion

Traditional GWA Mapping Identifies Differences in Host and Pathogen Genomic Architectures.

To determine the genomic architecture underlying FOI for devils and DFTD independently, we first used a traditional, single-genome GWA Bayesian model (45) and found that a substantial proportion of the variance in FOI could be explained by devil (0.501; 95% CI = 0.248–0.750) and DFTD (0.342; 95% CI = 0.083–0.818) genomic variation independently (Fig. 1C). However, the genomic architectures underlying FOI differed across the genomes. The devil component of FOI was characterized by few loci of large-effect (~19% of phenotypic variance explained by ~2 SNPs; 95% CI = 1–4 SNPs), consistent with previous work on devil survival following infection (36), whereas the DFTD component of FOI was much more polygenic (~16% of phenotypic variance explained by ~35 SNPs; 95% CI = 1–238 SNPs).

The simpler genomic architecture identified for devil FOI may be explained by selection acting on standing genetic variation and the strong selective pressures DFTD imposed on devils (27). The recency of DFTD emergence [≤10 devil generations (24, 25) makes it unlikely that novel beneficial mutations have arisen in devils. Indeed, a previous genomic study found evidence of soft selective sweeps (i.e., positive selection rapidly fixing beneficial alleles from standing variation) in devils at cancer-related genes (34). Given our sampling at or near disease arrival (≤10 devil generations; SI Appendix, Fig. S1) at each respective trapping site along with universal susceptibility to DFTD (37) and high devil mortality (26), the large effect size of devil alleles for FOI suggested that devils may have initially been far from their phenotypic optimum upon DFTD emergence; such a scenario would be expected for a novel pathogen (46) and also would have facilitated large-effect alleles initially outcompeting smaller-effect alleles (47). In contrast, because DFTD was founded by a single clone, subsequent evolution and adaptation only occurred through de novo mutations in DFTD. The lack of large-effect mutations in DFTD likely occurred for several reasons. First, DFTD may have been closer to its phenotypic optimum for FOI upon emergence relative to devils given host naivety, and new beneficial mutations thus conferred small fitness gains. Second, novel DFTD mutations which are mildly deleterious cannot be purged due to the cancer’s clonal nature (but see ref. 48). Here, purifying selection would have removed any large-effect deleterious mutations such that only deleterious mutations of small effect size were maintained. Therefore, the polygenic architecture of DFTD may be the result of mildly beneficial and mildly deleterious mutations accumulating in different cancer lineages.

Next, we used per SNP posterior inclusion probabilities (PIP; the probability of being a large-effect SNP; SI Appendix, Fig. S2A) to identify specific candidate loci associated with variation in FOI in devils and DFTD. Because of extensive LD in devils [~200 kb (34)] and DFTD (clonal), specific loci associated with FOI represent markers for larger genomic regions potentially linked with a causative variant. Here, we focused on SNPs located within cancer-implicated genes which were thus at a higher likelihood of being putatively associated with the phenotype. Because DFTD is a cancer, cancer-implicated genes likely play a large role in FOI variability, and we thus expected many candidate variants to be found within cancer genes. In devils, 94 of the 100 SNPs with the highest PIP were intergenic or located within unannotated genes, with 5 of the 6 genic SNPs found within the introns of cancer-implicated genes (rank = 24–97; PIP = 0.0104–0.0041). For example, the devil SNP with the highest PIP in a cancer-implicated gene (PIP = 0.0104; rank = 24) was located in CNTN4, the overexpression of which has been found to be associated with malignant neuroendocrine tumors (49). In DFTD, 82 of the 100 SNPs with the highest PIP were similarly located in intergenic regions or unannotated genes. Here, 13 of the 18 DFTD genic SNPs were found in cancer-implicated genes (rank = 8–99; PIP = 0.0620–0.0087), and all such SNPs were within introns (with one of these intronic SNPs in a splice region) other than a single 5′ UTR SNP. The cancer-implicated SNP with the highest PIP in DFTD (PIP = 0.0620; rank = 8) was located in the 5′ UTR of GLI2, an oncogene which has been found to be overexpressed in human glioblastomas (50). Given these specific cancer-implicated SNPs in both devils and DFTD had overall effect sizes (i.e., the sum of small and large effects) in the >99th percentile, these loci represent biologically relevant candidate genomic regions associated with FOI variance.

Intergenomic Interactions Explain More Trait Variance than Either Genome Alone.

Although the traditional GWA framework used above revealed that both devil and DFTD genomic variation explained significant variation in FOI, such approaches do not account for interactions between genomes. We next used ATOMM (20), a co-GWAS approach, to estimate the phenotypic variance explained (PVE) for FOI by devil–DFTD intergenomic interactions. We found that devil–DFTD intergenomic interactions explained more variation in FOI (median PVE = 0.317; 95% CI = 0.293–0.330) than either devil (median PVE = 0.122; 95% CI = 0.108–0.128) or DFTD (median PVE = 0.084; 95% CI = 0.0829–0.0848) genomes alone (Fig. 2A). The clear importance of devil–DFTD intergenomic interactions in explaining FOI implicates FOI as a probable joint phenotype between devils and DFTD, providing evidence that devil–DFTD coevolution has either occurred or minimally has the potential to occur.

Fig. 2. (A) Proportion of phenotypic variance for force of infection explained by devil genomes, DFTD genomes, and devil–DFTD intergenomic interactions obtained through ATOMM joint model fitting (n = 100 ATOMM runs). (B) Top 250 most significant devil–DFTD genotype-by-genotype interactions (based on P-value). Each line represents a devil–DFTD interaction, with blue ends showing the location of the interacting site in the devil genome, and red ends showing the location of the interacting site in the DFTD genome. Asterisks show the three regions with the highest interaction density for devils (blue) and DFTD (red). (C) Enriched gene ontology (GO) terms for SNPs falling within interaction hotspots. Chromosomes are shown on the bottom of the plot and GO terms on the top. Connecting lines indicate chromosome but not the position of the enriched SNPs for devils (Top) and DFTD (Bottom), independently. (D) Relationship of genotype-by-genotype interactions for two of the top cancer x cancer interactions. For each interaction, “0” indicates the reference allele and “1” indicates at least one alternate allele (i.e., heterozygous or homozygous alternative). Because all alternate alleles in DFTD were found in samples exclusive to Freycinet, only samples from Freycinet are shown (n = 71; all data are shown in SI Appendix Fig. S14). Asterisks indicate P < 0.05 based on Tukey HSD testing. (E) Differential expression of ANK3 in devils based on RNA-seq data. The “0/0” category indicates homozygous for the reference allele, “0/1” indicates heterozygous, and “1/1” indicates homozygous for the alternate allele. “ns” indicates nonsignificance. All boxplots show the following: median (black horizontal line), 25th percentile (bottom edge of the box), 75th percentile (top edge of the box), 1.5 * interquartile range (IQR: 75th percentile - 25th percentile, shown as the upper and lower whiskers), and outliers (shown as dots above and below the whiskers). TPM: transcripts per million. Devil GO term abbreviations: G1, regulation of lipid metabolism; G2, regulation of lipid synthesis. DFTD GO term abbreviations: G1, biosynthesis; G2, regulation of biosynthesis; G3, cellular biosynthesis; G4, organic substance biosynthesis; G5, heterocycle biosynthesis; G6, cellular nitrogen compound biosynthesis; G7, nucleobase-containing compound biosynthesis; G8, regulation of metabolism; G9, regulation of primary metabolism; G10, regulation of macromolecule metabolism; G11, organic cyclic compound metabolism; G12, cellular aromatic compound metabolism; G13, cellular nitrogen compound metabolism; G14, regulation of nitrogen compound metabolism; G15, nucleobase-containing compound metabolism; G16, nucleic acid metabolism; G17, RNA metabolism; G18, regulation of RNA metabolism; G19, regulation of RNA synthesis; G20, gene expression; G21, regulation of gene expression; G22, regulation of nucleic acid-templated transcription; G23, DNA-templated transcription; G24, regulation of DNA transcription; G25, G-protein-linked receptor activity.

Total PVE for FOI explained by ATOMM (0.523) was more than the PVE explained for devils (0.501) and DFTD (0.342) by the traditional single-genome models independently. However, the traditional, single-genome models described above explained more variance relative to the ATOMM devil PVE estimate (0.122 relative to traditional single-genome devil PVE = 0.501) and DFTD PVE estimate (0.084 relative to traditional single-genome DFTD PVE = 0.342), suggesting that the traditional single-genome model PVEs were inflated, and a portion of the PVE for this traditional GWA approach was likely attributable to intergenomic interactions that the traditional single-genome method could not model. Indeed, ATOMM partitioned variance between devil genomes, DFTD genomes, and devil–DFTD intergenomic interactions within a single model compared to the traditional genome model estimates which partitioned variance into just devil or just DFTD genomes in separate models. As such, the total PVE for the traditional genome models ranged from 0 to 2 relative to a range of 0 to 1 for the ATOMM model, indicating that the traditional single genome models explained less PVE = 0.423 (0.843/2) than the total PVE explained by ATOMM (0.523).

We further investigated the differences between the traditional GWA and co-GWA approaches by comparing the top 100 single-genome SNPs (based on the highest PIP) for devils (top 0.05% SNPs) and DFTD (top 0.5% SNPs) to the top 0.01% SNPs (25,242) participating in genotype-by-genotype interactions (based on lowest P-value). We used a higher percentage of single-genome SNPs to ensure an adequate number of SNPs were compared to co-GWAS SNPs. We found limited overlap between the single-genome and co-GWAS SNPs (30/100 overlapping SNPs in devils and 13/100 overlapping SNPs in DFTD; SI Appendix, Fig. S3). The low number of overlapping SNPs was likely attributable to genotype-by-genotype interactions that could only be explicitly modeled and detected in the co-GWA approach (SI Appendix, Fig. S4), but this fails to explain why any SNPs were found to overlap between the single-genome model and the co-GWA model. Given that devil–DFTD intergenomic interactions likely partially contributed to the single-genome model PVEs, it is expected that at least a small fraction of the top SNPs in the single-genome model would be represented in the co-GWA model. Indeed, some of the single-genome SNPs perhaps had inflated effect sizes partially attributable to intergenomic interactions with the other lineage which the single-genome approach could not explicitly model. Furthermore, few of the overlapping SNPs were significant (4/30 in devils and 3/13 in DFTD) under ATOMM single-genome testing (a traditional maximum likelihood GWA approach; SI Appendix, Fig. S2B), indicating that a portion of their effect sizes were likely attributable to intergenomic interactions. Thus, considering genotype-by-genotype interactions relative to host- or pathogen-only models is necessary to unravel the genomic underpinnings of joint phenotypes, emphasizing the importance of co-GWA approaches.

The Most Significant DFTD SNPs Are Uniquely Informative Despite Clonal Reproduction.

Asexually reproducing organisms present unique challenges to GWAS due to their lack of genetic recombination through sexual reproduction, often giving rise to extensive LD and, particularly in pathogens, strong population stratification (51–53). Given that DFTD is clonal, extensive and perhaps even genome-wide LD may have reduced the resolution with which potential regions possessing a candidate variant could be identified (but see Materials and Methods). In other words, if all variants were in LD within DFTD then nearly all variants should have been equivalently informative in estimating PVE, irrespective of their actual contributions to phenotypic variance. Hence, we investigated the contribution of our most significant SNPs to our PVE estimates within the ATOMM framework. First, to determine whether our most significant single-genome DFTD SNPs identified by ATOMM were informative in estimating PVE (Fig. 2A), we removed the most significant DFTD-only SNPs and reran ATOMM (SI Appendix, Fig. S5). After removing the top 5% most significant DFTD SNPs based on ATOMM single-genome testing (i.e., SNPs expected to be most important in estimating DFTD-only PVE), we found that the DFTD PVE estimate dropped from 0.084 to approximately zero (0.004; SI Appendix, Fig. S5A). Next, we removed the top 5% most significant DFTD SNPs based on ATOMM interaction testing (i.e., SNPs expected to be most important to the interaction PVE) and, consistent with our results above, identified an interaction PVE of zero (relative to an initial estimate of 0.317; SI Appendix, Fig. S5B). Similarly, removing both the top 5% most significant single-genome DFTD SNPs and top 5% most significant interaction DFTD SNPs yielded a DFTD PVE near zero (0.020) and an interaction PVE of zero (SI Appendix, Fig. S5C). The clear pattern of reduced PVE estimates after removal of significant DFTD SNPs indicated that these SNPs were not in complete LD with nonsignificant SNPs and, therefore, all variants were not equally informative. Finally, to ensure that the lower PVE estimates were not a consequence of simply removing SNPs from the model, we removed 10% of DFTD SNPs at random and found a negligible change in PVE estimates from our initial analysis shown in Fig. 2A (host PVE Δ = 0.005, pathogen PVE Δ = 0.006, interaction PVE Δ = 0.033; SI Appendix, Fig. S5D). Negligible PVE changes would be expected only if a small fraction of DFTD SNPs were informative, as removing SNPs at random would more likely remove uninformative SNPs than the few informative SNPs. Thus, our most significant DFTD SNPs were highly informative relative to our nonsignificant SNPs, representing candidate loci associated with variation in FOI. Indeed, although identifying causal variants was not one of our goals, these significant candidate SNPs were minimally yet uniquely in LD with putative causal variants and thus reliable genomic markers of variation in FOI.

Although the most significant DFTD SNPs were informative as described above, a previous study identified at least four distinct DFTD lineages using low coverage sequencing (40). To ensure that DFTD samples within specific lineages (i.e., genetically similar DFTD samples derived from the same clonal “parent”) were not in genome-wide LD and affecting the ATOMM results, we needed to assess the degree of genetic variability within tumor lineages. We first confirmed the existence of four distinct DFTD lineages within our dataset (SI Appendix, Fig. S6C) with a relatively homogenous spatial distribution (SI Appendix, Fig. S6D), similar to previous work (40). To quantify intralineage variation, we next observed how often all ATOMM samples within a lineage possessed the same allele. Across the four lineages we identified, one or more samples possessed a different allele relative to the other samples within the lineage at 79 to 88% (9,791 to 10,904) of the 12,420 ATOMM SNPs. For each lineage, a median of ~4.5% of the samples (n = 66–107) differed in their allele at each SNP (SI Appendix, Fig. S7), corresponding to a median of 3–5 samples possessing a different allele at each site. Finally, we quantified the number of SNPs that varied between the most similar tumors within a particular lineage and found that 12 to 17% (1,456 to 2,114; SI Appendix, Fig. S8) of alleles differed, indicating that even the most similar tumors differed by >10% of the SNPs included as input to ATOMM. Clearly, despite asexual reproduction in DFTD, sufficient genetic variation exists both within and between lineages to be informative to GWAS.

Intergenomic Interaction Hotspots Are Biased toward Cancer Genes and Enriched for Selection Signatures in Both Host and Pathogen.

The most significant devil–DFTD genotype-by-genotype interactions were nonrandomly distributed across the genomes of the two lineages, as interacting SNPs were often clustered together in the genome (i.e., “interaction hotspots”; Fig. 2B). Scanning of the devil genome revealed hotspots on all chromosomes, although there were no devil interaction hotspots on any unplaced scaffolds, and chromosomes 5 and 6 lacked DFTD interaction hotspots (SI Appendix, Fig. S9 A and B). SNPs within DFTD hotspots, however, may actually have been located on different chromosomes due to chromosomal rearrangements in the tumor (54) and using the devil genome as our reference (Materials and Methods), although the probability of this occurring was likely decreased due to small hotspot window size (1 Mb). To ensure hotspots were not simply biased toward regions with more genetic variation, we regressed interaction count on SNP density within hotspots and found that SNP density did not predict interaction count for devils (P = 0.335; SI Appendix, Fig. S9C) or DFTD (P = 0.946; SI Appendix, Fig. S9D). Additionally, LD within hotspots did not differ significantly from nonhotspot LD in devils (P = 0.259; SI Appendix, Fig. S10 A and B) or DFTD (P = 0.091; SI Appendix, Fig. S10 C and D). Gene ontology (GO) enrichment analysis found that interaction hotspot SNPs were significantly enriched for lipid metabolism and synthesis in devils (Fig. 2C). Altered lipid metabolism in host cells constituting the tumor microenvironment can facilitate tumor growth, invasion, and immunosuppression (55), suggesting that lipid metabolism in devils may be implicated in susceptibility and/or transmission. In DFTD, significant GO terms were associated with biosynthetic processes (G1–G19), gene expression (G20–G24), and G-protein-linked receptor activity (G25; Fig. 2C). Notably, GO terms associated with nucleotide metabolism were enriched in DFTD (e.g., G7; nucleobase-containing compound biosynthesis). In cancers, nucleotide metabolism plays an important role in the interaction between the host immune system and the cancer (56), suggesting these interaction hotspots may be involved in devil–DFTD immune interactions.

Similar to the candidate variants we identified using the single-genome Bayesian approach (see above), the candidate genotype-by-genotype interaction variants may have been in LD with putative causal variants and, therefore, may represent interactions between larger genomic regions associated with FOI variation rather than specific candidate genotypes. To increase the likelihood that at least some of our candidate interactions were associated with the phenotype, we first tested for cancer implications (described above) in our most significant genotype-by-genotype interactions. To test for an enrichment in our top interactions between candidate variants in putative cancer-implicated genes, we first obtained the most significant 0.001% (2,524 interactions; 5.53e-06 ≤ P ≤ 6.51e-04) of the 252 million devil–DFTD genotype-by-genotype interactions (SI Appendix, Fig. S11A). Next, we identified gene-by-gene interactions and found that 10.5% (264) of the most significant interactions occurred within genes for both devils and DFTD (i.e., within an intron, exon, or untranslated region; SI Appendix, Fig. S11B). Because we could only infer cancer implication in annotated genes, we retained interactions exclusively between annotated genes (96 of the 264 gene-by-gene interactions). We then used the DisGeNET database (57) to find genes putatively implicated in cancer, identifying 88 of the 96 interactions as occurring between two cancer-implicated genes. To test for significant enrichment, we used Monte-Carlo simulations and found that the most significant gene-by-gene interactions occurred significantly more frequently between two cancer-implicated genes than what was expected by chance (P < 0.0001; SI Appendix, Fig. S12). The significant enrichment we found for cancer-by-cancer gene interactions likely indicated that at least some of our candidate genotype-by-genotype interactions contributed to variation in FOI.

Next, we looked for signatures of selection in our candidate genotype-by-genotype interactions using pcadapt (58) and found a significant enrichment for selection in both devil (P < 0.0001; 22.1% candidates under selection versus 13.9% noncandidates; SI Appendix, Fig. S13B) and DFTD (P < 0.0001; 66.8% candidates under selection versus 24% noncandidates; SI Appendix, Fig. S13D) candidates before and after LD thinning relative to noncandidates (see SI Appendix, Materials and Methods for details). The total proportion of variants under selection was significantly lower in devils (14.1%; SI Appendix, Fig. S13A) relative to DFTD (26.1%; SI Appendix, Fig. S13C; P < 0.0001), consistent with their differing genomic architectures (Fig. 1C). Overall, the candidate genotype-by-genotype interactions likely represent biologically relevant markers for variation in FOI because these candidate intergenomic interactions 1) explained substantial variation in FOI (Fig. 2A and SI Appendix, Fig. S5B), 2) were significantly biased toward cancer-implicated genes (SI Appendix, Fig. S12), and 3) were significantly enriched for signatures of selection relative to noncandidates (SI Appendix, Fig. S13 B and D).

Allele-specific Effects in Candidate Genotype-by-genotype Interactions Are Associated with Phenotype Differences.

Only two devil–DFTD candidate cancer interactions occurred between genic SNPs outside of an intron: A devil SNP within the 5’ UTR of ANK3 interacted with a nonsynonymous DFTD SNP in ARHGAP5 (P = 3.84e-04), and the same ANK3 devil SNP interacted with a DFTD nonsynonymous SNP in ACAP1 causing a stop-gain mutation (P = 5.39e-04). ANK3 is important to cytoskeleton membrane stability (59), and ANK3 knockdown was found to promote cancer cell invasion in prostate cancer (60). Reduced expression of ARHGAP5, which may mediate cytoskeletal changes through RHO GTPases, has been shown to inhibit cancer cell migration and invasion (61). ACAP1, implicated with GTPase functioning, is up-regulated in ovarian cancers and associated with poor prognosis (62). Interestingly, these SNPs were nonsignificant when examining the effects of the devil and DFTD genomes separately using single-genome tests in ATOMM (ANK3 P = 0.602, ARHGAP5 P = 0.883, and ACAP1 P = 0.600) or GEMMA (ANK3 PIP < 1e-6, ARHGAP5 PIP = 0.0018, and ACAP1 PIP = 0.0020; SI Appendix, Fig. S2 A and B), again highlighting the importance of co-GWAS approaches.

We next determined the trapping site and tumor lineage for samples possessing an alternate allele for the ANK3, ARHGAP5, and ACAP1 SNPs. We found that DFTD samples with an alternate SNP belonged to lineages 1 and 4 (SI Appendix, Table S1), and all tumors were located exclusively in Freycinet, a single trapping location in eastern Tasmania (FR in Fig. 1B); we thus restricted our subsequent analyses to include only Freycinet samples to further account for population structure (Materials and Methods). For both the ANK3 x ARHGAP5 and ANK3 x ACAP1 interactions, the presence of alternate alleles in devils and DFTD was associated with a significantly longer time period before a devil became infected in Freycinet than when the alternate allele occurred in only DFTD (P = 0.006) or neither lineage (P = 0.045; Fig. 2D). The lack of significant difference between FOI associated with alternate alleles in both lineages and an alternate allele only in devils (P = 0.279) was likely due to low sample sizes (n = 3 in the “1 x 0” category; Fig. 2D). Running the same analysis with all devil–DFTD paired samples across all sites yielded qualitatively similar results (SI Appendix, Fig. S14), suggesting that ATOMM effectively accounted for population structure of the host and pathogen (Materials and Methods). Interestingly, the candidate alleles identified in both genes were found in two distinct DFTD genetic lineages (SI Appendix, Table S1); whether the alleles arose once in a common ancestor or multiple times independently was unclear and warrants exploration in a future study.

Overall, we detected significant intergenomic interactions underlying how quickly susceptible devils became infected at three different scales: 1) across the entire genome in the form of interaction hotspots that were enriched for GO terms associated with cancer (Fig. 2 B and C), 2) at the genic level with a significant enrichment toward cancer x cancer interactions and selection signatures, and 3) putatively between specific candidate alleles in cancer-related genes that showed a clear association with how quickly susceptible devils became infected in a single locality (Fig. 2D).

A Candidate Allele Is Associated with Differential Expression of Cancer Genes in the Tasmanian Devil.

To determine whether the SNP in the 5′ UTR of ANK3 in devils putatively affected gene expression, we used publicly available RNA-seq data (n = 20) (63) and found that the presence of the 5′ UTR alternate allele was associated with downregulation of ANK3 (Fig. 2E). Although ANK3 knockdown in human cancer cells promoted cancer progression (61), we found that decreased expression in devils was associated with an increase in the time it takes for a devil to become infected (when matched with an alternate allele in DFTD; Fig. 2D). Given that the SNP observed in ANK3 was in devils rather than the cancer, different effects on cancer progression relative to what was observed in a human cancer were unsurprising. Ultimately, the change in gene expression we found associated with the ANK3 SNP suggested that regulatory alleles in devils may have interacted with coding-sequence alleles in DFTD, indicating that genotype-by-genotype interactions can perhaps occur between different classes of mutations.

Conclusion

Collectively, our results revealed that devil–DFTD intergenomic interactions were vital to explaining variability in how quickly susceptible devils became infected, a key disease-related joint phenotype. Although GWA models are necessarily associative, we showed that the top associated genotype-by-genotype interactions were significantly enriched for cancer genes and selection, providing evidence that at least some of the intergenomic interactions we identified were putative signatures of coevolution (or minimally provide the substrate for coevolutionary interactions). Our findings further indicate that the genetics underlying coevolution may occur between antagonists with apparently differing genomic architectures (i.e., polygenic versus simple) and between different classes of mutation. Yet, as evidenced by the few overlapping SNPs we found between the traditional single-genome approach and the co-GWA model, independent genomic architectures of coevolving lineages incompletely represent the genetic basis of coevolutionary interactions; rather, coevolutionary architecture appears as an emergent property only when simultaneously modeling the genomes of both lineages. Taken together, our study represents evidence that devils and DFTD are likely coevolving, providing an important step forward in understanding the genomic architectures underlying host–pathogen coevolution in general. Future work investigating differing spatiotemporal dynamics of devil–DFTD coevolution among sites and lineages may uncover important aspects of ongoing host–pathogen coevolution as well as their effects on epidemiological progression.

Unlike human cancers, DFTD can outlive its host and progresses without treatment, facilitating the study of long-term cancer evolution in vivo and transgenerational host-cancer coevolutionary dynamics. For example, the devil resistance trait we investigated here is analogous to cancer resistance in nontransmissible cancers, and, through comparative oncology (64–66), the genes we identified (such as ANK3 in devils) may represent useful candidates in future cancer studies on resistance or susceptibility to carcinogenesis and metastasis. The clear importance of devil–DFTD intergenomic interactions underlying a disease-relevant trait also have implications for devil conservation. Devil breeding programs may be able to use genetic information from both devils and DFTD to improve disease resilience in captive devil populations. If translocations are proven to be beneficial, specific individuals can be selected by accounting for intergenomic interactions between translocated devils and local strain(s) of DFTD. Ultimately, accounting for devil–DFTD intergenomic interactions is necessary in predicting their evolutionary trajectories (67–69) and thus making scientifically informed management decisions.

Materials and Methods

Sample Collection.

Tasmanian devils were trapped from 2006–2020 at five primary sites across Tasmania (Fig. 1B; see SI Appendix, Fig. S1 for trapping date distribution and SI Appendix, Table S3 for DFTD arrival date at each site). Devils were trapped under a capture-mark-recapture framework, using PVC culvert pipe traps to catch devils and microchip transponders to permanently mark individuals. Standard measures, including sex, estimated date of birth, disease status, and tumor dimensions (if present), were collected (70). Devils were aged using a combination of teeth eruption and wear (71). Tissue samples for sequencing were collected as ear or lip biopsies for all devils upon first capture, and DFTD tissue samples were obtained from tumor margins each time a diseased devil was trapped. Tumor samples were confirmed as DFTD by either macroscopical examination, histopathological examination, and/or PCR (72). Sample collection and animal work was performed under University of Tasmania ethics approval A13326, WSU IACUC approval ASAF 6796, and USF IACUC approval IS00008158.

Capture Panel-design and Genome Library Preparation.

To effectively sample across the genome, we designed a target-capture panel of ~197k probes (Twist Bioscience, San Francisco, CA). Such an approach should map a significant proportion of both genomes owing to extensive LD in devils (~200 kb) and DFTD (clonal), similar to previous work which successfully applied GWAS to devils using a sparser capture panel than what was used here (36). Specifically, we ensured our target-capture panel was sufficiently dense to capture at least one SNP per linkage block in devils and DFTD, providing enough genetic variation to identify genotype–phenotype associations using GWAS (23). Probes were chosen to target a) loci under selection based on previous work in both devils (34, 36, 39, 73) and DFTD (74–76), b) cancer-related genes from the PanCancer atlas (77) and c) sufficient coverage of the genome with backbone tiling targeting a SNP every ~17 kb.

We extracted DNA from devil and tumor tissues using the Qiagen DNeasy Blood and Tissue 96 Kit, including an RNase A (ThermoFisher Scientific) treatment and elution in Qiagen Buffer EB. Prior to library preparation, we sheared DNA samples using a Covaris M220 focused ultrasonicator and microTUBE-50 AFA fiber screw cap tubes to achieve a modal fragment size of ~200 to 250 bp. Using 50 ng sonicated DNA per sample, we performed library preparation using the Twist Mechanical Fragmentation and Twist Universal Adapter System (Twist Bioscience). Following validation using an Agilent Bioanalyzer, libraries were multiplexed into groups of eight and enriched for our target capture panel using the Twist Target Enrichment Protocol. We used salmon sperm DNA (Invitrogen) as a blocking solution as our trials indicated marginal improvement in capture performance over Twist’s provided solution. Enriched libraries were again validated using a Bioanalyzer before further multiplexing into sequencing pools containing 24 capture reactions each (192 total libraries per sequencing pool).

Genomic Sequencing.

Libraries were sequenced 150-bp PE on an Illumina NovaSeq 6000 in batches of 192 individuals per S4 lane (5 total batches; 960 total libraries) to ~35× coverage per library (range of 28× to 64×). Sequencing was performed at the Genomics Sciences Laboratory at North Carolina State University. All raw data were deposited in the Sequence Read Archive (SRA) under BioProject PRJNA947341 and accessions SRR23935061–SRR23936020 (78).

Genomic Alignments and Variant Calling.

Raw read quality was assessed using FastQC (79), and reads were trimmed with TrimGalore! (80) at default settings. Trimmed reads were aligned to the reference genome mSarHar1.11 (25) using BWA MEM version 0.7.17 (81) with the -M flag and default settings, and Samtools (82) was used to sort the aligned reads. Picard version 2.25.0 (83) was used to mark PCR and optical duplicate reads via MarkDuplicates. GATK version 4.2.0.0 (84) was used to call SNPs. Here, HaplotypeCaller was run with the “-ERC GVCF” and “–do-not-run-physical-phasing” flags for each individual sample. To combine the output from HaplotypeCaller, GenomicsDBImport was run separately for devils and DFTD with intervals set to each chromosome of the devil reference assembly (including unplaced scaffolds). GenotypeGVCFs was then used to extract SNPs and indels from the DFTD and devil databases, generating lineage-specific variant call format (VCF) files across all samples. GenotypeGVCFs was run in parallel such that extraction was done per-chromosome and combined in the order of each chromosome once all runs had completed. Next, SelectVariants was used to generate SNP files separately for each lineage (indels were excluded), and VariantFiltration was run with filters QD < 2.0, FS > 60.0, MQ < 40.0, MQRanksSum < -12.5, and ReadPosRankSum < 8.0 as recommended by GATK developers. To remove potential host contamination in the DFTD SNP set (i.e., obtain a set of tumor-specific somatic mutations), tumor SNPs were further filtered using BCFtools isec (82) to remove any SNPs common to both the devil and tumor VCF files. The resulting VCF files contained 7,636,616 devil SNPs and 6,183,694 DFTD SNPs and were further filtered based on parameters specific to each analysis as described below (see SI Appendix, Table S2 for analysis-specific sample and SNP counts).

Estimating the Force of Infection Proxy.

Traditionally, force of infection measures the rate of infection over an entire population, representing an average value calculated over many individuals (85). Our force of infection proxy differed in two important ways from a typical force of infection calculation: 1) We calculated force of infection for specific individuals rather than as a population average, and 2) our force of infection proxy was calculated as the time until disease after exposure, making it analogous to the inverse of true force of infection (i.e., 1/FOI). Specifically, our force of infection proxy was calculated as the length of time (in days) it took for a devil to become infected with DFTD when infection was possible based on devil age (SI Appendix, Materials and Methods) and date of disease arrival (26, 27, 35) (SI Appendix, Table S3).

Because tumor size varies between individuals, it was necessary to estimate the date of first infection for individual sequenced tumors. A model of tumor growth, based on a logistic growth curve, had previously been fit to data for a single population (West Pencil Pine) also included in this study (38); this model was used to backward project an initial tumor volume and estimate the date of first infection based on the tumor volume of the first trapping event where DFTD was observed. Because initial tumor volume is too small to estimate, initial tumor volume was assumed to be initial tumor volume + 1, thus shifting the logistic growth curve down one unit of tumor volume. This alteration allowed for an initial tumor volume of zero to be used, making the back calculation dependent on time rather than initial tumor volume. The fit model had an averaged estimated maximum tumor volume (Mmax) of 202 cm3 (95% CI = 198–223 cm3) and a lag phase of ~60 d during which it is assumed that a tumor is not observable upon trapping (38).

To estimate the age at which a devil was first infected with DFTD, only the volume of the sequenced tumor at the earliest trap date was used (see SI Appendix, Materials and Methods for details). To minimize potential error from estimating devil age at first infection using the tumor growth model, data were filtered in several ways. First, if the devil’s estimated date of initial infection was before the arrival date of DFTD, the sample was removed. In West Pencil Pine, a tetraploid DFTD strain with differing transmission dynamics was initially present in 2006 until being replaced by the typical diploid strain in 2011 (33). Thus, devils in West Pencil Pine with an estimated initial infection date before 2011 were removed. Second, to maximize power, tumor depth was imputed for samples missing a depth measurement (~2% of samples; see SI Appendix, Materials and Methods for details). Finally, because the Mmax for the growth model averaged 202 cm3, accurate estimates could not be made for tumor volumes significantly larger than this maximum, and backprojections greater than a year (the backprojection obtained at Mmax) were unreliable. As such, devils exceeding a tumor volume of 223 cm3 (the 97.5% tumor load Mmax CI (38) on their first trap with an estimated infection date >60 d from a previous trap (i.e., when a tumor would be visible) were removed. Phenotype files containing FOI estimates can be found at https://github.com/D-gallinson/Devil-DFTD-FOI-Coevolution (86).

Characterizing Force of Infection Genomic Architecture in Devils and DFTD.

Genome-wide Efficient Mixed Model Association [GEMMA (87)] was used to implement a Bayesian Sparse Linear Mixed Model [BSLMM (45)] to explore the genomic architecture underlying FOI separately in devils and DFTD. Model outputs include the proportion of phenotypic variance explained (PVE; the total variance in the phenotype explained by both small- and large-effect SNPs), proportion of genotypic variance explained (PGE; the proportion of PVE explained only by large-effect SNPs), a posterior inclusion probability (the probability of being a large-effect SNP; PIP) for individual SNPs, the number of large-effect SNPs, and three other hypermeters related to PVE, PGE, and the number of large-effect SNPs.

Because a BSLMM was fit separately to devils and DFTD, FOI was calculated slightly differently for each lineage. Briefly, when fitting the model to DFTD genomes, age of infection was calculated using individual tumor volumes as input to the tumor logistic growth model (38) (see the previous section). For devil genomes, the age at first infection was found using tumor load (the sum of all tumor volumes when multiple tumors were present) at the earliest trap date as input (see SI Appendix, Materials and Methods for details). For both devils and DFTD, the VCF file was first filtered to remove samples with a missing phenotype. SNPs were subsequently filtered on 5% missingness, the removal of sites which were nonbiallelic, a minor allele frequency (MAF) of 0.05 for devils, and a MAF of 0.01 for DFTD. To remove putative repetitive genomic elements, SNPs with a mean max depth above the 90th percentile of the filtered devil VCF (89.5×) and the 95th percentile of the filtered DFTD VCF (69.6×) were removed. Genotypes were converted into genotype likelihoods by multiplying each genotype by its likelihood value calculated by HaplotypeCaller (fitting a BSLMM to discrete genotypes did not meaningfully change the results; SI Appendix, Fig. S15). A linear link function was used (-bslmm 1) with a centered relatedness matrix, default priors, and five independent Markov Chain Monte Carlo (MCMC) chains with 15,000,000 iterations per chain (with a 1,500,000 burn-in). Within and between chain convergence was assessed using Rhat (88), bulk effective sample size (ESS), and tail ESS (SI Appendix, Table S4) implemented in Rstan (89) in R version 4.1.0 (90). Convergence was also checked visually by manually inspecting hyperparameter distributions and trace plots. See SI Appendix, Materials and Methods for details.

Host–pathogen Joint Genome Modeling in ATOMM.

We used ATOMM to characterize putative signatures of devil–DFTD coevolution by estimating 1) genomic heritability (i.e., PVE) of the marginal host (i.e., devil-only), marginal pathogen (i.e., DFTD-only), interaction of the host and pathogen genomes, and noise, and 2) genotype-by-genotype interactions between individual host–pathogen SNPs (SI Appendix, Materials and Methods). The FOI data used as input for ATOMM were derived in a tumor-centric manner, whereby infection age was estimated based on the volume of the sequenced tumor (rather than tumor load as described above for devils). However, because ATOMM accounted for both devil and DFTD genomes, tumors which were not the first to infect the devil were retained and first-infection status was used as a binary fixed effect in the model. After removing samples which were unpaired and/or lacked FOI data, we retained 312 devil–DFTD pairs (SI Appendix, Table S2). At each primary sampling site, trapping began within ~5 devil generations (devil generation = 2 years) of DFTD arrival and proceeded for a similar duration at all sites (median = 4 devil generations; SI Appendix, Fig. S1) other than Black River (trapping duration = 2 devil generations) due to the recency of disease arrival (SI Appendix, Table S3). Given the rapid evolutionary response found in devils after only a few generations (34), we chose to include Black River samples to maximize power despite the short trapping duration. ATOMM assumes the response phenotype to be multivariate Gaussian or binomial (20), and the FOI values were thus standardized using the RankNorm function in RNOmni (91) in R version 4.1.0 (90). Because ATOMM is incapable of handling missingness and requires haploid genotypes as input, the host and tumor VCF files were first filtered to remove sites with missing genotypes (SI Appendix, Materials and Methods) and/or sites that were not biallelic. Next, genotypes were converted into a matrix of zeros and ones, where zero indicated two reference alleles (i.e., homozygous reference) and one indicated the presence of at least one alternate allele (i.e., heterozygotes and homozygotes for alternate alleles were considered identical). Both haploid genotype matrices were then filtered using a custom script based on a MAF of 0.05 for devils and MAF 0.01 for tumors as described above. Finally, both files were then filtered to remove SNPs with a mean max depth above 67.3× for devils and 80.2× for DFTD (the 90th and 95th percentiles of the MAF-filtered devil and DFTD files, respectively) as described above.

The final filtered devil haploid genotype matrix had nearly ten times the number of SNPs relative to the DFTD genotype matrix (115,018 devil SNPs versus 12,420 DFTD SNPs). Asymmetry in SNP number may bias the PVE estimates toward the host because ATOMM uses Fisher’s infinitesimal approach to model polygenic random effects (20). Thus, to reduce the additive effects of devil SNPs, a single SNP within the devil haploid genotype matrix was selected at random within 100 kb windows using BCFtools (82) with the prune plugin. The process of randomly selecting SNPs was then repeated 100 times (with replacement for each new file), generating 100 devil input files (20,324 SNPs per file). ATOMM model fitting was done for each of the 100 devil input files with the single DFTD input file, generating 100 estimates of a) PVE for the devil genome, b) PVE for the DFTD genome, c) PVE for the interaction between the two genomes, and d) noise (i.e., variation unaccounted for in the model). Random initializations were used for each ATOMM run (SI Appendix, Materials and Methods) and the ATOMM convergence delta parameter was changed from the default value (0.01) to 0.0001 to prevent premature convergence before reaching the optimum. Because genotype-by-genotype interaction testing becomes computationally expensive for even a modest number of SNPs (e.g., 240 million tests for 20,000 host SNPs and 12,000 pathogen SNPs), genotype-by-genotype testing was suppressed when fitting ATOMM to the 100 devil inputs for PVE estimation.

Marginal Testing of Host and Pathogen Variants in ATOMM.

To look at the effects of individual SNPs from only the host genome or only the pathogen genome on FOI (i.e., marginal effect testing), ATOMM tests each genome using a traditional GWA approach: the host test constrains pathogen effects and host–pathogen interaction effects to zero, and the pathogen test constrains host effects and host–pathogen interaction effects to zero. Because ATOMM was run using 100 devil input files, there were 100 marginal effect p-values for each devil and DFTD SNP, respectively. These p-values were summarized using the harmonic mean p-value (SI Appendix, Fig. S2B), which corrects for false-positives generated from multiple testing while maintaining statistical power for dependent tests (92). The harmonic mean p-value was used because p-values across ATOMM runs for a given SNP are expected to be dependent.

Host–pathogen Genotype-by-genotype Interactions in ATOMM.

To identify genotype-by-genotype interactions significantly associated with variation in FOI, significance estimates were performed for all pairwise comparisons of devil and DFTD SNPs using ATOMM. For each individual interaction test, a P-value is generated which indicates the probability that the combined effects of two SNPs are modified by an additional interaction (20). Because the devil input was split into linkage groups, interaction tests were run ten times using ten random devil input genotype matrices (selected without replacement from the 100 devil input genotype matrices described above) crossed with the single DFTD genotype matrix. To summarize the results of these interaction tests across SNPs within a single linkage group in devils, the harmonic mean p-value (see above) was estimated across all ten interaction tests performed within that linkage group. Here, the harmonic mean p-value accounted for the dependent nature of SNPs in linkage disequilibrium.

Devil and DFTD Population Structure.

To ascertain devil and DFTD population structure, we separately clustered the 456 devil samples and 504 DFTD samples using discriminant analysis of principal components (DAPC; ref. 93) on a set of SNPs with equivalent filters to the BSLMM filters (see above). First, we used find.clusters (which clusters samples using k-means) to find the number of genetic clusters, selecting the number of clusters with the lowest Bayesian Information Criterion. Next, clusters were described using the dapc function, where we retained principal components (PCs) explaining 90% of the variance for both devils (350 PCs) and DFTD (410 PCs) and retained one discriminant axis for devils and three discriminant axes for DFTD. See SI Appendix, Materials and Methods for details.

Testing Whether the most Significant DFTD SNPs Were Informative to ATOMM PVE Estimates.

Despite the difficulties implicit to conducting GWAS on asexual organisms (i.e., potential genome-wide LD due to clonal reproduction; but see SI Appendix, Materials and Methods), many GWAS have been successfully carried out on asexually reproducing organisms (52, 53, 94–96). If DFTD is in genome-wide LD, however, then the most significant SNPs would be expected to be no more informative than the nonsignificant SNPs, and we would thus expect to see very little change in the PVE estimates after removing the significant SNPs. As such, to determine whether the most significant DFTD SNPs (as identified by ATOMM) were more informative than the nonsignificant SNPs (i.e., DFTD is not in genome-wide LD), ATOMM was rerun after removal of the top 5% most significant DFTD SNPs. DFTD SNPs were removed under four conditions: 1) top 5% lowest harmonic mean p-value based on ATOMM single-genome testing (SI Appendix, Fig. S4A), 2) top 5% lowest harmonic mean p-value based on ATOMM interaction testing (SI Appendix, Fig. S4B), 3) both the top 5% lowest p-value based on ATOMM marginal testing and the top 5% lowest p-value based on ATOMM interaction testing (9.9% of SNPs were removed due to overlap between the single-genome test and interaction test SNPs), and 4) randomly removing 10% of DFTD SNPs to ensure any change to PVE above was not simply due to a reduction in the number of SNPs included in the model. After removal of the DFTD SNPs, PVE estimates for devils, DFTD, and devil–DFTD intergenomic interactions were obtained by rerunning ATOMM as described above.

To quantify genetic dissimilarity between DFTD samples of the same lineage, we used the DFTD haploid matrix input to ATOMM, which consisted of 12,420 SNPs across 312 DFTD samples. We first calculated how often tumors differed in their allele at each SNP within a lineage, corresponding to a per-lineage MAF. We also ensured that within-lineage minor alleles were unlikely to be due to sequencing error by extracting genotype quality scores for all minor alleles (see SI Appendix, Results and Discussion for quality metrics). We next assessed within-lineage variation between all SNPs for each sample, a metric which characterized differences at the sample-level rather than the metric characterizing differences at the SNP-level described above. Here, each sample received an individual score based on 1) selecting the sample that it had the most alleles in common with and 2) counting how many alleles differed between the two samples (see the caption of SI Appendix, Fig. S8 for further details).

Identifying Devil–DFTD Interaction Hotspots.

To find regions in the devil and DFTD genome with a high number of interacting SNPs (i.e., “interaction hotspots”), the top 0.1% interaction tests based on the harmonic mean p-value (252,000 interactions; 5.53e-06 ≤ P ≤ 1.30e-02) were selected. Chromosome and position were obtained for devil and DFTD SNPs of the top interaction tests, and the top interacting SNPs were counted within 1-Mb sliding windows separately for devils and DFTD. For devil interactions, a single SNP per linkage group was used when counting top interacting SNPs within each window. However, for both devils and DFTD, a single SNP in one lineage interacting with many SNPs in the other lineage was counted for each interaction (e.g., a single devil SNP which interacted with 100 DFTD SNPs was counted 100 times). Hotspots were determined by taking the windows with the top 1% SNP counts (32 hotspots) for devils (SI Appendix, Fig. S9A) and DFTD (SI Appendix, Fig. S9B), and the reference allele and alternate allele of all SNPs falling within these hotspots were selected. To ensure that SNP density did not predict interaction density, separate linear regressions were fit as interaction density ~ SNP density by counting interacting SNPs (described above) and ATOMM input SNPs in 1 Mb sliding windows. For devils, all ATOMM input SNPs within a linkage group were counted when calculating SNP density. The regressions were fit separately with SNP density as an explanatory variable against the following four response variables: 1) interaction count in devil nonhotspot windows, 2) interaction count in devil hotspot windows, 3) interaction count in DFTD nonhotspot windows, and 4) interaction count in DFTD hotspot windows. For devils, neither nonhotspot SNP density (P = 0.758) nor hotspot SNP density (P = 0.335; SI Appendix, Fig. S9C) predicted interaction density. For DFTD, nonhotspot SNP density predicted interaction density (R2 = 0.330; P < 0.0001) but hotspot SNP density did not (P = 0.946; SI Appendix, Fig. S9D), indicating that SNP density did not predict interaction density within DFTD hotspots. To determine whether LD was elevated in hotspots, and thus a potential confounder to hotspot identification, pairwise R2 was calculated between all SNPs within each 1 Mb sliding window separately for devils and DFTD (only windows with >1 SNP were included). The mean R2 was then obtained for each window, and R2 distributions were plotted separately for hotspot and nonhotspot windows for devils and DFTD (SI Appendix, Fig. S10). Because the R2 distributions were not normally distributed (SI Appendix, Fig. S10), a Mann–Whitney U test was used to determine whether the mean of the hotspot R2 distribution differed from the mean of the nonhotspot R2 distribution.

SNP gene Ontology Enrichment Analysis.

To determine whether devil–DFTD interacting sites were significantly enriched for particular pathways, a gene ontology (GO) enrichment analysis was run using SNP2GO (97) version 1.0.6, setting all ATOMM input SNPs falling within the 32 interaction hotspots as candidate and all other ATOMM input SNPs as noncandidate (a single SNP was used per devil linkage group). Settings were otherwise as follows: FDR to 0.05, runs to 100,000, and extension to 0. The enriched GO terms and their corresponding loci were then obtained and plotted using circlize (98) version 0.4.15 within R.

Characterizing Devil–DFTD Cancer Gene Interactions.

To determine whether the significant interactions identified by ATOMM were enriched for cancer-related genes, the top 0.001% interactions based on the harmonic mean p-value (2,524 interactions; P-value range: 5.53e-06 ≤ P ≤ 6.51e-04) were selected and further characterized using Ensembl’s variant effect predictor [VEP (99)]. To identify the candidate locus of a single SNP in each devil linkage group, the locus of the SNP yielding the lowest interaction P-value for a given linkage group was used. Next, VEP was run separately for devils and DFTD using the mSarHar1.11 reference genome (25) (database version 106), and the arguments “–symbol” and “–pick” (all other settings were default). For devils, the arguments “–distance 100000” (to account for devil LD) and “–pick_order rank” were also used. Only annotated genes were retained because cancer-implication could not be inferred for genes lacking an annotation. See SI Appendix, Materials and Methods for details.

The final list of genic interactions for devils and DFTD was then annotated with known gene–disease associations (GDAs) using disgenet2r (100). Cancer-associated genes were identified by searching for a set of keywords in the disease class name, semantic type, and name fields output by DisGeNET (57). Next, the count of cancer gene x cancer gene interactions (i.e., both host gene and tumor gene had a cancer-related GDA) was calculated for the top 0.001% ATOMM interactions. To determine whether the number of cancer x cancer interactions was significantly greater than expected, Monte-Carlo simulations were performed with 100,000 iterations. Because the ATOMM input SNPs were a subset of all SNPs in the probe-capture targets, the enrichment for cancer-implicated genes in the probe-capture targets may be present in the ATOMM input SNPs, biasing these SNPs (and hence genotype-by-genotype interactions) toward cancer genes. Thus, only genes containing a SNP from the devil and DFTD ATOMM input SNPs were used to construct the simulated dataset. By sampling only from genes containing an ATOMM input SNP, the null distribution created through Monte-Carlo simulations captures any bias toward cancer genes in the initial ATOMM input SNP set. For each iteration, the genes containing an ATOMM input SNP were randomly sampled 96 times (with replacement) for devils and DFTD separately and then paired to simulate the number of annotated gene x gene interactions found in the top 0.001% ATOMM interactions. Next, a count of the cancer x cancer gene interactions (annotated using DisGeNET) out of the 96 interactions was obtained. The final P-value was calculated as the proportion of simulations with a cancer x cancer interaction count greater than or equal to the observed cancer x cancer interaction count (i.e., 88 interactions).

Signatures of Selection in Devil–DFTD Genotype-by-genotype Interactions.

To detect SNPs that exhibited signatures of selection, a genome scan was performed using pcadapt (58). Separate scans were run on the devil and DFTD SNPs used as input to ATOMM (see above), where population structure recovered from pcadapt was equivalent to DAPC population structure (SI Appendix, Fig. S6 A and C). Next, a p-value was calculated for all SNPs corresponding to the probability of a SNP being under selection, and FDR adjustment was used to correct for multiple testing (SI Appendix, Fig. S13 A and C). To determine whether SNPs in the most significant genotype-by-genotype interactions were enriched for selection, the number of SNPs under selection (FDR < 0.05) within devils and DFTD were separately compared between candidates (top 5% lowest interaction p-value within the ATOMM input matrices for devils and DFTD) and noncandidates (all other ATOMM input SNPs). Separate binomial tests (with probability of success as the frequency of significant noncandidates for devils and DFTD) were used to determine whether candidate SNPs were significantly enriched for signatures of selection (SI Appendix, Fig. S13 B and D). Additionally, we used a single binomial test (with probability of success equal to the frequency of devil SNPs significantly under selection) to determine whether the total number of DFTD SNPs under selection was significantly different than the number of devil SNPs under selection (SI Appendix, Fig. S13 A and C). See SI Appendix, Materials and Methods for details.

Genotype–phenotype Association of ANK3 x ARHGAP5 and ANK3 x ACAP1 Devil–DFTD Interactions.

To determine whether there was an association with FOI variation and different categories of interacting alleles for the ANK3 x ARHGAP5 and ANK3 x ACAP1 devil–DFTD interactions, alleles and FOI values were extracted for the paired devil–DFTD samples which were used as input to ATOMM. Each devil and DFTD sample had either a value of “0” (reference allele) or “1” (alternate allele), and each devil–DFTD pair thus existed as one of four interaction categories (shown as devil x DFTD respectively): reference x reference (“0 × 0”), reference x alternate (“0 × 1”), alternate x reference (“1 × 0”), or alternate x alternate (“1 × 1”). DFTD alternate allele quality metrics were also verified (see SI Appendix, Results and Discussion for quality metrics). To account for the possibility of spurious genotype–phenotype associations due to population structure, the spatial distribution of samples within the four interaction categories was determined for the ANK3 x ARHGAP5 and ANK3 x ACAP1 interactions. All interaction categories were found in each major sampling site (Fig. 1B) other than category “0 × 1” and “1 × 1”, which were found exclusively in Freycinet, a single sampling site (FR in Fig. 1B). Because all DFTD samples with an alternate allele for ARHGAP5 (n = 48) and ACAP1 (n = 49) were located exclusively at Freycinet, a comparison of devil–DFTD interaction category versus FOI was made only for Freycinet pairs (Fig. 2D). Comparisons between interaction categories were made with all Freycinet pairs (i.e., all samples shown in Fig. 2D) with an ANOVA and Tukey-HSD post hoc testing. However, because devils in Freycinet did not have a significantly higher FOI relative to the other major sampling sites based on an ANOVA (0.44 ≤ P ≤ 0.99), a comparison was also made which included all devil–DFTD pairs (SI Appendix, Fig. S14).

Differential Expression Analysis of ANK3 in Devils.

To determine whether the devil SNP in the 5′ UTR of ANK3 conferred a change in the expression of ANK3, RNA-seq reads were downloaded from SRA BioProject PRJNA693818 (63). These samples were present in our dataset and possessed all three genotypes for the ANK3 SNP in devils (six homozygous reference allele, nine heterozygotes, and five homozygous alternate allele). The quality of the raw reads was assessed using FastQC (79). TrimGalore! (80) was used to trim adapter sequences, ends of reads with a Phred quality score <10, and to remove reads <30 bp in length (63). Next, trimmed reads were aligned to the mSarHar1.11 (25) reference genome using HISAT2 version 2.1.0 (101) with the –dta flag and sorted using Samtools (82). Stringtie version 1.3.4b (102) and the mSarHar1.11 version 108 (25) devil reference annotation were used to assemble transcripts and estimate gene expression in transcripts per million (TPM).

Supplementary Material

Appendix 01 (PDF)

Dataset S01 (XLSX)

We are grateful to The Save the Tasmania Devil Program from the Department of Natural Resources and Environment—Tasmania for facilitating data and sample collection. The study was funded by the NSF DEB 2027446 to M.J.M., A.S., M.E.J., and H.I.M. and DEB 2324456 to M.J.M., A.S., and R.H., the Sarah and Daniel Hrdy Visiting Fellowship in Conservation Biology at Harvard University to M.J.M., the Australian Research Council (ARC) Future Fellowship (FT100100250) and ARC Large Grants (A00000162) to M.E.J., the ARC grants (DE 170101116; LP 170101105), the Agence Nationale de la Recherche (ANR Blanc Project TRANSCAN), and a Centre National de la Recherche Scientifique (CNRS) International Associated Laboratory Grant to R.H., and the University of South Florida to M.J.M. and R.M.

Author contributions

H.I.M., R.H., M.E.J., A.S., and M.J.M. designed research; D.G.G., C.P.K., R.M.R., M.A.B., M.R.-A., S.C., D.G.H., D.H.K., H.I.M., R.H., M.E.J., A.S., R.M., and M.J.M. performed research; D.G.G., R.M.R., D.H.K., and R.M. analyzed data; and D.G.G., R.M., and M.J.M. wrote the paper.

Competing interests

The authors declare no competing interest.

Data, Materials, and Software Availability

Genomic sequencing data are deposited under BioProject PRJNA947341 (https://dataview.ncbi.nlm.nih.gov/object/PRJNA947341?reviewer=abr2clfa1fv6s02ntrndueqd63) and accessions SRR23935061–SRR23936020 (78). Phenotype data are available at https://github.com/D-gallinson/Devil-DFTD-FOI-Coevolution (86). RNA-seq data were previously deposited under BioProject PRJNA693818 accessions SRR13765773, 774, 794, 796, 798–800, 802, 804, 805, 807, 809, 810, 817, 819–824 (https://www.ncbi.nlm.nih.gov/bioproject/?term=PRJNA693818) (63).

Supporting Information

This article is a PNAS Direct Submission.
==== Refs
1 P. Daszak, A. A. Cunningham, A. D. Hyatt, Emerging infectious diseases of wildlife–threats to biodiversity and human health. Science 287 , 443–449 (2000).10642539
2 J. Castiblanco, J.-M. Anaya, Genetics and Vaccines in the Era of Personalized Medicine. Curr. Genomics 16 , 47–59 (2015).25937813
3 V. Sintchenko, E. C. Holmes, The role of pathogen genomics in assessing disease transmission. BMJ 350 , h1314 (2015).25964672
4 O. Anacleto , Genetic differences in host infectivity affect disease spread and survival in epidemics. Sci. Rep. 9 , 4924 (2019).30894567
5 E. K. Karlsson, D. P. Kwiatkowski, P. C. Sabeti, Natural selection and infectious disease in human populations. Nat. Rev. Genet. 15 , 379–393 (2014).24776769
6 C. P. Kozakiewicz , Pathogens in space: Advancing understanding of pathogen dynamics and disease ecology through landscape genetics. Evol. Appl. 11 , 1763–1778 (2018).30459828
7 M. J. Newport, C. Finan, Genome-wide association studies and susceptibility to infectious diseases. Brief. Funct. Genomics 10 , 98–107 (2011).21436306
8 J. P. Andras, P. D. Fields, L. Du Pasquier, M. Fredericksen, D. Ebert, Genome-wide association analysis identifies a genetic basis of infectivity in a model bacterial pathogen. Mol. Biol. Evol. 37 , 3439–3452 (2020).32658956
9 A. Gori , Pan-GWAS of Streptococcus agalactiae highlights lineage-specific genes associated with virulence and niche adaptation. mBio 11 , e00728-20 (2020).32518186
10 T. A. Manolio, Genomewide association studies and assessment of the risk of disease. N Engl. J. Med. 363 , 166–176 (2010).20647212
11 V. Tam , Benefits and limitations of genome-wide association studies. Nat. Rev. Genet. 20 , 467–484 (2019).31068683
12 E. Uffelmann , Genome-wide association studies. Nat. Rev. Methods Primers 1 , 1–21 (2021).
13 S. I. Rothstein, A model system for coevolution: Avian brood parasitism. Ann. Rev. Ecol. Syst. 21 , 481–508 (1990).
14 T. G. Whitham , Community and ecosystem genetics: A consequence of the extended phenotype. Ecology 84 , 559–573 (2003).
15 A. Shimizu , Fine-tuned bee-flower coevolutionary state hidden within multiple pollination interactions. Sci. Rep. 4 , 3988 (2014).24496444
16 C. Rafaluk , Rapid evolution of virulence leading to host extinction under host-parasite coevolution. BMC Evol. Biol. 15 , 112 (2015).26070343
17 C. Bonneaud , Rapid antagonistic coevolution in an emerging pathogen and its vertebrate host. Curr. Biol. 28 , 2978–2983.e5 (2018).30197084
18 P. H. Thrall , Rapid genetic change underpins antagonistic coevolution in a natural host-pathogen metapopulation. Ecol. Lett. 15 , 425–435 (2012).22372578
19 M. A. Ansari , Genome-to-genome analysis highlights the effect of the human innate and adaptive immune systems on the hepatitis C virus. Nat. Genet. 49 , 666–673 (2017).28394351
20 M. Wang , Two-way mixed-effects methods for joint association analysis using both host and pathogen genomes. Proc. Natl. Acad. Sci. U.S.A. 115 , E5440–E5449 (2018).29848634
21 A. MacPherson, S. P. Otto, S. L. Nuismer, Keeping pace with the red queen: Identifying the genetic basis of susceptibility to infectious disease. Genetics 208 , 779–789 (2018).29223971
22 H. Märkle, A. Tellier, S. John, Cross-Species association statistics for genome-wide studies of host and parasite polymorphism data. BioRxiv [Preprint] (2021). 10.1101/72616.
23 M. Kardos, A. Husby, S. E. McFarlane, A. Qvarnström, H. Ellegren, Whole-genome resequencing of extreme phenotypes in collared flycatchers highlights the difficulty of detecting quantitative trait loci in natural populations. Mol. Ecol. Resour. 16 , 727–741 (2016).26649993
24 A. H. Patton , A transmissible cancer shifts from emergence to endemism in Tasmanian devils. Science 370 , eabb9772 (2020).33303589
25 M. R. Stammnitz , The evolution of two transmissible cancers in Tasmanian devils. Science 380 , 283–293 (2023).37079675
26 C. E. Hawkins , Emerging disease and population decline of an island endemic, the Tasmanian devil Sarcophilus harrisii. Biol. Conserv. 131 , 307–324 (2006).
27 C. X. Cunningham , Quantifying 25 years of disease-caused declines in Tasmanian devil populations: Host density drives spatial pathogen spread. Ecol. Lett. 24 , 958–969 (2021).33638597
28 A.-M. Pearse, K. Swift, Transmission of devil facial-tumour disease. Nature 439 , 549–549 (2006).16452970
29 H. McCallum , Transmission dynamics of Tasmanian devil facial tumor disease may lead to disease-induced extinction. Ecology 90 , 3379–3392 (2009).20120807
30 S. Lachish, H. McCallum, M. Jones, Demography, disease and the devil: Life-history changes in a disease-affected population of Tasmanian devils (Sarcophilus harrisii). J. Anim. Ecol. 78 , 427–436 (2009).19021786
31 K. Wells , Individual and temporal variation in pathogen load predicts long-term impacts of an emerging infectious disease. Ecology 100 , e02613 (2019).30636287
32 H. McCallum , Distribution and impacts of Tasmanian devil facial tumor disease. EcoHealth 4 , 318 (2007).
33 R. K. Hamede , Transmissible cancer in Tasmanian devils: Localized lineage replacement and host population response. Proc. Biol. Sci. 282 , 20151468 (2015).26336167
34 B. Epstein , Rapid evolutionary response to a transmissible cancer in Tasmanian devils. Nat. Commun. 7 , 12684 (2016).27575253
35 B. T. Lazenby , Density trends and demographic signals uncover the long-term impact of transmissible cancer in Tasmanian devils. J. Appl. Ecol. 55 , 1368–1379 (2018).30089931
36 M. J. Margres , Large-effect loci affect survival in Tasmanian devils (Sarcophilus harrisii) infected with a transmissible cancer. Mol. Ecol. 27 , 4189–4199 (2018).30171778
37 H. V. Siddle , Transmission of a fatal clonal tumor by biting occurs due to depleted MHC diversity in a threatened carnivorous marsupial. Proc. Natl. Acad. Sci. U.S.A. 104 , 16221–16226 (2007).17911263
38 K. Wells , Infection of the fittest: Devil facial tumour disease has greatest effect on individuals with highest reproductive output. Ecol. Lett. 20 , 770–778 (2017).28489304
39 M. J. Margres , The genomic basis of tumor regression in Tasmanian devils (Sarcophilus harrisii). Genome Biol. Evol. 10 , 3012–3025 (2018), 10.1093/gbe/evy229.30321343
40 Y. M. Kwon , Evolution and lineage dynamics of a transmissible cancer in Tasmanian devils. PLoS Biol. 18 , e3000926(2020).33232318
41 R. Dawkins, The Extended Phenotype (Oxford University Press, 1982).
42 M. E. Jones, D. Paetkau, E. Geffen, C. Moritz, Genetic diversity and population structure of Tasmanian devils, the largest marsupial carnivore. Mol. Ecol. 13 , 2197–2209 (2004).15245394
43 W. Miller , Genetic diversity and population structure of the endangered marsupial Sarcophilus harrisii (Tasmanian devil). Proc. Natl. Acad. Sci. U.S.A. 108 , 12348–12353 (2011).21709235
44 S. Hendricks , Conservation implications of limited genetic diversity and population structure in Tasmanian devils (Sarcophilus harrisii). Conserv. Genet. 18 , 977–982 (2017).28966567
45 X. Zhou, P. Carbonetto, M. Stephens, Polygenic modeling with Bayesian sparse linear mixed models. PLoS Genet. 9 , e1003264 (2013).23408905
46 J. M. Alves , Parallel adaptation of rabbit populations to myxoma virus. Science 363 , 1319–1326 (2019).30765607
47 D. R. Rokyta, P. Joyce, S. B. Caudle, H. A. Wichman, An empirical test of the mutational landscape model of adaptation using a single-stranded DNA virus. Nat. Genet. 37 , 441–444 (2005).15778707
48 M. Ní Leathlobhair, R. E. Lenski, Population genetics of clonally transmissible cancers. Nat. Ecol. Evol. 6 , 1077–1089 (2022).35879542
49 L. Evenepoel , Expression of contactin 4 is associated with malignant behavior in pheochromocytomas and paragangliomas. J. Clin. Endocrinol. Metab. 103 , 46–55 (2018).28938490
50 D. Huang , GLI2 promotes cell proliferation and migration through transcriptional activation of ARHGEF16 in human glioma cells. J. Exp. Clin. Cancer Res 37 , 247 (2018).30305138
51 D. Falush, R. Bowden, Genome-wide association mapping in bacteria? Trends Microbiol. 14 , 353–355 (2006).16782339
52 T. D. Read, R. C. Massey, Characterizing the genetic basis of bacterial phenotypes using genome-wide association studies: A new direction for bacteriology. Genome Med. 6 , 109 (2014).25593593
53 P. E. Chen, B. J. Shapiro, The advent of genome-wide association studies for bacteria. Curr. Opin. Microbiol. 25 , 17–24 (2015).25835153
54 J. E. Deakin , Genomic restructuring in the Tasmanian devil facial tumour: Chromosome painting and gene mapping provide clues to evolution of a transmissible tumour. PLoS Genet. 8 , e1002483 (2012).22359511
55 S. Vasseur, F. Guillaumond, Lipids in cancer: A global view of the contribution of lipid pathways to metastatic formation and treatment resistance. Oncogenesis 11 , 1–15 (2022).35013097
56 J. Ma , Emerging roles of nucleotide metabolism in cancer development: Progress and prospect. Aging (Albany NY) 13 , 13349–13358 (2021).33952722
57 J. Piñero , The DisGeNET knowledge platform for disease genomics: 2019 update. Nucleic Acids Res. 48 , D845–D855 (2020).31680165
58 F. Privé, K. Luu, B. J. Vilhjálmsson, M. G. B. Blum, Performing highly efficient genome scans for local adaptation with R Package pcadapt Version 4. Mol. Biol. Evol. 37 , 2153–2154 (2020).32343802
59 E. Kordeli, S. Lambert, V. Bennett, Ankyrin: A new ankyrin gene with neural-specific isoforms localized at the axonal initial segment and node of ranvieR (∗). J. Biol. Chem. 270 , 2352–2359 (1995).7836469
60 T. Wang , Ankyrin G expression is associated with androgen receptor stability, invasiveness, and lethal outcome in prostate cancer patients. J. Mol. Med. 94 , 1411–1422 (2016).27534968
61 G. Dong , SIRT1 suppresses the migration and invasion of gastric cancer by regulating ARHGAP5 expression. Cell Death Dis. 9 , 977 (2018).30250020
62 J. Zhang, Q. Zhang, J. Zhang, Q. Wang, Expression of ACAP1 is associated with tumor immune infiltration and clinical outcome of ovarian cancer. DNA Cell Biol. 39 , 1545–1557 (2020).32456571
63 C. P. Kozakiewicz , Spatial variation in gene expression of Tasmanian devil facial tumors despite minimal host transcriptomic response to infection. BMC Genomics 22 , 698 (2021).34579650
64 J. D. Schiffman, M. Breen, Comparative oncology: What dogs and other species can teach us about humans with cancer. Philos. Trans. R Soc. B Biol. Sci. 370 , 20140231 (2015).
65 F. Riccardo, L. Aurisicchio, J. A. Impellizeri, F. Cavallo, The importance of comparative oncology in translational medicine. Cancer Immunol. Immunother 64 , 137–148 (2015).25548094
66 A. M. Boddy, T. M. Harrison, L. M. Abegglen, Comparative oncology: New insights into an ancient disease. iScience 23 , 101373 (2020).32738614
67 R. M. Anderson, R. M. May, Coevolution of hosts and parasites. Parasitology 85 , 411–426 (1982).6755367
68 R. M. May, R. M. Anderson, Epidemiology and genetics in the coevolution of parasites and hosts. Proc. R. Soc. Lond. B Biol. Sci. 219 , 281–313 (1983).6139816
69 M. A. Brockhurst , Running with the Red Queen: The role of biotic conflicts in evolution. Proc. R Soc. B Biol. Sci. 281 , 20141382 (2014).
70 R. Hamede , Reduced effect of Tasmanian devil facial tumor disease at the disease front. Conserv. Biol. 26 , 124–134 (2012).21978020
71 S. James , Tracing the rise of malignant cell lines: Distribution, epidemiology and evolutionary interactions of two transmissible cancers in Tasmanian devils. Evol. Applications 12 , 1772–1780 (2019).
72 R. Loh , The pathology of devil facial tumor disease (DFTD) in Tasmanian Devils (Sarcophilus harrisii). Vet. Pathol. 43 , 890–895 (2006).17099145
73 A. K. Fraik , Transcriptomics of Tasmanian devil (Sarcophilus harrisii) ear tissue reveals homogeneous gene expression patterns across a heterogeneous landscape. Genes 10 , 801 (2019).31614864
74 E. P. Murchison , Genome sequencing and analysis of the tasmanian devil and its transmissible cancer. Cell 148 , 780–791 (2012).22341448
75 B. Wright , Variants in the host genome may inhibit tumour growth in devil facial tumours: Evidence from genome-wide association. Sci. Rep. 7 , 423 (2017).28341828
76 M. R. Stammnitz , The origins and vulnerabilities of two transmissible cancers in Tasmanian devils. Cancer Cell 33 , 607–619.e15 (2018).29634948
77 J. N. Weinstein , The Cancer Genome Atlas Pan-Cancer Analysis Project. Nat. Genet. 45 , 1113–1120 (2013).24071849
78 D. G. Gallinson, C. P. Kozakiewicz, M. A. Beer, A. Storfer, M. J. Margres. Tasmanian devil and devil facial tumor disease probe-capture sequencing. Sequence Read Archive (SRA). https://dataview.ncbi.nlm.nih.gov/object/PRJNA947341?reviewer=abr2clfa1fv6s02ntrndueqd63. Deposited 16 January 2023.
79 S. Andrews, FastQC: A Quality Control tool for High Throughput Sequence Data (2010) (16 November 2022).
80 F. Krueger, F. James, P. Ewels, E. Afyounian, B. Schuster-Boeckler, FelixKrueger/TrimGalore: v0.6.7–DOI via Zenodo (2021). 10.5281/zenodo.5127899. 24 November 2022.
81 H. Li, Aligning sequence reads, clone sequences and assembly contigs with BWA-MEM. arXiv:1303.3997 (2013) (19 February 2022).
82 P. Danecek , Twelve years of SAMtools and BCFtools. Gigascience 10 , giab008 (2021).33590861
83 Broad Institute, Picard toolkit. GitHub repository. https://broadinstitute.github.io/picard/. Accessed 19 February 2022.
84 R. Poplin , Scaling accurate genetic variant discovery to tens of thousands of samples. BioRxiv [Preprint] (2018). 10.1101/201178 (Accessed 24 November 2022).
85 H. Muench, Derivation of rates from summation data by the catalytic curve. J. Am. Stat. Associat. 29 , 25–38 (1934).
86 D. G. Gallinson, Devil-DFTD-FOI-Coevolution. Github. https://github.com/D-gallinson/Devil-DFTD-FOI-Coevolution. Deposited 1 April 2023.
87 X. Zhou, M. Stephens, Genome-wide efficient mixed-model analysis for association studies. Nat. Genet. 44 , 821–824 (2012).22706312
88 A. Vehtari, A. Gelman, D. Simpson, B. Carpenter, P.-C. Bürkner, Rank-normalization, folding, and localization: An Improved Rˆ for Assessing Convergence of MCMC (with Discussion). Bayesian Analysis 16 , 667–718 (2021).
89 Stan Development Team, RStan: the R interface to Stan, R package version 2.21.2. https://mc-stan.org/. Accessed 26 November 2022.
90 R Core Team, R: A language and environment for statistical computing (2022). 8 December 2022.
91 Z. R. McCaw, J. M. Lane, R. Saxena, S. Redline, X. Lin, Operating characteristics of the rank-based inverse normal transformation for quantitative trait analysis in genome-wide association studies. Biometrics 76 , 1262–1272 (2020).31883270
92 D. J. Wilson, The harmonic mean p-value for combining dependent tests. Proc. Natl. Acad. Sci. U.S.A. 116 , 1195–1200 (2019).30610179
93 T. Jombart, S. Devillard, F. Balloux, Discriminant analysis of principal components: A new method for the analysis of genetically structured populations. BMC Genet. 11 , 94 (2010).20950446
94 J. A. Lees , Genome-wide identification of lineage and locus specific variation associated with pneumococcal carriage duration. eLife 6 , e26255 (2017).28742023
95 D. A. Smith , Viral genome wide association study identifies novel hepatitis C virus polymorphisms associated with sofosbuvir treatment failure. Nat. Commun. 12 , 6105 (2021).34671027
96 G. Hahn , Genome-wide association analysis of COVID-19 mortality risk in SARS-CoV-2 genomes identifies mutation in the SARS-CoV-2 spike protein that colocalizes with P.1 of the Brazilian strain. Genetic Epidemiol. 45 , 685–693 (2021).
97 D. Szkiba, M. Kapun, A. von Haeseler, M. Gallach, SNP2GO: Functional analysis of genome-wide association studies. Genetics 197 , 285–289 (2014).24561481
98 Z. Gu, L. Gu, E. Roland, M. Schlesner, B. Brors, circlize implements and enhances circular visualization in R (2014) (November 28, 2022).
99 W. McLaren , The Ensembl Variant Effect Predictor. Genome Biology 17 , 122 (2016).27268795
100 IBI Group, disgenet2r: An R package to query, visualize, and expand DisGeNET data (2021).
101 D. Kim, J. M. Paggi, C. Park, C. Bennett, S. L. Salzberg, Graph-based genome alignment and genotyping with HISAT2 and HISAT-genotype. Nat. Biotechnol. 37 , 907–915 (2019).31375807
102 S. Kovaka , Transcriptome assembly from long-read RNA-seq alignments with StringTie2. Genome Biol. 20 , 278 (2019).31842956
