Expanding the phylogeography and connectivity of Perkinsus species across North and Central America

Parasites in the genus Perkinsus infect marine molluscs globally, with novel detections expanding and reshaping our knowledge of their biogeographic patterns and the factors influencing those patterns. Here, we aimed to characterize the phylogeography and genetic connectivity of Perkinsus spp. in bivalves across North and Central America, which included infection hot spots (e.g., the Gulf of Mexico, Chesapeake Bay) and areas where these parasites had not been previously reported (e.g., along the California coast).

However, parasite biogeography remains elusive, which is often attributed to challenges in detecting them due to their small size, resulting in data limitations and their lack of inclusion in most biodiversity studies (Carlson et al., 2020).Our limited knowledge of marine parasite biogeography hinders our ability to detect and assess explanations for reported range expansions, which could result from changing environmental conditions, anthropogenic dispersal via shipping or aquaculture (Bax et al., 2003;Goedknegt et al., 2016;Pagenkopp Lohan et al., 2018).Additionally, without prior knowledge of parasite presence in an area, it is difficult to distinguish between the emergence and re-occurrence of disease within host populations (Harmon et al., 2019;Hewson et al., 2014).
Previously, we examined the diversity and host specificity of four Perkinsus spp.found for the first time in Panama (Pagenkopp Lohan et al., 2016), then assessed their global connectivity (Pagenkopp Lohan et al., 2018).We found identical parasite haplotypes across broad geographic scales (e.g., across ocean basins) indicating a high level of global genetic connectivity for some species (Pagenkopp Lohan et al., 2018).Our prior results led to the hypothesis that Perkinsus spp.geographic ranges result from a complex interaction between abiotic (i.e., salinity, temperature), biotic (i.e., host specificity, availability and immunity) and anthropogenic (e.g., maritime shipping, aquaculture) factors influencing the dispersal and survivability of these parasites (Pagenkopp Lohan et al., 2018).The data reported herein explored the phylogeography of these parasites across North America and Central America, including locations where mortality events in molluscs have not been previously reported, to expand our understanding of biogeographic patterns across the region.
We sampled bivalve populations along both coasts of the USA and the Gulf Coast of Mexico, then compared these data to our prior results from bivalves collected in Panama (Figure 1).Our goals for this study were to: (1) assess the biogeographic ranges of P. chesapeaki, P. beihaiensis and P. marinus across North and Central America, and (2) examine the genetic connectivity of haplotypes across regional, continental and global scales.Our results indicate that Perkinsus species inhabit broader geographic areas than previously reported, that genetic connectivity within species is high, and that these species occur in relatively high prevalence in areas where mass mortality events have not been reported.

| Sample collection
We collected bivalves from nine genera from different habitats within six regions in North and Central America between December 2012 and November 2015 (Table S1).On average, 64 individuals (range = 9-175) per species per region were collected (Table 1).Collections from Panama, Mexico, Virginia and Florida (Atlantic), USA were done systematically as described in Pagenkopp Lohan et al. ( 2016) and processed on-site.Collections from the US Gulf of Mexico (Florida, Mississippi and Texas) and California, USA were opportunistic and were shipped live on ice overnight to the Smithsonian Environmental Research Center (Edgewater, MD, USA).Bivalves were collected by removing individuals or clusters of individuals from the hard substrate, using oyster tongs (when submerged), or via dredge (Texas only).Upon collection, oysters were placed in coolers on ice and kept at ∼4°C in the laboratory for no more than 72 h.The oysters were tentatively identified based on morphology.Shells were thoroughly cleaned, dried, labelled with unique identification numbers and retained as vouchers.Pieces of gill, mantle and digestive gland from each individual were removed and preserved in 95% ethanol for subsequent molecular analysis.

