Whole genome resequencing reveals signatures of rapid selection in a virus-affected commercial fishery.

Infectious diseases are recognized as one of the greatest global threats to biodiversity and ecosystem functioning. Consequently, there is a growing urgency to understand the speed at which adaptive phenotypes can evolve and spread in natural populations to inform future management. Here we provide evidence of rapid genomic changes in wild Australian blacklip abalone (Haliotis rubra) following a major population crash associated with an infectious disease. Genome scans on H. rubra were performed using pooled whole genome resequencing data from commercial fishing stocks varying in historical exposure to haliotid herpesvirus-1 (HaHV-1). Approximately 25,000 single nucleotide polymorphism loci associated with virus exposure were identified, many of which mapped to genes known to contribute to HaHV-1 immunity in the New Zealand pāua (Haliotis iris) and herpesvirus response pathways in haliotids and other animal systems. These findings indicate genetic changes across a single generation in H. rubra fishing stocks decimated by HaHV-1, with stock recovery potentially determined by rapid evolutionary changes leading to virus resistance. This is a novel example of apparently rapid adaptation in natural populations of a nonmodel marine organism, highlighting the pace at which selection can potentially act to counter disease in wildlife communities.

INTRODUCTION

The spread of infectious diseases is recognized as one of the most pressing global threats to biodiversity and ecosystem function (Cunningham et al., 2017; Daszak et al., 2000; Tompkins et al., 2015). In recent decades, infectious diseases have devastated a range of wildlife groups (Berger et al., 1998; Hansen et al., 2005; Kim & Harvell, 2004; Lorch et al., 2016), often exacerbating species declines in ecosystems already stressed by climate change and habitat destruction (Bosch et al., 2018; Brearley et al., 2013; Harvell et al., 2002). The persistence of many species will probably depend on their ability to adapt to environmental changes associated with increased disease prevalence, although selection for disease resistance or tolerance may not keep pace with rates of pathogen evolution and the emergence and turnover of novel diseases (Hawley et al., 2013; Ujvari et al., 2014).

Detecting rapid evolutionary changes in populations impacted by environmental disturbance is often challenging. Disease‐affected populations provide ideal models given disease exposure is an easily characterized selective pressure and studies of this nature have been greatly assisted by modern genomic technologies (Blanchong et al., 2016; Storfer et al., 2020). These technologies now allow for rapid and cost‐effective estimates of genome‐wide variation among populations spanning disease infection gradients and individuals with distinctive phenotypes related to disease response (Auteri & Knowles, 2020; Elbers et al., 2018; Grogan et al., 2018; Margres et al., 2018). Importantly, a number of studies using these technologies have reported rapid evolutionary changes across several generations in natural populations of nonmodel organisms impacted by disease, including Tasmanian devils (Sarcophilus harrisii) (Epstein et al., 2016; Hubert et al., 2018; Margres et al., 2018) and North American house finches (Carpodacus mexicanus) (Bonneaud et al., 2011). Additionally, recent studies have reported evidence of rapid selection for disease‐resistant genotypes across a single generation in North American sea stars (Pisaster ochraceus) and little brown bats (Myotis lucifugus), following rapid and severe population crashes due to infectious diseases (Auteri & Knowles, 2020; Schiebelhut et al., 2018). Such studies are pivotal in highlighting the pace at which selection can act to counter disease in wildlife communities and opening up opportunities for interventions, such as deliberate translocations of adaptive phenotypes, that can increase the adaptability of threatened populations (Hoffmann et al., 2020; Hohenlohe et al., 2019). Despite this progress, the number of studies demonstrating rapid evolutionary responses to infectious diseases in natural populations remains limited and biased towards terrestrial systems.

Marine infectious diseases are responsible for incremental and mass mortalities in a variety of wildlife groups, including keystone and habitat‐forming taxa (Clemente et al., 2014; Harvell & Lamb, 2020; Harvell et al., 2007; Martin et al., 2016; Montecino‐Latorre et al., 2016), and species supporting wild commercial fisheries (Cawthorn, 2011; Crosson et al., 2020; Lafferty et al., 2015; Marty et al., 2010). The Australian blacklip abalone (Haliotis rubra), a species targeted by the world’s largest wild abalone fisheries and a rapidly expanding aquaculture industry (FAO FishStat, 2021), was heavily impacted by disease between 2006 and 2010 (Mayfield et al., 2012). Beginning in 2006, abalone viral ganglioneuritis (AVG) caused by the haliotid herpesvirus‐1 (HaHV‐1) spread along the western coastline of Victoria in southeastern Australia, causing rapid and severe population collapses (>90% mortality in some areas) and devastating both wild and farmed abalone stocks (Hooper et al., 2007). Despite the impact of AVG, abalone stocks in the Western Zone fishery have seen significant recovery (Western Abalone Divers Association, 2020). It is possible that rapid evolutionary responses to the virus have contributed to this recovery, facilitated by the abalone’s short generation time (~4 years; Andrews, 1999), large population sizes (Mayfield et al., 2012) and high genetic variability that contributes to existing patterns of adaptation across the fishery (Miller et al., 2019).

