Exploration of the 2016 Yellowstone River ﬁ sh kill and proliferative kidney disease in wild ﬁ sh populations

. Proliferative kidney disease (PKD) is an emerging disease that recently resulted in a large mortality event of salmonids in the Yellowstone River (Montana, USA). Total PKD ﬁ sh mortalities in the Yellowstone River were estimated in the tens of thousands, which resulted in a multi-week river closure and an estimated economic loss of US$500,000. This event shocked scientists, managers, and the public, as this was the ﬁ rst occurrence of the disease in the Yellowstone River, the only reported occurrence of the disease in Montana in the past 25 yr, and arguably the largest wild PKD ﬁ sh kill in the world. To understand why the Yellowstone River ﬁ sh kill occurred, we used molecular and historical data to evaluate evidence for several hypotheses: Was the causative parasite Tetracapsuloides bryosalmonae a novel invader, was the ﬁ sh kill associated with a unique parasite strain, and/or was the outbreak caused by unprecedented environmental conditions? We found that T. bryosalmonae is widely distributed in Montana and have documented occurrence of this parasite in archived ﬁ sh collected in the Yellowstone River prior to the ﬁ sh kill. T. bryosalmonae had minimal phylogeographic population structure, as the DNA of parasites sampled from the Yellowstone River and distant water bodies were very similar. These results suggest that T. bryosalmonae could be endemic in Montana. Due to data limitations, we could not reject the hypothesis that the ﬁ sh kill was caused by a novel and more virulent genetic strain of the parasite. Finally, we found that single-year environmental conditions are insuf ﬁ cient to explain the cause of the 2016 Yellowstone River PKD outbreak. Other regional rivers where we documented T. bryosalmonae had similar or even more extreme conditions than the Yellowstone River and similar or more extreme conditions have occurred in the Yellowstone River in the recent past, yet mass PKD mortalities have not been documented in either instance. We conclude by placing these results and unresolved hypotheses into the broader context of international research on T. bryosalmonae and PKD, which strongly suggests that a better understanding of bryozoans, the primary host of T. bryosalmonae , is required for better ecosystem understanding.