| DNA extraction, amplification, sequencing
Sequences for the bivalve hosts and parasites were generated as de- The mitochondrial cytochrome oxidase I (COI) gene of the bivalves was amplified using primers jgLCOI490Fand jgHCO2198R (Geller et al., 2013).PCR reagents consisted of 1× GeneAmp 10× PCR Gold Buffer (150 mM Tris-HCl, pH 8.0; 500 mM KCl; Applied Biosystems), 1.5 mM MgCl 2 , 0.2 mM each nucleotide, 0.4 μM each primer, 0.2 mg/mL bovine serum albumin (BSA; New England Biolabs), and 0.03 units/μL of AmpliTaq Gold with water to a final volume of 20 μL.Thermocycling was carried out using a S1000 Thermal Cycler (Bio-Rad) with an initial denaturation of 94°C for 10 min, 35 cycles of 94°C for 30 s, 48-52°C for 30 s, 72°C for 1 min and a final elongation of 72°C for 5 min.The PCR product (5 μL) was electrophoresed on an agarose gel (2% w/v) stained with GelRed (Phenix Research) and visualized under UV light.Any sample that did not amplify for this assay (n = 3) was not used for parasite screening.
To screen for Perkinsus species, we used the genus-specific primers PerkITS85FNEW forward and PerkITS750R reverse (Casas et al., 2002;Moss et al., 2007) that target the first internal transcribed spacer region (ITS1) of the ribosomal gene complex (rDNA).Additionally, as this was the fragment used in our prior studies (Pagenkopp Lohan et al., 2016, 2018), this allowed us to ensure that we could compare these new amplicons to those we had generated previously.PCR reagents consisted of 1× AmpliTaq Gold PCR Buffer (150 mM Tris-HCl, pH 8.0, 500 mM KCl; Applied Biosystems), 1.5 mM MgCl 2 , 0.2 mM each nucleotide, 0.5 μM each primer, 0.2 mg/mL BSA (New England Biolabs), 0.025 units/ μL of AmpliTaq Gold, with water to a final volume of 20 μL.Thermocycling was carried out using an S1000 Thermal Cycler (Bio-Rad) with an initial denaturation of 94°C for 10 min 40 cycles of 94°C for 30 s, 55°C for 1 min, 72°C for 1 min and a final elongation of 72°C for 5 min.Positive control samples, consisting of extracted genomic DNA from infected bivalves that had successfully amplified in the past, were obtained from our own collection for Perkinsus species.Amplicons were visualized as described above.
All ITS1 amplicons were purified with a 1:10 dilution of Exo-SAP IT (Affymetrix, Santa Clara, CA).Thermocycling was carried out for 30 min at 37°C and 15 min at 80°C.Cleaned amplicons were then bi-directionally cycle-sequenced using the Big Dye Terminator Cycle Sequencing Kit v3.1 (Applied Biosystems, Inc.) on an ABI 3130 Sequencer (Applied Biosystems, Inc.).Sequences were edited using Sequencher 5.4.6 (Gene Code Corporation, Ann Arbor, MI).Only unique sequences, defined as sequences differing by either a single base pair or gap were used in the phylogenetic analyses.

| Bivalve phylogenetic analyses
For the bivalve hosts, COI sequences were checked for open reading frames and only unique sequences were used to create phylogenetic trees.Sequences were aligned in Geneious Prime® 2020.1.1 (Biomatters Ltd., San Diego, CA) with the MAFFT v7.450 (Katoh et al., 2002;Katoh & Standley, 2013) plug-in, using the default parameters.We used the Gblocks server (Castresana, 2000) with both more stringent (not allowing many contiguous nonconserved positions) and less stringent (allowing smaller final blocks, gap positions within the final blocks, and less strict flanking positions) parameters to remove the less conservative portions of the alignment.As the alignments and resulting trees were highly similar or identical, we report the alignment that was used to generate the reported tree.The COI alignment for the Ostreidae dataset was 623 bp with the relaxed setting.The COI alignment for the Geukensia spp.sequences was 524 bp with the stringent TA B L E 1 Total number of oyster species collected and screened for Perkinsus species from across North and Central America sampled.corrected values.This process allowed for Bayesian analyses to use the appropriate number of available substitution models.The Bayesian analyses for each dataset using the suggested substitution model were performed in Geneious Prime® 2020.1.1 using MrBayes 3.2.6 (Ronquist & Huelsenbeck, 2003).