Evolving immunity to HaHV‐1 depends on the availability of specific, heritable genetic variants within abalone populations that lead to resistant phenotypes. Indeed, previous research has demonstrated heritable genetic variation relating to herpesvirus immunity in sister Haliotid species, highlighting the possibility of evolved immunity in H. rubra. Challenge tests performed on New Zealand paua (Haliotis iris) and Japanese black abalone (Haliotis discus), involving controlled exposure to HaHV‐1, indicated complete immunity to AVG (Chang et al., 2005; Corbeil et al., 2017), with complementary transcriptomic analyses helping to characterize the genetic basis of the resistance (Bai, Zhang, et al., 2019; Neave et al., 2019). Similar tests on H. rubra yielded no evidence of resistance to AVG (Corbeil et al., 2016; Crane et al., 2013); however, these experiments were performed on animals from a limited number of locations affected by AVG. While complete immunity may not occur in H. rubra, the presence of AVG immunity in sister taxa hints at the potential for some level of resistance developing through standing genetic variation following AVG exposure.

For the first time in a decade AVG was recorded in 2021, leading to abalone mortalities at a few proximal fishing locations heavily impacted by AVG in the early 2000s (Agriculture Victoria, 2021). Unlike the first outbreak, animal mortality and disease spread has been minimal. While environmental and epidemiological factors may be contributing to the suppression of the disease (Bai, Li, et al., 2019; Corbeil, 2020), it is also possible that the mass mortality event from 2006–2010 selected for adaptive phenotypes which is reducing the number of susceptible animals and overall viral load within affected fishing stocks. To test this hypothesis, we investigate potential signatures of evolutionary changes in recovering H. rubra fishing stocks devastated by AVG. Specifically, we performed genome scans using pooled whole genome resequencing data on H. rubra specimens from fishing stocks varying in disease exposure. Our findings point to rapid changes in population‐level allele frequencies over a single generation timescale in virus‐affected fishing stocks, with stock recovery potentially determined by rapid evolutionary changes leading to virus resistance. This study highlights the pace at which adaptive phenotypes can potentially evolve and spread in wildlife communities to counter threats from infectious diseases. We discuss these findings in the context of future biosecurity management of Australian abalone fisheries and wildlife conservation more generally.

MATERIALS AND METHODS

Sample collection and DNA sequencing

Tissue biopsies were collected from 343 individual Haliotis rubra from 14 locations spanning the western and central coastline of Victoria. Locations were selected based on their known virus exposure history according to records held by the Victorian wild fishing sector and the Victorian Fisheries Authority. During the outbreak event, rigorous PCR (polymerase chain reaction) testing and diver surveys were conducted across the fishery to determine the spatial extent of the virus and rates of local mortality at affected sites, estimated from the abundance of dead and dying abalone and recently vacated shells (Gorfine et al., 2008). Based on these data we selected 10 AVG‐affected locations that showed >70% animal mortalities, and four AVG‐unaffected locations (Table 1, Figure 1). Sampling for five of the locations was coordinated in 2009 by the Department of Economic Development, Jobs, Transport and Resources (DEDJTR). It is expected that animals from these locations during this sampling period survived the virus event and represent the genomic variation preserved post‐AVG. Sampling of the remaining nine locations was performed between 2015 and 2020. To avoid the potential swamping effects of intergenerational gene flow since the disease outbreak, sampling was biased towards fishing stocks expected to be largely self‐recruiting based on biophysical connectivity models (Young et al., 2020), and toward large adult animals (expected to be either direct survivors or first‐generation post‐virus survivors). However, it is important to note that these fishing stocks are not completely self‐recruiting entities given previous studies have demonstrated a lack of genetic structure across the fishery (Miller et al., 2016, 2019). This sampling was performed by contract divers, commercial fishermen and our research team. At each location, individual abalone were collected within a 100‐m2 area, with tissue biopsies consisting of 20 mg of muscle tissue from the abalone lip obtained using sterile dissection tools to avoid cross‐contamination. Biopsied material was transferred to 2‐ml microcentrifuge tubes containing 80%–100% ethanol and stored at 4°C until required for genomic analysis.

TABLE 1

Site location details and corresponding codes for 14 collection locations of Haliotis rubra used for genomic analyses. Sample sizes and AVG exposure history are also provided

