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.
Transmissible 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)(11)(12), ignoring joint or extended phenotypes (e.g., refs. [13][14][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)(17)(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)(19)(20)(21)(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 pnas.org than devil extinction (31), suggesting coevolutionary interactions will be key drivers of disease dynamics and thus important for long-term disease management.
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)(33)(34)(35)(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)(43)(44) and/or distinct DFTD lineages (24,40). (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.
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.
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 pnas.org 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 genotypeby-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)(52)(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 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.
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 singlegenome 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.
and Pathogen. The most significant devil-DFTD genotypeby-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 andB). 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 andB) or DFTD (P = 0.091; SI Appendix, Fig. S10 C andD). 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-proteinlinked 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 pnas.org 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 andD).
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 andB), 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).
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.
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.
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)(75)(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 cm 3 (95% CI = 198-223 cm 3 ) 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 cm 3 , 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 cm 3 (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 largeeffect 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 95 th 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 hostpathogen 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.
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. 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 withinlineage 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 (R 2 = 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 R 2 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 R 2 was then obtained for each window, and R 2 distributions were plotted separately for hotspot and nonhotspot windows for devils and DFTD (SI Appendix, Fig. S10). Because the R 2 distributions were not normally distributed (SI Appendix, Fig. S10), a Mann-Whitney U test was used to determine whether the mean of the hotspot R 2 distribution differed from the mean of the nonhotspot R 2 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). Cancerassociated 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).
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 andC). 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 andC). 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 andD). 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
Downloaded from https://www.pnas.org by "UNIVERSITY OF SOUTH FLORIDA, LIB SERIALS/ACQUISITIONS" on March 12, 2024 from IP address 97.97.52.6.