| Perkinsus phylogenetic analyses
A subset of sequences for all Perkinsus spp.were obtained from GenBank and used to generate a phylogenetic tree to determine parasite species.Sequences were aligned in Geneious Prime® 2020.1.1 (Biomatters Ltd., San Diego, CA) with the MAFFT v7.450 (Katoh et al., 2002;Katoh & Standley, 2013) plug-in, using the default parameters.Alignments were trimmed to the longest of our sequences, and then terminal gaps were filled with N bases for those sequences that were shorter than the alignment.We used the Gblocks server v0.91b (Castresana, 2000) as described above to remove the less conservative portions of the alignment.The resulting alignment was 624 bp.IQTree Model Finder (Trifinopoulos et al., 2016) was used as described above.We used MrBayes 3.2.6 plug-in (Ronquist & Huelsenbeck, 2003) with the following parameters: GTR Evolutionary Model, Gamma Variable rate and 4 variable sites and PhyML 3.3.20 (Guindon et al., 2010) using TPM3 model Custom code 012012 option with 1000 bootstrap replicates.P.
quagwadi was chosen as the outgroup as prior studies have demonstrated that this species is genetically distinct from others in the genus (Moss et al., 2008;Villalba et al., 2004).

| Perkinsus haplotype networks
To generate global haplotype networks for the three parasite species, we included only those sequences available in GenBank™ that contained information on the geographic origin of the host.
Because many GenBank™ records did not include specific regions of a country where the host bivalve was collected from, we were unable to build networks with geographic specificity below the country level.For each Perkinsus spp., we aligned sequences as above.Alignments were trimmed to the longest of our sequences, and then terminal gaps were filled with N bases for those sequences that were shorter than the length of the alignment.Doing this allowed us to use the greatest amount of sequence information for the regional and global haplotype networks to maximize the power of the analyses for examining genetic connectivity at regional and global scales.The three alignments used to generate global networks included 126 sequences for P. beihaiensis and was 665 bp, 95 sequences for P. chesapeaki and was 578 bp and 148 sequences for P. marinus and 863 bp.Haplotype networks were generated using these alignments with the TCS algorithm (Clement et al., 2002) in PopArt (Leigh & Bryant, 2015) Additional alignments were created as described above using only the sequences we generated, which were then used to calculate Tajima's D (Tajima, 1993)
All the bivalves collected had morphological identifications verified with phylogenetic methods.We confirmed identifications for two species of Isognomon, including an unnamed species described in Pagenkopp Lohan et al. ( 2015) and I. bicolour (Figure S1).Additionally, we confirmed identifications for Mercenaria mercenaria (Figure S2), Geukensia demissa (Figure S3), Crassostrea virginica, C. gigas, Ostrea stentina and O. lurida (Figure S4).Thus, we confirmed the identifications of 16 species from three different countries across North and Central America and tested them via molecular methods for the presence of Perkinsus spp.

Perkinsus beihaiensis
Furthermore, 15 sequences across all regions except Panama could only be identified as Perkinsus spp., which we suspect represent multiple detections of two or more Perkinsus species or haplotypes due to the apparent double peaks noted in the sequence chromatograms.

| Global diversity and distribution of Perkinsus haplotypes
The P. beihaiensis global network contained 27 haplotypes from 7 countries (Figure 3).Sequences used to generate the haplotype network had a range in genetic similarity of 82%-100%, with a range of single nucleotide polymorphisms (SNPs) from 1 to 11.
For P. beihaiensis, the greatest number of haplotypes (n = 23) were found in China, with 21 of those being unique to China.
The P. marinus global network contained 34 haplotypes from five countries (Figure 4).Sequences used to generate the haplotype network had a range in genetic similarity of 94.54% -100%, with a range of SNPs from 1 to 15.For the P. marinus network, the greatest number of haplotypes (n = 11) were found in Brazil.The most common haplotype was found in Panama, Brazil, Mexico, and multiple locations in the USA.The lone haplotype found in India was most like a haplotype found in Mexico but had a 15 bp difference.Most of the P. marinus F I G U R E 3 Global network analysis for Perkinsus beihaiensis generated with the PopArt software (Leigh & Bryant, 2015) using sequences from North America (this study) and additional sequences from GenBank that contained the sample location in their records.Colours indicate locations.Each crossbar indicates a single nucleotide polymorphism or a gap.The proportion of colour in pie charts with >1 sequence indicates the proportion of unique sequences obtained from that country.CA_USA, California, USA; AtlFL_USA, Atlantic Coast of Florida, USA.