Zone and location	Code	Year sampled	Sample size	GPS location		AVG status
				Latitude	Longitude
Port Macdonell	PMC	2020	25	–38.054	140.881	Unaffected
Inside Murrels	MUR	2009	25	–38.407	141.524	Affected
Inside Nelson	ISN	2009	25	–38.409	141.558	Affected
Lady Julia Percy	LJP	2009	23	–38.422	141.993	Unaffected
The Crags	CRG	2009	25	–38.390	142.135	Affected
Killarney	KIL	2015	20	–38.363	142.321	Affected
Levies	LEV	2009	25	–38.385	142.235	Affected
Childers Cove	CHC	2019	25	–38.490	142.672	Affected
Bay of Islands	BIP	2019	25	–38.582	142.827	Affected
Cat Reef	CAT	2015	25	–38.741	143.188	Affected
White Cliffs	WCF	2015	25	–38.758	143.330	Affected
Castle Cove	CCV	2020	25	–38.783	143.422	Unaffected
Parker River	PKR	2020	25	–38.855	143.538	Affected
Blanket Bay	BLK	2015	25	–38.827	143.586	Unaffected

FIGURE 1

Sampling sites selected for population genomic analysis from southeastern Australia. Figure legend and colour coding of mapped sites indicate history of virus exposure. Refer to Table 1 for sample codes

Total genomic DNA was extracted from 10 mg of tissue using a DNeasy Blood and Tissue Kit (Qiagen) following the manufacturer’s instructions. Resulting DNA extracts were quantified using a Qubit version 2 fluorometer (Life Technologies). To obtain population genomic data, we applied the Pool‐Seq approach (Futschik & Schlotterer, 2010), which involves pooling the DNA of a large number of individuals from the same population and then sequencing the “population variability genome.” This was achieved by pooling equimolar amounts of individual DNA extracts from each sample location, splitting the 25 individuals per location into two pools per location consisting of DNA from 12 and 13 individuals, respectively, to account for potential sequencing bias. The resulting 28 pooled libraries were prepared for sequencing using the Nextera DNA Sample Preparation kit, pooled into a single Illumina NovaSeq S4 flowcell (Illumina) and sequenced across all four lanes with the 150‐bp paired‐end protocol. Sequencing was performed allowing for 3× genome coverage per individual per pool, equating to 80–100× genome coverage per population.

Data preparation

The Illumina NovaSeq sequencing yielded a total of 25 × 109 assigned 150‐bp reads, and a total of 45–100 Gb of sequence data for each of the 28 pooled DNA libraries. Raw DNA sequence reads from the two separate pooled libraries per sample location were pooled for processing purposes. Raw sequences were processed using trimmomatic version 0.36 (Bolger et al., 2014) by removing Nextera adaptors and discarding all reads that had a Phred score below 20. All retained reads were subsequently aligned to the H. rubra reference genome (NCBI RefSeq QXJH00000000.1; Gan et al., 2019) using the ppalign package in the PoolParty pipeline (Micheletti & Narum, 2018) with default parameters. Single nucleotide polymorphisms (SNPs) were called using poolfstat (Hivert et al., 2018) where sites were required to have a read depth of 40–200 reads to be called. SNPs with a minor allele frequency of ≥0.05 were used for downstream genomic analysis.

| Estimating overall genetic structure

SNP frequencies over all loci were initially contrasted between all 14 sample locations to determine patterns of overall genetic structure and population connectivity. The software poolfstat implemented in R (Hivert et al., 2018) was used to calculate global and pairwise measures of population differentiation (F ST; Weir & Cockerham, 1984). Estimates of global and pairwise F ST were also generated using a data set consisting of only candidate loci identified from association analyses.

| Genotype × environment association analysis

To identify SNPs associated with virus exposure status, we performed a genotype × environment association analyses using baypass 2.1 (Gautier, 2015). Analyses were performed under the auxiliary (AUX) covariate mode (‐covmcmc and ‐auxmode flags), after scaling the variables with the ‐scalecov flag, using input files containing genotypic data for each sampling location and corresponding virus exposure history (exposed vs. unexposed). The underlying models explicitly account for the covariance structure among the population allele frequencies that originates from the shared history of the populations through estimation of the population covariance matrix Ω, which removes the variation associated with demography (Bonhomme et al., 2010; Günther & Coop, 2013). The auxiliary covariate model specifically involves the introduction of a binary auxiliary variable to classify each locus as associated or not associated. This allows computation of posterior inclusion probabilities (and Bayes Factors) for each locus while explicitly accounting for multiple testing issues. The auxiliary covariate model was applied with default parameters, a 5000 burn‐in of iterations in the Markov chain Monte Carlo (MCMC) simulation, followed by 25,000 iterations. To reduce artefacts due to potential variability between runs, we performed three independent baypass simulations. We then calculated the average Bayes Factor (BF), expressed in deciban units (dB), for each SNP as a quantitative estimate of the strength of association with virus exposure and the standardized allele frequency. For each SNP, the level of effect was assessed based on the BF models according to Jeffrey’s rule (Jeffreys, 1961). SNPs with BF scores ≥50 were regarded as decisive associations with virus exposure and were retained as potential candidate loci.