INTRODUCTION
The Upper Yellowstone River in Montana, USA, supports a wide variety of uses, including tourism, agriculture, industry, and recreation. Its vital river fishery stimulates over US$150 million annually to the local economy (Sage 2016). This river is also a major headwater for the upper Missouri River and a gateway into Yellowstone National Park, which is annually visited by millions of people from across the world. The health of the aquatic ecosystems in the area is therefore of national and global interest, and the emergence of fish mass mortality events in the upper Yellowstone River has attracted wide attention and concern (Lamborn and Smith 2019).
In early August of 2016, Montana Fish Wildlife and Parks (FWP) began receiving angler reports of dead and moribund mountain whitefish (Prosopium williamsoni) on the Upper Yellowstone River (Fig. 1). FWP initiated surveys near Emigrant and Livingston, Montana to determine the cause and the extent of mortalities. During the surveys, over 6000 mortalities were recorded (Opitz and Rhoten 2017), resulting in estimates of total fish mortalities in the tens of thousands for the event. Mountain whitefish accounted for >99% of observed mortalities, but population estimates were not available prior to the outbreak so the proportion of the mountain whitefish population that died is unknown. Mortalities were also observed for brown trout (Salmo trutta), rainbow trout (Oncorhynchus mykiss), Yellowstone cutthroat trout (Oncorhynchus clarkii bouvieri), longnose sucker (Catostomus catostomus), and longnose dace (Rhinichthys cataractae). A subsample of both dead and apparently healthy fish from these six species was sent to the U.S. Fish and Wildlife Services Bozeman Fish Health Center (BFHC; Bozeman, Montana, USA) for diagnosis. Microscopic and molecular methods identified the cause of mortality in the salmonid fish as Proliferative Kidney Disease (PKD). The cause of mortality in non-salmonid species (i.e., longnose suckers and dace) was inconclusive and assumed to be heat-related stress. Healthy salmonids submitted for diagnosis lacked the hallmark kidney inflammation of PKD. However, the BFHC observed exceptionally high densities of Tetracapsuloides bryosalmonae, the causative parasite of PKD, in the kidney interstitial spaces of mountain whitefish using impression smears. The presence of T. bryosalmonae was confirmed in the kidneys of apparently healthy fish and those that were dead or moribund upon collection using conventional polymerase chain reaction (PCR; Kent et al. 1998) T. bryosalmonae DNA was also detected by PCR in several other tissues, including gill, spleen, and liver.
Tetracapsuloides bryosalmonae has a complex life cycle that involves freshwater bryozoans (Phylum Bryozoa: Class Phylactolaemata) and salmonid fish hosts. Freshwater bryozoans are sessile suspension-feeding invertebrates that develop as spreading colonies attached to hard surfaces (Wood and Okamura 2005) and act as primary hosts (supporting sexual reproduction). T. bryosalmonae infections have been detected in bryozoan species in the genera Fredericella and Plumatella (Canning and Okamura 2003). In Europe, the most common host of T. bryosalmonae is Fredericella sultana and most studies of bryozoan-parasite interactions focus on that system . The parasite suppresses growth and reproduction of bryozoan hosts during brief periods of spore development but otherwise causes no apparent effect on growth or mortality (Tops et al. 2009, Hartikainen and. This low impact on bryozoan hosts and the ability of the parasite to spread along with host propagation (via colony fragmentation and the production of dormant propagules) is likely to contribute to the persistence of infections over space and time in bryozoan populations (Fontes et al. 2017b). Fish-to-bryozoan and bryozoan-to-fish transmission are achieved by spores released into the water, and no fish-to-fish transmission takes place (D'Silva et al. 1984). Spores released from bryozoans are broadly infectious to salmonids (see Appendix S1: Table S1 Initial infection of fish involves extrasporogonic proliferation of the parasite in the kidney. Clinical PKD is caused when the myxozoan parasite T. bryosalmonae elicits an inflammatory response during this proliferation in the kidneys of salmonid fishes, resulting in renal swelling, accompanied by severe kidney pathology and anemia (Schmidt-Posthaus et al. 2012). The development of clinical PKD is temperature dependent-in studies of rainbow and brown trout, infections are generally subclinical and mortality is low when water temperatures remain under 15°C (Okamura et al. 2011), but such thresholds are not well described for other salmonid species, such as mountain whitefish. Warmer temperatures and potentially linked multiple stressors (e.g., reduced water flow/ levels, eutrophication) may exacerbate infection severity leading to clinical PKD (Schmidt-Posthaus and Wahli 2015). In addition, infection dynamics or disease severity may vary depending on the host-parasite coevolutionary history. For example, PKD can be highly problematic in Europe, when farmed, non-native rainbow trout are infected by local strains of T. bryosalmonae (Morris and Adams 2006b). Interestingly, these non-native rainbow trout and other exotic fish hosts do not support spore production and may represent dead-end hosts for the parasite Adams 2006a, Grabner andEl-Matbouli 2008). In permissive fish hosts, parasite stages migrate to kidney tubules where spores infectious to bryozoans are produced.
The first documented case of PKD in North America occurred at the Hagerman State Fish Hatchery in Idaho, USA in 1981. Since that time, cases have also been reported in Washington, California, Alaska, and Montana, USA and British Columbia, Canada (Hedrick et al. 1993, Kent et al. 1995, Gorgoglione et al. 2020). The first recorded cases of PKD affecting free-ranging fish in North America occurred in a population of Yellowstone cutthroat trout in Middle Creek Reservoir in Montana during the summers of 1990 and 1991 (MacConnell and Peterson 1992). Another documented case of free-ranging fish infection soon followed in Chinook salmon (Oncorhynchus tshawytscha) in British Columbia (Kent et al. 1995). T. bryosalmonae was again observed in Montana in 1995, this time in rainbow trout sampled from Cherry Creek, a tributary to the Madison River (Bozeman Fish Health Center, unpublished data). Proliferative kidney disease-related mortalities were reported in mountain whitefish in the South Fork Snake River in Idaho in 2012 and in both the South Fork Snake and Yellowstone Rivers in 2016. Outbreaks of PKD have been documented in other wild fish populations in Norway (Sterud et al. 2007), Iceland (Kristmundsson et al. 2010), and Idaho (Bruun 2015). Declines in wild fish populations in Switzerland (Borsuk et al. 2006) and Germany (Arndt et al. 2019) have been linked with PKD-induced mortality. It is becoming increasingly clear that T. bryosalmonae infections are widespread, with sampling confirming infection in wild fish populations across Europe (Lewisch et al. 2018) and western North America, and new cases continue to be reported (Gorgoglione et al. 2020).
To reduce stress imposed on fish in the Yellowstone River during the 2016 PKD outbreak and to prevent the potential spread of T. bryosalmonae to other water bodies, Montana FWP closed a total of 295 km of the Yellowstone River and its tributaries to all recreational activities for several weeks. This closure resulted in the loss of an estimated USD$500,000 to local businesses, which comprise a thriving river-based recreational economy (Sage 2016). As mortalities subsided, the river was re-opened section by section beginning in early September until no more mortalities were observed during surveys in late September (Opitz and Rhoten 2017). In August of 2017, anglers again reported mountain whitefish mortalities. Subsequent Montana FWP initiated surveys documented 146 mountain whitefish mortalities near Livingston, Montana (Opitz 2017).
Several questions loomed in the wake of the 2016 and 2017 PKD outbreaks on the Yellowstone River. What is the distribution of T. bryosalmonae in Montana? Were the mortalities the result of a recent introduction of T. bryosalmonae into a na€ ıve population, or was the parasite endemic and the outbreak is explained by environmental conditions alone? Why were mortalities almost exclusively associated with mountain whitefish? This unexpected outbreak of PKD in the Yellowstone River underscores a lack of understanding of the ecology and epidemiology of the causative parasite in wild fish and in bryozoan hosts.
The aim of our work following the 2016 and 2017 Yellowstone River PKD outbreaks was to gain ecological and epidemiological insights on T. bryosalmonae in the Yellowstone River and the broader region, which has many local economies that rely on healthy fisheries. Here, we use population genetic studies, pathological investigations, transmission experiments, and environmental DNA (eDNA) detection to assess the following hypotheses: (1) T. bryosalmonae has recently invaded the Yellowstone River and encountered a na€ ıve host population, (2) a new genetic variant of T. bryosalmonae was responsible for the outbreaks, and (3) outbreaks were driven entirely by extreme environmental conditions.

Prevalence and distribution of T. bryosalmonae in Montana
To characterize the current distribution of T. bryosalmonae in Montana, we first used molecular diagnostic methods to screen wild-caught salmonid kidneys for presence of parasite DNA (non-salmonid fish were also occasionally tested as tissue-negative controls; Appendix S1: Table S1). This included kidney tissue collected as part of several large sampling efforts enabled by collaboration between state, federal, and international partners: (1) ad hoc archived samples from the Yellowstone River from 2011 and 2012, (2) ad hoc samples collected intermittently from 2016 to 2019 that were analyzed at the BFHC according to state and federal agency needs, (3) Montana FWP collections from the Yellowstone River and surrounding area in response to the 2016 PKD outbreak, and (4) a follow-up effort by Montana FWP in 2017 to characterize T. bryosalmonae distribution in 15 water bodies that are considered ecologically and economically important trout waters to Montana. Samples from 2016 through 2019 were collected either by electrofishing or angling. Yellowstone River samples from 2011 and 2012 were extracted from archived histology blocks of kidney tissue preserved in either Dietrich's or v www.esajournals.org Davidson's solutions. For each sample, a single 30-µm scroll was cut from a 4-cm bedding block and DNA was extracted with the same method used for fresh tissue.
In total, 1256 wild fish (summarized in Appendix S1: Table S1) and 238 histological samples were screened via endpoint PCR and/or quantitative (q)PCR for T. bryosalmonae between August of 2016 and September of 2019 (Table 1). Infections were inferred in endpoint PCR with visible bands using gel electrophoresis and in qPCR reactions that reach threshold fluorescence prior to cycle 40. Total lengths for sampled fish were not recorded, but, anecdotally, fish were generally >150 mm in length. Kidney tissue was removed from each fish, homogenized in sterile bags, and preserved at À20°C. DNA extractions of fish kidneys were carried out using Qiagen DNeasy Blood and Tissue Kits following the manufacturer's instructions (Qiagen.com, catalog 69506). An extraction blank (EB), with no tissue, was included in each extraction run to detect possible contamination from extraction equipment and reagents.
Histology samples that amplified T. bryosalmonae DNA and an equal number of putative negative samples were fixed and examined microscopically.
Paraffin-embedded kidney tissues were sectioned twice, each at a thickness of 5 µm. Paired sections were stained with hematoxylin and eosin (H&E) or Giemsa, respectively. H&E-stained kidney sections were scanned completely at 4009 magnification. Matching Giemsastained sections were examined at 1009 power to screen for areas of interest. When areas of interest were identified, those areas were examined under 4009 power as well.

T. bryosalmonae phylogeography and genetic diversity
To assess the phylogeographic structure of T. bryosalmonae in Montana and gain insights on whether the parasite may have recently invaded the system or might be a new strain, phylogeographic diversity assessments were made with internally transcribed spacer region I (ITS1), the cytochrome oxidase subunit I (COI), and a panel of simple sequence repeat (SSR) markers on a subset of screened fish tissue samples. We considered the use of multiple lines of evidence as imperative for drawing inferences from genes and genomic regions that have no known linkage to virulence. We selected 95 DNA extracts that were positive for T. bryosalmonae to gain insights on genetic variation associated with a diversity of salmonid species from across Montana (see Appendix S1: Tables S2-S4, for samples used in ITS1, COI, and SSR analyses, respectively).
Tetracapsuloides bryosalmonae-specific ITS1 primers were designed (Table 1) to make an initial assessment of parasite genetic variation within Montana waters. These primers show no crossreactivity with fish host DNA and produced a single amplicon of~500 base pairs (bp). Primers were validated in silico using NCBI Primer Blast to ensure that no similar binding regions exist in other parasites of salmonid fishes. Primer Blast did not return any non-target matches within the entire nucleotide database for products between 100 and 1000 bp. Amplified product for sequencing was prepared using ExoSAP-IT PCR Product Cleanup Reagent (Thermofisher.com, Catalog 78201.1.ML), followed by PCR with a BigDye Terminator v3.1 Cycle Sequencing Kit (Thermof isher.com, Catalog #4337455), and precipitation with sodium Acetate and EDTA, all according to manufacturer instructions. Sequencing was carried out on an Applied Biosystems 3500 Genetic Analyzer (Applied Biosystems, Foster City, California, USA). Samples with a high proportion of ambiguities were discarded prior to alignment. A total of 102 new sequences generated as part of this work were aligned with ClustalW using Bioedit (mbio.ncsu.edu/BioEdit/bioedit) and then manually edited. Consensus trees were created from posterior distributions generated using MrBayes v3.1.2 (nbisweden.github.io/ MrBayes). The trees for ITS1 sequences are unrooted, as no sufficiently similar outgroup sequence is currently available. A Metropolis Coupled Markov Chain Monte Carlo with six chains (one cold and five hot chains) was used to generate posterior probabilities after 1 9 10 6 generations, excluding the first 25% for burn-in.
v www.esajournals.org Modeltest v3.06 (evomics.org/resources/softwa re/molecular-evolution-software/modeltest) was used to select the general time-reversible model of nucleotide evolution with a gamma distribution (GTR+G). Graphics for phylogenetic trees were generated using the ggtree package (Yu et al. 2017) for R v3.6 (R Core Team 2014). In addition to new sequences generated as part of this work, T. bryosalmonae ITS1 sequences were also obtained from NCBI's nucleotide database. Sequences were trimmed to a homologous region of~500 bp prior to alignment and tree generation as described above. Our methods differ from a previous analysis of ITS1 diversity in T. bryosalmonae in North America and Europe (Henderson and Okamura 2004) in three notable ways. First, we did not clone sequences into plasmids and, thus, the influence of coinfecting strains on subsequent phylogenetic analysis is not accounted for. Second, we did not have the appropriate source material to assess intragenomic variation. However, Henderson and Okamura (2004) report intragenomic variation that is one order of magnitude less than intergenomic variation in T. bryosalmonae. Lastly, Henderson and Okamura (2004) weighted polymorphisms associated with SSRs to reflect single mutational events. However, the weighted and unweighted results of Henderson and Okamura (2004) were nearly identical. In our analysis, phylogenetic trees were constructed using unweighted nucleotide sequences, which yields slightly higher per nucleotide divergence estimates than weighting SSRs. (Table 1. Continued.)

Assay by application
Sequence Sequencing was also performed on the cytochrome oxidase I (COI) gene. Nested PCR was used to first amplify a 1320 bp region of the COI gene using primers Tbryo_COI_F28 and Tbryo_-COI_R1348, followed by an additional PCR using Tbryo_COI_F729 with the same reverse primer (Table 1). The same procedures used for sequencing, aligning, and phylogenetic analysis of ITS1 sequences were also used for COI sequences. Only 71 sequences were of sufficiently high quality for alignment and tree construction. COI trees are rooted with Buddenbrockia plumatellae as the outgroup.
Sixteen SSR markers were tested in pilot trials and 12 that amplified most consistently across preliminary test extracts were selected for use in screening 97 kidney extracts that had already tested positive for T. bryosalmonae DNA. SSR primers and cycling conditions can be found in Table 1. Each SSR assay was amplified in singleplex with untailed, fluorescently labeled forward primers and 5 0 -GTTT tailed reverse primers. Multiple SSR amplicons for the same sample were then combined to form three panels of four SSR assays that were analyzed with a LIZ600 sizing standard in multiplex on an Applied Biosystems 3130xl (Applied Biosystems) at the Idaho State University's Molecular Research Core Facility.
Genetic diversity estimates were made for several geographic scales using average nucleotide diversity (p) calculated in the pegas package (Paradis 2010) for R v3.6 (R Core Team 2014). Average nucleotide diversity is the number of pair-wise polymorphisms per nucleotide for all sequences in a defined group. The analysis of molecular variance (AMOVA; Excoffier et al. 1992) method and the / ST statistic were used to measure genetic differentiation of DNA sequences for ITS1 and COI genes using the poppr package (Kamvar et al. 2014). The / ST statistic is a representation of the proportion of DNA sequence diversity that can be attributed to differences between groups of sequences. SSR diversity was analyzed visually with discriminant analysis of principle components (DAPC) using the adegenet package (Jombart 2008).

eDNA sampling
We established an eDNA surveillance campaign focusing on southwestern Montana to increase understanding of spatiotemporal patterns relating to parasite distribution within two water bodies: the Yellowstone and Madison Rivers (Appendix S1: Table S5). During 2018, water samples were collected weekly from six sites spanning 64 km of the Yellowstone River ( Fig. 1) from April 20th to September 25th. During 2019, weekly water samples were collected from four of the previous six sites on the Yellowstone River and three sites spanning 32 km of the Madison River from June 4th to September 18th. At each site, six 250-mL replicate environmental water samples were collected using sterile Whirl-Pak bags. A field blank of reverse osmosis water was poured prior to sampling at each site and transported along with environmental samples. Samples were filtered in the laboratory through 4.5 cm diameter, 1.5-µm pore nitrocellulose filters using sterile, single-use polypropylene filter cups, and a peristaltic pump. A filtration blank, which involved filtering 250 mL of reverse osmosis water according to the same procedures as environmental samples, was performed at the end of each filtering session to detect possible contamination from filtration equipment. Filters were extracted as described above for fish tissue except that lysis was carried out in Qiagen Investigator Lyse & Spin Baskets (Qiagen.com, catalog 19597). An EB, with no filter of any kind, was included in each extraction to detect possible contamination from extraction equipment and reagents. T. bryosalmonae eDNA detection was carried out using qPCR according to Hutchins et al. (2018b). Cycle 40 was used as the quantification cycle cutoff for a positive detection.

Bryozoan surveys
In 2018 and 2019, visual and molecular attempts to document the occurrence of bryozoans were carried out concurrently with eDNA sampling. Visual surveys were performed by one person at each sample site by haphazardly inspecting exposed hard substrates every three steps along a transect parallel to the stream bank for ten minutes. Survey transects were positioned downstream of eDNA sampling but were carried out by moving in an upstream direction. Surveys were restricted to water depths of less than~25 cm for the sake of safety and water clarity. An endpoint PCR assay designed to amplify an approximately 156 bp region of the 12s rRNA gene was used to screen for bryozoans in each eDNA sample collected in 2018 (PLFR_12s_001f: AGTACTAGGGGAACCTGTG and PLFR_12s_ 156r: GCTGCACCTTGACCTG). This assay was designed specifically to amplify all species in the Plumatella and Fredericela genera for which sequence data were publicly available because these are the two groups most often associated with PKD outbreaks (Okamura and Wood 2002). This assay was validated in silico by ensuring that all plumatellids and fredericellids sequences had no mismatches with the primers and that there were no non-bryozoan species with fewer than three total mismatches, but the primer set will also amplify species from other genera. Specificity down to the two target genera was not essential because inclusivity within Phylactoleamata is beneficial for our intended exploratory use. In vitro validation was performed on tissue samples of cultured Plumatella vaihiriae and Fredericella indica (Bryo Technologies, Beavercreek, Ohio, USA) and on tissue from wild-growing Plumatella rugosa and F. indica collected from the Columbia River, near Grand Coulee, Washington, USA. The cycling conditions used for the PLFR assay were identical to that of the PKX18s assay (Table 1). Amplicon presence was inferred using gel electrophoresis with SYBR green stain.

Sentinel fish exposures
Na€ ıve, disease-free sentinel fish exposures were employed in the Yellowstone and Madison rivers to correlate presence of infectious waterborne T. bryosalmonae spores with patterns of eDNA parasite detections. Sentinel cage locations coincided with eDNA sample locations. Two galvanized steel-mesh fish cages were placed at each of four sample sites along the Yellowstone River (2018 and 2019) and at three sample sites along the Madison River (2019 only). Cages were affixed to the substrate with steel rebar and positioned such that at least 50% of the 120-L cages were submerged in an area where the velocity of moving water was less than 0.25 m/s. As the water level in the rivers lowered over the course of the experiment, cages were repositioned as needed at each site to adhere to the criteria above. Disease-free Arlee (2018) and Fish Lake strains (2019) of rainbow trout were acquired from the U.S. Fish and Wildlife Service's Ennis National Fish Hatchery for use in the sentinel cage exposures. Methods of animal handling were approved by the Institutional Animal Care and Use Committee of Montana State University (permit #2019-108). While it would have been ideal to use mountain whitefish for the study, these fish are notoriously fragile in rearing environments, and there are no established diseasefree brood stocks available. Rainbow trout were chosen because they have previously demonstrated infection and mortality from T. bryosalmonae (Grabner and El-Matbouli 2009), are an important model species for disease research, have been introduced across the world, and account for 18% of all global salmonid aquaculture production (FAO 2020). Sentinel fish total lengths in 2018 varied between~100 and 250 mm at the time of exposure, while in 2019, fish lengths varied between~100 and 150 mm.
On the first day of the exposure study (mid-July) and every three weeks thereafter, 20 fish were placed in each cage. After three-week exposures, all fish were collected from cages at each site for T. bryosalmonae screening. New fish were then loaded for the next three-week exposure. Four rounds of exposures were performed, and the study was terminated in late September. Caged fish collections occurred immediately following upstream eDNA sampling at the cage site. Upon collection, fish were euthanized with an overdose of pH-buffered MS-222 (300 mg/L tricaine methanesulfate, SigmaAldrich.com, catalog #E10521). Collected fish were transported back to the lab and frozen at À20°C until further processing. Gill and kidney tissue were removed from each sampled fish and DNA extracted as described above. All kidney tissues and a subset of gill tissues were screened for T. bryosalmonae using qPCR as described above. A negative control cage was maintained for a single three-week exposure at the Ennis National Fish Hatchery in both 2018 and 2019 whereas positive control fish were held for two months in the summer of 2019 at a location in the Columbia River (Washington, USA) known to have annual PKD outbreaks.

T. bryosalmonae infections in Montana
We detected T. bryosalmonae DNA in 501 of the 1214 salmonid fish that were screened, indicating an infection prevalence of 41%. T. bryosalmonae had only been documented twice before in Montana (MacConnell and Peterson 1992; USFWS Bozeman Fish Health Center, unpublished data), although it has never been routinely monitored for in Montana. While one of those historical detections was just 70 km west of the Yellowstone River (in a tributary to the Madison River), hydrologically, both observations were more than 1600 river km apart. Because of the temporal and spatial disparity of these observations, they were generally regarded as rare and isolated instances. This notion was what ultimately lead to speculation that the 2016 Yellowstone River PKD outbreak may have been due to a recent introduction. But our data paint a different picture, as T. bryosalmonae infections were present in every watershed in which sampling occurred in 2016 and 2017 (Fig. 2). The highest prevalence of T. bryosalmonae infection was in mountain whitefish (50%, n = 541), and the lowest in Yellowstone cutthroat trout (24%, n = 127), but these prevalence estimates were biased by the location and timing of sampling. For example, many mountain whitefish samples were collected in the Yellowstone River during the outbreak, and a large proportion of Yellowstone cutthroat trout were collected from the Lamar River in Yellowstone National Park. Infection prevalence estimates for brown trout, brook trout, and rainbow trout were 35% (n = 257), 38% (n = 34), and 39% (n = 255), respectively. Given that sample collection sites were not randomly selected water bodies, we cannot make general inferences regarding the hostspecific prevalence of T. bryosalmonae, but our results do confirm a broad host range and spatial distribution in the studied area. T. bryosalmonae DNA was also detected by PCR in six of the 238 archived histology samples from 2011 to 2012. Objects meeting the diagnostic criteria of T. bryosalmonae were not detected in the sections of either the PCR-positive or PCR-negative fish by histology. However, PCR is likely a more sensitive diagnostic tool and the material sampled for histology (10 µm thickness in total) was approximately one third that of what was sampled for PCR, so the discrepancy between the two diagnostic methods is not altogether surprising.

T. bryosalmonae genetic variation
The average nucleotide diversity of ITS1 sequences was highest at the global scale (p = 0.014) and similar for water bodies across Montana for both ITS1 and COI (p = 0.006 for each). However, compared with average nucleotide diversity from the Yellowstone River (ITS1 p = 0.010, COI p = 0.008), nucleotide diversities associated with dead or moribund fish were lower for ITS1 (p = 0.006) and identical for COI (p = 0). The ITS1 and COI genes and the SSR panel exhibit little evidence of sorting into distinct geographic clades for our samples at the scale of the water body (ITS1 / ST = 0.060, COI / ST = 0.401, SSR / ST = 0.325; Fig. 3A, B) but continental-scale groups accounted for a greater proportion of observed genetic variation in ITS1 and COI (ITS1 / ST = 0.100, COI / ST = 0.600; Fig. 3C, D). The North American parasite ITS1 sequences sorted into a single clade, that includes parasite sequences isolated from Italian and French populations (Fig. 3C); this pattern is identical to the findings of Henderson and Okamura (2004). This general pattern that European parasite diversity appears to be a subset of North American diversity is mirrored in the COI tree, where sequences from the UK, Germany, and France form a monophyletic group within the USA sequences (Fig. 3D). There was no evidence that sequence similarity was related to the host species (data not shown). Poor sorting of the samples may be due to any combination of high gene flow, very recent geographic isolation, low mutation rates, or, in the case of ITS1, the confounding effects of high intragenomic variation (Henderson and Okamura 2004).
None of the 12 SSR markers amplified across all 97 samples, and six of the markers showed no samples that amplified (Appendix S1: Table S4). After removing nontyped individuals (those which did not amplify for any SSR marker) and invariant loci, six markers and 72 samples remained. DAPC was used to identify clusters of genetically similar individuals from SSR data based on the river from which the sample was collected (Fig. 4A) and the condition of the fish host (Fig. 4B). As with the ITS1 sequence data, there was no evidence of clustering according to river of origin. However, distinct clusters distinguished samples taken from healthy fish vs. mortality/moribund fish suggesting a genetic difference for parasite strains associated with fish mortalities in the Yellowstone River. Larger sample sizes would strengthen inferences on genetic variation in relation to river basin and fish host condition.
The population genetics tools do not point toward geographic differentiation of the parasite. Taken in isolation, this low diversity could be the result of recent and extremely rapid spread of the parasite throughout the region. However, when taken in concert with other data collected as part of this study, recent and rapid expansion is an unlikely scenario, which is expounded upon further below. Rather, the SSR and, to a lesser degree, COI evidence suggest that parasites that caused mortalities may have been differentiated from those found in apparently healthy fish both during and after the 2016 mortality event. This interpretation, however, is based on a small sample of temporally correlated mortality samples assessed with DNA sequences that, as far as we know, do not influence virulence. Until a greater number of mortalities can be collected and/or functional genomics methods can be deployed, alternative explanations should be considered.

Environmental conditions
If strong evidence that an introduced species or strain caused the Yellowstone River mortality events is lacking, how well do environmental conditions, alone, explain the events? Under this hypothesis, we would expect environmental conditions in the Yellowstone River in 2016 to be anomalous in both time and space. It is well established that PKD is a temperature-related disease in salmonid fish, with clinical disease signs often appearing above 15°C for cool-water trout species (Bettge et al. 2009). While in 2016, there were 70 summer (June-August) days above 15°C recorded at the Yellowstone River USGS streamgage (06192500) near Livingston, Montana, the average AE standard error from 2000 to 2019 was just 59.8 AE 2.1 d. However, similar or hotter conditions have been observed in previous years in the Yellowstone River (Fig. 5A) and in other southwestern Montana rivers where T. bryosalmonae has been documented in 2016 (Fig. 5B). River discharge may also influence the extent and severity of PKD through density-dependent interactions of fish proximity to infected bryozoans and fish stress (Bruneaux et al. 2017). Lower flows effectively reduce habitat and limit resources for competing fish (Hakala and Hartman 2004) and may simultaneously increase time periods when fish linger in the vicinity of bryozoans. Low flow conditions, in terms of mean daily, peak, and cumulative discharge, have occurred in the Yellowstone River multiple times in the recent past (Fig. 6). Thus, it would appear that, prior to 2016, PKD outbreaks were not provoked despite the occurrence of T. bryosalmonae infections in fish (confirmed by sampling archived tissues) and environmental conditions coinciding with those of the 2016 outbreak. This presence in fish outside of mortality events suggests that an additional variable or dimension is required to trigger an outbreak, such as multi-year effects or complex interactions involving hosts, parasites, and environmental conditions. For instance, it is noteworthy that 2016 was the second consecutive year of warm and low-flow hydrological conditions in the region, suggesting interplay between outbreaks and multi-year conditions. Alternatively, previous outbreaks may have occurred undetected. This propelled us to examine associations between environmental parameters (water temperature, flow, substrate, etc.) and the probability of detecting T. bryosalmonae in water and in sentinel fish as described below.

Detection of eDNA
Tetracapsuloides bryosalmonae DNA was first detected as early as April in 2018, when sampling commenced ( Fig. 7; Appendix S1: Table S5). T. bryosalmonae detections in at least one field v www.esajournals.org replicate occurred most frequently at the most downstream site (Carter's Bridge) and least frequently at the most upstream site (Brogan's Landing) of the Yellowstone River. Detection at the latter site was somewhat surprising as it is separated from the area of high fish mortality in 2016 by a section of high gradient rapids. In 2018, detections were more frequent at the lower three sites on the Yellowstone River, but not in 2019. Of the 24 weeks spanning the sampling period, there were seven occasions when T. bryosalmonae was not detected at any sample site. In 2019, sampling did not begin until June, and the first detections for both the Madison and Yellowstone River occurred on June 11th (Fig. 7). Across both years in the Yellowstone River, peak T. bryosalmonae detection occurred in late June or early July.
T. bryosalmonae DNA detection was too sporadic through space and time to infer any associations between environmental parameters (water temperature, depth, velocity, and substratum type) and detection (data not shown).
Our data do not distinguish T. bryosalmonae DNA in spores infective to bryozoans and fish thus complicating the significance of eDNA detections in river water samples with respect to fish infection. Previous investigations do, however, provide evidence that T. bryosalmonae spores are released from bryozoans in the UK largely during late spring/early summer and again in autumn (Tops et al. 2004). The timings of first spore detection (April) and peak spore concentrations in water (June/July) in the Yellowstone River broadly agree with the spring/summer production of infective spores in bryozoans. The occurrence and identity of bryozoan hosts of T. bryosalmonae in the Yellowstone River require investigation. Yearly variation in temperature and flow regimes may influence the amount and timing of infectious spore release from bryozoans that may vary in identity, abundance, and distributions, contributing to the myriad of potential factors leading to an outbreak.

Bryozoan surveys
A total of 204 visual surveys (six sites and 34 events) and 1224 endpoint PCR reactions (six sites, 34 events, and six field replicates) were performed for the Yellowstone River in 2018. In 2019, 64 visual surveys were performed on the Yellowstone River (four sites and 16 events) and 48 times on the Madison River (three sites and 16 events). No endpoint PCRs for bryozoans were performed in 2019. Neither visual surveys nor endpoint PCR revealed the presence of bryozoans in any of the eDNA samples at and any time or sampling location (data not shown). Bryozoans are obligate hosts of T. bryosalmonae, and so must have been present in the river at some point before the 2016 outbreak. Possible explanations for our inability to find physical evidence of bryozoans beyond their true absence from the system are numerous and include, while not being limited to, inappropriate sample/survey locations or methods. The locations were limited by logistical constraints of rapidly moving water and by our ability to thoroughly inspect only that which we could either lift out of the water or clearly see under the water. These locations also varied based on drastically different water levels based on the river stage, and such variable conditions are not conducive to large colony formation by bryozoans. eDNA surveys may not be an effective tool in bryozoan surveillance because they may simply not shed enough DNA or disperse it widely enough to be easily detected. Lastly, because we do not know with certainty what bryozoan species inhabit waters of the region, it is possible that the species responsible for hosting T. bryosalmonae in these waters is not amplified by the assay designed for this study.

Sentinel fish study
There were no detections of T. bryosalmonae in any gill or kidney tissue sampled from caged fish except for three fish held at the known positive control site. These results are consistent with the lack of wild fish mortalities during the study. Notably, hydrologic and thermal conditions during the two years in which cages were deployed (2018 and 2019), starkly contrast with conditions during 2016 and 2017 when PKD mortalities occurred. Given the strong link between PKD prevalence and severity and thermal conditions, it is therefore crucial to interpret our results in the context of the conditions that fish hosts experience. In general, our results provide further field-based evidence that environmental conditions contribute to PKD outbreaks.
The lack of infection in sentinel fish is not, however, consistent with the geographically broad 41% infection prevalence that we found from screening wild fish as described above. Furthermore, we know from concurrent sampling that infected wild fish were present in the study area during the sentinel fish study period. However, fish malacospores are only infective to bryozoans, so direct fish-to-fish transmission is not possible (D'Silva et al. 1984). Non-infective spores released from fish could explain the simultaneous eDNA detection of T. bryosalmonae near cages but not in sentinel fish. Alternatively, spores may have difficulty in infecting fish, especially if they are present in low concentrations. Wild fish likely move more broadly than the caged fish and may preferentially use habitats where they come into closer contact with bryozoans, such as shelter near large boulders, crevices, riparian tree roots, deadwood, and banks of aquatic vegetation. The concentrations of infective T. bryosalmonae spores between habitats may vary, resulting in different likelihoods of exposing wild-ranging vs. caged fish. This inference is supported by fact that thousands of T. bryosalmonae spores develop within multiple parasite sacs inside bryozoan colonies (Okamura et al. 2011) and that spores may not disperse equally within a river network (Carraro et al. 2018).

Synthesis and future avenues of research
Screening of free-ranging fish has revealed widespread T. bryosalmonae infection throughout salmonid waters in Montana, whereas the parasite was not found during every sampling effort, it has been detected in every surveyed water body. Collectively, the widespread distribution, historical observations of T. bryosalmonae in other regional water bodies, and the lack of phylogeographic structure revealed by ITS1, COI, and microsatellite data suggest that the outbreak in 2016 did not result from a recent T. bryosalmonae introduction. This hypothesis is further refuted by the positive molecular detections in archived histology tissue of salmonid species from the Yellowstone River collected in 2012. Henderson and Okamura (2004) suggest that T. bryosalmonae was introduced to Europe from North America sometime since the most recent glaciation. This likely phylogeographic history and the widespread detection of T. bryosalmonae in western regions of North America, support a scenario of long coevolutionary history between host and parasite, rather than recent parasite introduction. Nevertheless, T. bryosalmonae's long-term presence in the Yellowstone River does not preclude the possibility that a novel and more virulent genetic strain of the parasite was the cause of the 2016 outbreak. Our SSR and ITS1 data suggest that this is a possibility in view of the distinct grouping of T. bryosalmonae samples that were collected from dead or moribund fish. However, given our low sample sizes, this correlation between fish condition and T. bryosalmonae genetic diversity requires further investigation. Future genomic, transcriptomic, and functional v www.esajournals.org analyses will be crucial to determine the T. bryosalmonae virulence factors and disease dynamics.
Though we cannot exclude the possibility of host-specific virulence factors, the relatively high mountain whitefish mortalities indicate that this species may be particularly valuable as an indicator of temperature-related disease. Measures of thermal tolerances, such as critical thermal maximum (Becker and Genoway 1979), of mountain whitefish (26.7°C; Brinkman et al. 2013) are generally lower than for other salmonid species found in the Yellowstone River: rainbow trout (29.1°C; Currie et al. 1998), brown trout (29.8°C;Elliott and Elliott 1995), and cutthroat trout (28.1°C; Wagner et al. 2001). Further work is required to determine whether patterns of mortality may exemplify species-specific thresholds in susceptibility to temperature-related diseases like PKD.
The upper Yellowstone River watershed has experienced several years that were both hotter and dryer than the conditions experienced in the outbreak years. However, 2016 was the only warm and low-flow year in the past two decades of this region that immediately followed another year with abnormally high temperatures and low discharge (Figs. 5, 6). Multi-year hydrological conditions such as these are rarer than single-year instances and may play a vital role in cultivating the conditions that can produce PKD outbreaks. Multi-year environmental conditions of warm summers and low water flows may simultaneously influence the abundance, distribution, and infection state of bryozoans. Higher temperatures are associated with increased growth of bryozoans and higher intensity of infections in bryozoans in culture (Tops et al. 2009). Consecutive years favoring bryozoan development may thus increase infective spore concentrations in the river and demonstrate as mass outbreaks. Clinical PKD in Europe is predominantly observed in young-of-the-year trout, however, during the 2016 outbreak, whitefish of different age classes were affected (although no data were collected on this specifically). Rainbow trout in +1 yr age classes previously exposed to T. bryosalmonae either avoid infection or show earlier and more effective immune response (Foott andHedrick 1987, Bailey et al. 2017). It is not known whether only na€ ıve fish succumbed during the Yellowstone PKD outbreaks, or if mortalities also occurred amongst previously exposed fish. Resilience of populations to PKD outbreaks may thus jointly depend on previous exposure history and environmental factors that can overwhelm the natural immunity (e.g., increased stress for longer periods of time during abnormally hot/dry years). Various data thus suggest that sustained environmental conditions that promote infectious spore release and disease development in fish may be the most important drivers of infection prevalence and PKD outbreaks.
Greater understanding of the ecology of bryozoans and their interactions with T. bryosalmonae will reveal drivers of PKD epidemiology in Montana rivers systems, as bryozoans are, ultimately, the source of the disease. Notably, bryozoans in southwestern Montana have proven to be elusive. In spite of weekly efforts to detect bryozoans with visual surveys and molecular tools (e.g., eDNA sampling), we were unable to locate or detect bryozoans at sampling sites along the Yellowstone or Madison Rivers. A potential explanation for this is that bryozoan populations may themselves have been challenged in some unknown way in the years following the outbreaks. Intriguingly, a recent study has demonstrated that poor bryozoan host condition provokes T. bryosalmonae to produce spores for horizontal transmission to fish (Fontes et al. 2017a). This form of terminal investment may provide a means of bailing out of hosts on the precipice of death and could simultaneously explain the sudden fish outbreaks and the subsequent difficulty in finding bryozoans in the system.
Our current hypotheses for the causes of the PKD mortality events in the Yellowstone River based on the insights gained in this study are (1) multi-year conditions of low discharge and high temperature promote the growth of endemic infected bryozoan populations that would otherwise remain at low densities and simultaneously increase T. bryosalmonae proliferation within bryozoan hosts due to warming of river waters during warm and low-flow years; (2) warmer water temperatures (≥20°C) promote T. bryosalmonae spore release during summer months; (3) lower and warmer water conditions constrain fish to deeper pools where they will be subject to sustained exposure to relatively higher concentrations v www.esajournals.org of parasite spores; (4) mountain whitefish are particularly sensitive to temperature-dependent diseases because of a lower thermal tolerance than other salmonid species in the Yellowstone River.