TA B L E 2
Statistics for sequences generated from this study to assess the genetic diversity of Perkinsus spp.across North America calculated in MEGA11 (Tamura et al., 2021).
paramount for understanding these host-parasite interactions and preventing further anthropogenic dispersal of these parasites.

| Perkinsus prevalence and host specificity
In this study, we report the expansion of the potential host ranges of Perkinsus species.For example, we report novel detections of P.
beihaiensis in O. lurida, O. stentina, and C. virginica.P. marinus was primarily associated with C. virginica, particularly in areas where C.
virginica is highly abundant (e.g., Gulf of Mexico, Chesapeake Bay).This is in line with previous research where P. marinus was hypothesized to be more likely to infect oysters compared to other sympatric bivalves (Reece et al., 2008).Interestingly, we did not detect P. olseni outside of Panama, even though P. olseni is widespread across other continents with the widest range of host usage of all the Perkinsus species, having previously been detected in 27 host species (Pagenkopp Lohan et al., 2016;Villalba et al., 2004).
As all bivalve hosts were screened using PCR methods and infections were not confirmed with culturing using Ray's Fluid Thioglycollate Medium (RFTM) or histopathology, we follow the guidelines described in Burreson (2008)

| Global genetic diversity and connectivity
In these analyses, we used more sequence information in longer align- and globally dispersed through maritime traffic, a major mechanism for transporting goods between North and Central America (Ruiz et al., 2015) as well as globally.While natural dispersal (e.g., via ocean currents, rafting with hosts) is possible and has been hypothesized as a potential long-distance dispersal method for another bivalve parasite (Hill-Spanik et al., 2015), there are multiple lines of evidence suggesting maritime traffic as the primary anthropogenic dispersal mechanism.First, previous research demonstrated that hypnospores (i.e., pre-zoosporangia) can survive harsh conditions (e.g., low temperature, salinity or pH) for extended periods (Klein et al., 2010;McCollin et al., 2007;Pagenkopp Lohan et al., 2018;Villalba et al., 2004)

| Drivers of biogeographic patterns
While biogeographic patterns of parasites are likely shaped by the intersection of biotic and abiotic drivers, more research is needed to examine the main drivers for influencing biogeographic patterns for Perkinsus spp.Previously, three of the seven known Perkinsus species (P.marinus, P. beihaiensis, P. olseni) were detected in tropical waters, with P. beihaiensis being the only one exclusively found in tropical waters (Moss et al., 2008;Pagenkopp Lohan et al., 2016;Villalba et al., 2004).Thus, it was initially suspected that seawater temperature and salinity were the leading drivers affecting the biogeography of these three Perkinsus spp.
across regional and global scales.In this study, P. beihaiensis was recovered from native oysters in San Francisco, California, USA, so it does not solely reside in tropical waters.There are two potential explanations for this finding.The first is that the limited sampling of Perkinsus species has hampered our understanding of the biogeographic patterns for this parasite.As we find new Perkinsus associations along regional and global scales, it is certainly possible that our understanding of the factors shaping biogeographic patterns will simply continue to change with new distribution data.Additionally, it is possible that the impacts of global climate change on the abiotic conditions of coastal waters are allowing these parasites to expand their biogeographic boundaries.Most notably, this was reported for P. marinus, which expanded its biogeographic range north along the eastern seaboard of the USA, as surface water temperatures increased (Ford & Chintala, 2006).
Changing environmental conditions could influence parasite distributions in multiple ways, including increasing the susceptibility of the host to infection and/or increasing the proliferation of the parasite (Burrell Jr et al., 1984;Powell et al., 1996;Soniat, 1985).
Hosts in geographic areas experiencing warmer winters and drier summers could suffer increased incidence and severity of disease, in addition to the parasites continuing to expand their biogeographic ranges, potentially encountering novel and more susceptible host species.Perkinsus species do not cause overt signs of disease (e.g., dehydrated flesh pulled in from the edge of the shell) in all locations where they have been detected (i.e., Panama, San Francisco), suggesting that continued surveillance of natural communities for the presence of these pathogens is vital to monitor their presence and potential future impacts with changing environmental conditions in the world's oceans.
scribed in Pagenkopp Lohan et al. (2016).Briefly, following overnight K E Y W O R D S geographic range, haplotype, oyster, parasites, Perkinsus, phylogeography digestion with proteinase K, genomic DNA was extracted from the gill, mantle and digestive gland for each bivalve, which were pooled into a single extraction, using a Qiagen Biosprint Kit (Qiagen, Germantown, MD) following the manufacturer's protocols for animal tissues.All extractions completed within the same day included a blank extraction, which served as a negative extraction control for PCR.Aliquots of the extracted DNAs (50 μL) were made to avoid contamination of stock DNA elutions and stored at 4°C while stock DNA elutions were stored at −20°C.