Post hoc analyses including functional annotations

An analysis of principal components (PCA) was implemented in the adegenet package for R (Jombart, 2008; Jombart & Ahmed, 2011) to obtain a graphical depiction of patterns of genetic structure among virus‐affected and unaffected stocks based on all candidate SNPs identified by baypass (BFs ≥ 50).

Total linkage disequilibrium (LD) among all candidate loci was calculated using ldx, a package which uses an approximate maximum‐likelihood approach from pooled resequencing data (Feder et al., 2012). LD was calculated as r 2, the square of the correlation between alleles of SNP pairs within the paired sequence reads of each population. We subsequently calculated the average LD for each pairwise SNP comparison across sample sites. Next, we assessed the distribution of candidate loci and signatures of selection across the reference H. rubra genome consisting of 2854 annotated scaffolds varying between 1830 and 1.1 × 107 bp in length (Gan et al., 2019). This was achieved by regressing the total number of candidate loci against scaffold length using the ggpubr package for R (Kassambara & Kassambara, 2020). Scaffolds with exceptionally large numbers of candidate loci relative to scaffold length (deviating from a linear distribution) were interrogated further using the package ldblockshow (Dong et al., 2021) to measure pairwise linkage disequilibrium and haplotype blocks using the default ‐SeleVar option to calculate D′ (the ratio of the difference between the observed and expected frequency of a haplotype, and its maximum value when considering total allele frequencies).

snpeff version 2.0.3 (Cingolani et al., 2012) was used to map candidate SNP loci to the H. rubra genome and predict variant impacts: high (highly disruptive impact on protein function), moderate (nonsynonymous mutations, possible change in protein effectiveness), low (unlikely to change protein behaviour) or a modifier (synonymous mutations, noncoding or intergenic variant). Functional classification of candidate genes was achieved by aligning the peptide sequences for mapped candidate H. rubra genes with the annotated genomes for human (NCBI RefSeq IDs NC_000001–NC_000024), pacific oyster (RefSeq IDs NC_047559–NC_047568), scallop (RefSeq ID NC_007234.1) and the blue mussel (RefSeq ID NC_006161.1) using diamond (Buchfink et al., 2015). A maximum e‐value of 1e−40 was set to conservatively estimate the likelihood of similar gene functions between taxonomies. Protein GI accessions from the top hit of diamond alignments were imported into the web‐based version of the david bioinformatics tool (Huang et al., 2009a, 2009b), where corresponding annotations were generated. Given the challenges of functional annotations when dealing with large numbers of loci, we took the conservative approach of annotating only gene homologues known to be associated with virus–host interactions, in particular herpesvirus response pathways, including those in response to HaHV‐1 in the AVG‐resistant H. iris.

RESULTS

Genotyping and overall population structure

Pooled whole genome resequencing of 384 Haliotis rubra specimens from 14 locations yielded a total of 7,745,655 SNPs that were used for population genomic analyses. Estimates of overall genetic structure indicated a lack of structure and panmixia across the sampling distribution. Specifically, global F ST did not differ significantly from zero (F ST = 0.00, p > .05), consistent with reports of panmixia from previous population genomic studies on H. rubra (Miller et al., 2016, 2019). Additionally, no pairwise estimate of F ST between sampling locations was found to differ significantly from zero (Figure 2a). While some sample locations included in this study were chosen based on modelled dependencies on self‐recruitment, our results and those of Miller et al. (2016, 2019) indicate that there is sufficient gene flow across all locations to supress signatures of genetic structure.

FIGURE 2

Heatmap of pairwise estimates of genetic differentiation (F ST) among sample locations based on (a) all 7,745,655 SNPs and (b) the 25,854 SNP loci associated with AVG exposure. *Virus‐unaffected sample locations

Genotype × environment association analysis

Our genome scans found 25,854 candidate SNPs with strong associations with virus exposure (BF > 50), with the PCA based on candidate loci revealing clear patterns of genetic structuring between locations varying in historical AVG exposure (Figure 3). Estimates of F ST based on candidate loci also indicated significant genetic structure among sample locations (F ST = 0.06, p < .05), with pairwise estimates suggesting significant genetic structure between virus‐affected and ‐unaffected sites (Figure 2b). All pairwise estimates of F ST including a single virus‐unaffected location (LJP) were significantly different from zero, but estimates involving this location and other unaffected locations were notably weaker. To account for potential distorting effects associated with LJP, association analyses were repeated including all sites except LJP. Repeat analyses still recovered 21,039 candidate SNPs, with a PCA based on these candidate loci revealing consistent patterns of genetic structuring between locations varying in historical AVG exposure (Figure S1).

FIGURE 3

Plots of eigenvalues from the principal components analysis: (a) plot of axis 1 and 2 eigenvalues, and (b) density plot of axis 1 eigenvalues. Plots are based on candidate SNP genotypes from each of the 28 pooled whole genome resequencing libraries representing virus‐affected (red) and unaffected (black) fishing stocks

Estimates of LD were high at all locations (mean r 2 = .61 ± 0.01 SD) indicating nonrandom association of alleles, while comparisons of r 2 between affected and unaffected stocks did not differ significantly (p > .05). Analyses suggest a genome‐wide pattern, with regression analyses indicating a strong linear relationship between number of candidate loci and scaffold length (R 2 = .83; Figure 4a). However, scaffolds QXJH01000030.1 and QXJH01000212.1 exhibited a higher number of candidate loci relative to scaffold length (deviating from the linear distribution), with candidate SNPs comprising 0.0062% and 0.0108% of the total scaffold nucleotides, respectively. Pairwise linkage among SNPs across the entirety of these scaffolds was high (D′ ≈ 1), with the detection of large haplotype blocks indicating large sections of the genome linked to virus exposure (Figure 4b).

FIGURE 4

(a) Regression analysis indicating a positive linear relationship between number of candidate SNPs (BF > 50) and scaffold length. Outlier scaffolds with a greater frequency of candidate SNPs relative to scaffold length are plotted in red. (b) Linkage disequilibrium heatmaps of scaffolds with the greatest number of candidate SNPs (QXH01000030.1 and QXJH01000212.1) generated with the package ldblockshow (Dong et al., 2021). Heatmaps depict the pairwise linkage disequilibrium measure of D′ (refer to colour key) between each SNP with a BF ≥ 50, while green lines link the relative position of the candidate SNPs to the heatmap. In addition, black triangle sections represent detected haplotype blocks; these are genomic regions of low recombination

Functional annotations

SNP loci showing significant associations with virus exposure were successfully mapped to the annotated H. rubra genome. snpeff analyses predicted 333 candidate loci to have moderate effect on protein function (involving nonsynonomous mutations), while 489 candidates were predicted to have low effect, and 24,722 candidates were recognized as noncoding or intergenic variants. Candidate loci that successfully mapped to H. rubra genome peptide sequences were found to correspond with gene homologues in other animal systems including haliotids, nonhaliotid marine molluscs, crustaceans and humans. These include 13 gene homologues linked to HaHV‐1 immunity in New Zealand pāua (Haliotis iris) and 13 genes associated with herpes virus response pathways in Japanese disk abalone (Haliotis discus hannai), decapod crustaceans (Penaues monodon and Procambarus clarkia) and humans (Table 2). An additional 10 peptides mapped to gene homologues associated with host–virus interactions in various haliotids (H. discus hannai,H. laevigata and H. ruscefens), decapod crustaceans (Pe. monodon and Pr. clarkia) and humans (Table 2). All gene homologues and known functions are provided in Table 2. Notable findings include several genes linked to chitin‐binding peritrophin‐A domain, and the cytochrome P450 (CYP) 3A family, which have recognized associations with immune responses in aquatic molluscs (Badariotti et al., 2007; Zhao et al., 2017) and humans (Fattahi et al., 2018), respectively. Also, the CREB‐binding protein (CBP) was identified, which is associated with herpesvirus responses in humans (Chen et al., 2020; Gwack et al., 2001), as well as immune pathways for C‐type lectins that are important contributors to innate immune responses in invertebrates (Nam et al., 2016; Qin et al., 2019; Zhang et al., 2018).

TABLE 2

List of predicted genetic variant impacts, and genes that candidate loci mapped to. The table also includes gene functions, as well as the species from which these functions have been reported and their respective references