F
Map of sampling regions and locations in North and Central America.Colours indicate which Perkinsus parasite species were identified in each location.The proportion of colour in each pie chart indicates the proportion of each species that was identified by PCR in each location.The Perkinsus spp.infections appear to be multiple infections, based on multiple peaks in chromatograms (see text for more details).
settings.The COI alignment for the Mercenaria spp.sequences were 511 bp with the ends trimmed and relaxed settings.Finally, the COI alignment for the Isognomon spp.sequences was 425 bp for the relaxed settings.The IQTree Model Finder(Trifinopoulos et al., 2016) was used to determine the best substitution model for the alignments based on Akaike Information Criterion (AIC) in MEGA11(Tamura et al., 2021).For this statistic, a negative value indicates a recent population expansion or selective sweep, a positive value indicates either a balancing selection or a population contraction and a value of zero indicates no selection.
The number of sequences, nucleotide diversity (π), number of segregating sites, and the result of Tajima's D Neutrality test are shown.These statistics could not be calculated for Perkinsus chesapeaki due to too few novel sequences for this parasite (n = 2).
and reference PCR positives or host associations, not infections.In most cases, we recovered more than one PCR-positive individual per host per location, providing evidence that the parasites are present in those locations and likely infecting that host population, but further confirmation is required to assess actual prevalence and intensity values.Given the apparent double peaks in the sequence chromatograms, which we interpreted as dual infections of either multiple Perkinsus species or haplotypes, we likely underestimated the distribution of all Perkinsus species and haplotypes found.Future studies should utilize either genomic capture techniques or high throughput sequencing efforts to better assess the species and genetic diversity of Perkinsus species across hosts and geographic areas.
compared to the analyses in Pagenkopp Lohan et al.(2018), which allowed us to compress some of the haplotypes and see genetic connectivity that was not previously obvious.Our results expand current knowledge regarding the genetic connectivity and biogeography of these Perkinsus spp., particularly across North and Central America.For all the parasites examined, identical haplotypes occur across vast geographic spreads, including across continents and ocean basins.Due to the rapidly evolving nature of the loci used in this study, these findings indicate either rapid natural or anthropogenically assisted dispersal over rather short periods of time.Additionally, most of these networks also include at least some haplotypes that differ from others by many base pairs, indicating that we have likely undersampled the full genetic diversity of haplotypes across these species.Our results provide further evidence that Perkinsus spp.are likely being regionally F I G U R E 5 Global network analysis for Perkinsus chesapeaki generated with the PopArt software using sequences from North America (this study) and additional sequences from GenBank that contained the sample location in their records.Colours indicate locations.Each cross-bar indicates a single nucleotide polymorphism or a gap.Numbers above or inside pie charts indicate the number of consensus sequences that are represented.The proportion of colour in pie charts with >1 sequence indicates the proportion of unique sequences obtained from that country.Pie charts without numbers represent a single consensus sequence.VAMD_USA, Virginia & Maryland, USA.
Phylogram for Perkinsus spp.constructed using ITS1 sequences generated in this study (in bold) and GenBank sequences.The GTR + G substitution model for the Bayesian analysis was used to compute this topology, and posterior probabilities (>90) are included at the nodes.See text for additional details.
While the majority of P. behaiensis haplotypes were from bivalves located in tropical waters, one haplotype spanned temperate and tropical locations and was recovered from Atlantic Florida, USA, Japan, India, Mexico, China and California, USA.We found little genetic diversity in P. beihaiensis sequences from across North America andF I G U R E 2