Number of candidate SNPs	Predicted variant impact	Associated gene	Associated gene function	Species	Reference(s)
Genes involved H. iris HaHV−1 immune response
17	Moderate, low, modifier	SLC1A2	Excitatory amino acid transporter 2, response to HaHV−1 exposure	Haliotis iris	Neave et al., 2019
35	Moderate, low, modifier	CYP3A4	Response to HaHV−1 exposure	Haliotis iris	Neave et al., 2019
6	Moderate, modifier	ACE	Response to HaHV−1 exposure	Haliotis iris	Neave et al., 2019
2	Modifier	POU6F2	Response to HaHV−1 exposure	Haliotis iris	Neave et al., 2019
3	Modifier	NLGN4X	Response to HaHV−1 exposure	Haliotis iris	Neave et al., 2019
2	Modifier	CYP3A7	Response to HaHV−1 exposure	Haliotis iris	Neave et al., 2019
1	Modifier	Peritrophin 44 like (LOC105317660)	Chitin‐binding peritrophin‐A domain, response to HaHV−1 exposure	Haliotis iris	Neave et al., 2019
1	Modifier	Uncharacterized LOC105326593 (LOC105326593)	Chitin‐binding peritrophin‐A domain, response to HaHV−1 exposure	Haliotis iris	Neave et al., 2019
17	Modifier	CHIT1	Chitin‐binding peritrophin‐A domain, response to HaHV−1 exposure	Haliotis iris	Neave et al., 2019
15	Modifier	CYP3A5	Response to HaHV−1 exposure	Haliotis iris	Neave et al., 2019
1	Modifier	CYP3A43	Response to HaHV−1 exposure	Haliotis iris	Neave et al., 2019
14	Modifier	Ganglioside GM2 activator like (LOC105346019)	Chitin‐binding peritrophin‐A domain, response to HaHV−1 exposure	Haliotis iris	Neave et al., 2019
1	Modifier	Ganglioside GM2 activator like (LOC105348613)	Chitin‐binding peritrophin‐A domain, response to HaHV−1 exposure	Haliotis iris	Neave et al., 2019
Gene homologues associated with herpesvirus response pathways
2	Modifier	CASP8	Herpesvirus response pathway	Haliotis discus hannai	Nam et al., 2016
1	Modifier	TAF10	Herpesvirus response pathway	Homo sapiens	Wagner & DeLuca, 2013
1	Modifier	ARNTL	Herpesvirus response pathway	Homo sapiens	Edgar et al., 2016
7	Modifier	EEF1D	Herpesvirus response pathway	Homo sapiens	Boulben et al., 2003
2	Modifier	SRSF7	Herpesvirus response pathway	Homo sapiens	Tang et al., 2019
6	Modifier	Ube2r2	Herpesvirus response pathway	Homo sapiens	Beard et al., 2015
2	Modifier	CREBBP	Herpesvirus response pathway, C‐type lectin containing or involved in c‐type lectin pathway (invertebrate immune response)	Haliotis discus hannai,Procambarus clarkii,Penaeus monodon	Nam et al., 2016, Zhang et al., 2018, Qin et al., 2019, Wang & Wang, 2013
2	Modifier	SRPK1	Herpesvirus response pathway	Homo sapiens	Souki & Sandri‐Golden, 2009
2	Modifier	TAF6	Herpesvirus response pathway	Homo sapiens	Wagner & DeLuca, 2013
1	Modifier	TAB1	Herpesvirus response pathway	Homo sapiens	Jahanban‐Esfahlan et al., 2019
1	Modifier	Csnk2b	Herpesvirus response pathway	Homo sapiens	Carter, 2011
4	Modifier	Hnrnpk	Herpesvirus response pathway	Homo sapiens	Schmidt et al., 2010
1	Modifier	PPP1CA	Herpesvirus response pathway	Homo sapiens	Silva et al., 2015
Gene homologues associated with host–virus interactions
4	Modifier	HIST1H2AA	Abalone immune response	Haliotis discus hannai	Nam et al., 2016
1	Low	COTL1	Tropomyosin, abalone immune response	Haliotis discus hannai	Nam et al., 2016
1	Low, modifier	PDIA3	Protein disulphide isomerase activity, abalone immune response	Haliotis discus hannai	Nam et al., 2016
1	Modifier	HSP90AB1	Abalone immune response	Haliotis discus hannai	Nam et al., 2016
4	Modifier	Histone H2A (LOC105320412)	Abalone immune response	Haliotis discus hannai	Nam et al., 2016
8	Low, modifier	AVIL	Gelsolin domain, abalone immune response	Haliotis discus hannai	Nam et al., 2016
1	Modifier	Hsp70 member 12A (HSPA12A)	Abalone immune response	Haliotis laevigata,Haliotis ruscefens,Haliotis discus hannai	Shiel et al., 2015, Brokordt et al. 2015, Nam et al., 2016
3	Modifier	PSMA8	C‐type lectin containing or involved in c‐type lectin pathway (invertebrate immune response)	Haliotis discus hannai,Procambarus clarkii,Penaeus monodon	Nam et al., 2016, Zhang et al., 2018, Qin et al., 2019, Wang & Wang, 2013
1	Modifier	IRAK4	Toll‐like receptor activity, innate immunity	Crassostrea gigas,Haliotis diversicolor	Tang et al., 2017, Ge et al., 2011
2	Modifier	TIA1	Viral translation inhibition	Homo sapiens	McCormick & Khaperskyy, 2017

DISCUSSION

We provide evidence of rapid genetic changes in as little as a single generation in wild Haliotis rubra populations decimated by HaHV‐1. Specifically, our genome scans identified SNP loci associated with virus exposure, many of which mapped to genes known to contribute to HaHV‐1 immunity, herpes virus response pathways, and host–virus interactions in haliotids and other animal systems. These findings require experimental validation but are consistent with rapid evolutionary changes in H. rubra fishing stocks impacted by disease, with stock recovery potentially influenced by evolved resistance. This study highlights the pace at which selection can act to counter disease in wildlife communities by leading to an increased frequency of potentially adaptive genotypes. The implications of these findings are discussed in the context of future infectious diseases management in abalone fisheries and wildlife conservation more generally.

Evidence of rapid evolution in H. rubra

Our analyses identified 25,854 SNP loci associated with AVG exposure in H. rubra. Multiple lines of evidence point to selection being responsible for driving specific genetic variants to higher frequencies in virus‐affected populations, rather than random demographic factors such as bottleneck effects. Random demographic factors are capable of causing shifts in allele frequencies (particularly rare alleles) in populations that have suffered major declines, depending on the number of surviving individuals and the influence of random genetic drift in subsequent generations (Willi et al., 2022). However, the number of surviving abalone at virus‐impacted locations still consisted of many 100s to 1000s of individuals (Western Abalone Divers Association, 2020) expected to maintain genetic diversity. Also, we have shown that the observed patterns of differentiation between virus‐affected and unaffected fishing stocks are not being driven by rare alleles (Figure S2). Additionally, opportunities for allele frequencies to be influenced by random genetic drift are limited given the animals included in this study are either survivors or F1/F2 progeny of those survivors. Furthermore, we observed consistent genetic changes across a spatially replicated and stratified sampling distribution, which are unlikely to arise randomly in the absence of selection. Finally, we provide evidence of functional associations of candidate SNPs to genes known to contribute to disease adaptation in sister haliotid taxa, including a number of genes that contribute to HaHV‐1 immunity.

Despite the large number of candidate loci found to be associated with AVG exposure, it is likely that only some are directly influenced by selection and potentially contributing to AVG immunity. Evidence of LD between loci, and the noncoding and intergenic nature of most SNPs indicates many indirect associations due to physical linkage with adaptive loci that are embedded in large haploblocks and responding to selection (Uffelmann et al., 2021). Further, it is possible that some candidate loci are in statistical linkage, sharing similar allele frequency distributions across haploblocks and chromosomes (i.e. genetic indistinguishability; Skelly et al., 2016), and inflating the overall number of candidate loci. Despite the potential for linkage, a strong linear relationship between genome scaffold length and density of candidate SNPs was observed, suggesting a genome‐wide response to virus exposure. Analyses indicate a higher incidence of candidate SNPs on two scaffolds, which could represent important gene‐rich adaptive regions (Hohenlohe et al., 2012). However, further work is needed to test their functional significance as the genes annotated within these regions are not known to contribute to virus resistance. Nevertheless, our analyses indicate that some candidate loci may be directly involved in disease adaptation in H. rubra.

In particular, functional annotations of candidate loci point to associations with genes and protein domains that contribute to HaHV‐1 immunity in the New Zealand pāua (Haliotis iris). Neave et al. (2019) first characterized genes associated with HaHV‐1 immunity in H. iris through transcriptome analyses of animals subject to HaHV‐1 immersion challenge tests. Their study was the first to characterize the molecular basis of HaHV‐1 immunity in a haliotid species, and our study complements these findings by identifying a common set of genes involved in haliotid host response to HaHV‐1 exposure. Functional annotations of candidate loci also point to associations with genes and protein domains contributing to herpes simplex virus responses and immune responses in various haliotids and other animal systems. Importantly, our analyses indicate that a large number of candidate loci involve nonsynonymous mutations that are expected to have an effect on protein structure and function. It is also possible that some candidate SNPs recognized as being noncoding and intergenic variants might provide functionally important regulatory roles in genomic processes (Gusev et al., 2014; Li et al., 2012; Wei et al., 2020). Overall, these findings strongly support the notion that divergent adaptation, involving a polygenic response and possible selection for disease resistance phenotypes, has occurred in H. rubra fishing stocks impacted by AVG.

The results of this study point to rapid genetic changes in H. rubra fishing stocks impacted by disease. However, experimental validation will be needed to link any genetic changes to disease resistance. Challenge tests involving the exposure of animals with putatively resistant genotypes to HaHV‐1 will help to determine if, and how much, HaHV‐1 resistance is determined by the candidate genotypes (Corbeil et al., 2016; Crane et al., 2013). Although previous challenge tests performed by Crane et al. (2013) showed no, or very low, resistance to HaHV‐1 in H. rubra, the animal source locations differ from those included in our study. This points to potential spatial variation in adaptive responses across the fishery, and possible genotype–environment interactions affecting the expression of resistant phenotypes in different environmental settings (Hoffmann et al., 2020). Nevertheless, a large number of loci appear to be responding to the virus in our study, suggesting that changes in resistance will represent a polygenic response that can be followed by controlled breeding studies (Guarna et al., 2017; Gutierrez et al., 2018). The response to selection in these studies will depend on levels of heritable variation as well as the intensity of selection, which will determine the impact of the phenotypic change. Breeding studies and experimental trials will also be essential to assess trait heritability and genotype–environment interactions, particularly if industry intends to control for resistance traits in a culture environment for breeding purposes (discussed below).

Implications for fisheries management

Re‐emergence of AVG remains a significant threat to the economic viability of H. rubra fisheries in southeastern Australia (Corbeil, 2020; Lafferty et al., 2015). Therefore, characterizing the spatial distribution and prevalence of disease‐resistant genotypes will help managers identify stocks expected to be either resilient or vulnerable to AVG re‐emergence. Previous population genomic research has indicated a lack of biological stock structure in these fisheries (Miller et al., 2016), suggesting that gene flow could contribute to the spread of adaptive genotypes and resilience of naïve fishing stocks. Gene flow from unaffected parts of the fishery could also eventually reduce the frequency of adaptive genotypes over time in the absence of ongoing selection, but such effects are yet to become apparent. Whether selection for resistance will be a recurring process is still unclear, although the recent outbreak recorded in 2021 at several fishing locations previously impacted by AVG in the early 2000s indicates that this is a possibility. Notably, the recent outbreak has resulted in minimal animal mortality and disease spread, supporting the notion of evolved resistance following previous exposure to disease. Contrasting the genetic profiles of survivors and those that have succumbed to the virus at these locations as a result of the most recent outbreak will help to reinforce the findings of the current study. Similarly, the recent outbreak provides a unique opportunity to genomically select putatively resistant and vulnerable animals for challenge tests aimed at providing functional validation of AVG resistance.

Evidence of panmixia in H. rubra (Miller et al., 2016) suggests that standing genetic variation is likely to persist within disease‐naïve populations allowing for in situ adaptation to HaHV‐1. However, strategic stock augmentation activities, involving the translocations of animals with AVG‐resistant genotypes, could potentially assist the spread of genotypes to reduce risks of vulnerability across wild fisheries. Also, there may be future opportunities to biosecure farm fisheries through the establishment of AVG‐resistant breeding programmes, similar to disease‐related breeding programmes in other farmed mollusc, crustacean and finfish fisheries around the world (Kjøglum et al., 2008; Moss et al., 2012; Potts et al., 2021; Ragone Calvo et al., 2003). Overall, these results add to those of Miller et al. (2019) demonstrating patterns of genetic adaptation across environmental gradients and the adaptability of H. rubra populations to new environmental conditions. This is pertinent in southeastern Australia where rapid changes in the physical marine climate are threatening commercial fisheries through shifts in species distributions (Johnson et al., 2011; Ling, 2008), changes in habitat and trophic interactions (Holland et al., 2021), and risks of infectious diseases (Oliver et al., 2017).

CONCLUSIONS

While it has been proposed that selection in some species may fail to keep pace with rates of pathogen evolution and the emergence and turnover of novel diseases (Hawley et al., 2013; Ujvari et al., 2014), our study demonstrates that rapid evolutionary responses within a single generation may be possible in large populations under extreme selective pressure. This study highlights the value of genome scans for identifying signatures of potential adaptation among natural populations differing in virus exposure and characterizing putative genes that contribute to disease‐resistant phenotypes. Future studies of this nature will be critical for understanding the potential for rapid evolution in other species threatened by infectious diseases, informing the management of new outbreaks, and securing populations with little or no resistance.

AUTHOR CONTRIBUTIONS

This project was conceived by A.D.M., while O.J.H., M.T., A.D.M. and L.C. were responsible for generating and analysing the data, with assistance from C.A. and A.A.H. Writing of the manuscript was led by O.J.H. and A.D.M. with assistance from all authors.

CONFLICT OF INTEREST

The authors declare no conflict of interest.

OPEN RESEARCH BADGES

This article has earned an Open Data Badge for making publicly available the digitally‐shareable data necessary to reproduce the reported results. The data is available at https://doi.org/10.5061/dryad.b5mkkwhfb.

Fig S1‐2

Click here for additional data file.