Continent-wide recent emergence of a global pathogen in African amphibians

Introduction Emerging infectious diseases are increasingly recognized as a global threat to wildlife. Pandemics in amphibians, caused by the fungal pathogen Batrachochytrium dendrobatidis (Bd), have resulted in biodiversity loss at a global scale. Genomic data suggest a complex evolutionary history of Bd lineages that vary in pathogenicity. Africa harbors a significant proportion of global amphibian biodiversity, and multiple Bd lineages are known to occur there; yet, despite the decline of many host species, there are currently no described Bd-epizootics. Here, we describe the historical and recent biogeographical spread of Bd and assess its risk to amphibians across the continent of Africa. Methods We provide a 165-year view of host-pathogen interactions by (i) employing a Bd assay to test 4,623 specimens (collected 1908–2013); (ii) compiling 12,297 published Bd records (collected 1852–2017); (iii) comparing the frequency of Bd-infected amphibians through time by both country and region; (iv) genotyping Bd lineages; (v) histologically identifying evidence of chytridiomycosis, and (vi) using a habitat suitability model to assess future Bd risk. Results We found a pattern of Bd emergence beginning largely at the turn of the century. From 1852–1999, we found low Bd prevalence (3.2% overall) and limited geographic spread, but after 2000 we documented a sharp increase in prevalence (18.7% overall), wider geographic spread, and multiple Bd lineages that may be responsible for emergence in different regions. We found that Bd risk to amphibians was highest in much of eastern, central, and western Africa. Discussion Our study documents a largely overlooked yet significant increase in a fungal pathogen that could pose a threat to amphibians across an entire continent. We emphasize the need to bridge historical and contemporary datasets to better describe and predict host-pathogen dynamics over larger temporal scales.


Introduction
Amphibians are one of the most vulnerable groups of vertebrates with 41% of all species currently threatened with extinction (Stuart et al., 2004;Wake & Vredenburg, 2008;The IUCN Red List of Threatened Species, 2021). An emerging infectious disease, chytridiomycosis, caused by the fungal pathogen Batrachochytrium dendrobatidis (Bd), has been a major driver of amphibian declines (Skerratt et al., 2007;Wake and Vredenburg, 2008;Fisher et al., 2009). Invasion and spread of the Bd pathogen are linked with declines and extinctions of hundreds of amphibian species worldwide and Bd has been documented on every continent where amphibians occur (Skerratt et al., 2007;Fisher et al., 2009;Olson et al., 2013;Scheele et al., 2019;Fisher and Garner, 2020). However, the pathogen's specific historical distribution and geographic spread remain unclear, hampering efforts to predict which amphibian species are at greatest risk of disease outbreak (Berger et al., 2016).
While effects of Bd on vertebrates are greater than those of any other observed pathogen (Scheele et al., 2019;Fisher and Garner, 2020), there is still much that is not understood. For example, while Bd-caused epizootics resulting in amphibian species decline or extirpation have been documented in detail in various parts of the world (Lips et al., 2006;Fisher et al., 2009;Vredenburg et al., 2010;Walker et al., 2010), it is still not known why epizootics are found in some areas but not others (Berger et al., 2016). Delineating the geographic origin and spread of the pathogen, as well as its evolutionary history and population genetic structure, are necessary to describe the dynamics of this pathogen-host system. Recent studies reveal a complicated geographic and phylogenetic history for Bd, including multiple lineages that vary in pathogenicity and may have spread across and between continents over a period of hundreds or even thousands of years. The genetic structure of Bd populations suggests that competition and hybridization among and between Bd lineages has occurred (Farrer et al., 2011;Schloegel et al., 2012;Bataille et al., 2013;Rosenblum et al., 2013;O'Hanlon et al., 2018). The global panzootic lineage (Bd-GPL; Farrer et al., 2011) is regarded as the most virulent lineage, has the widest geographic distribution, and is associated with epizootic events (Schloegel et al., 2012;O'Hanlon et al., 2018). Although Bd is believed to have originated in Asia, it is unclear when, where, or how Bd-GPL diverged from older lineages (O'Hanlon et al., 2018). Other lineages thought to be less virulent are associated with areas that generally lack epizootics (Berger et al., 2016;O'Hanlon et al., 2018), including several identified in Asia (Bd-ASIA-1, believed to have originated in Korea, which includes Bd-CH, originally thought to have originated in Switzerland), South America (Bd-ASIA-2/Bd-BRAZIL) and southern Africa (Bd-CAPE, first described from Cape Province, South Africa) (Farrer et al., 2011;Bataille et al., 2013;Rodriguez et al., 2014;O'Hanlon et al., 2018).
Africa harbors nearly 16% of known living species of amphibians (AmphibiaWeb, 2021) and multiple lineages of Bd have been found to infect these hosts in the wild (including Bd-CAPE and Bd-GPL). There are few relevant Bd host studies in Africa and no described Bdepizootics in the region (Berger et al., 2016;O'Hanlon et al., 2018;Doherty-Bone et al., 2020;Zimkus et al., 2020). However, this is likely due to a lower Bd sampling effort in Africa relative to other continents rather than a true absence of such events. Several studies have reported a small number of Bd-infected amphibians collected in the early 20 th century in Africa (Weldon et al., 2004;Soto-Azat et al., 2010). More recent African amphibian studies have revealed much higher Bd infection prevalence that may be associated with host declines (Goldberg et al., 2007;Kielgast et al., 2010;Bell et al., 2011;Gower et al., 2012;Tarrant et al., 2013;Seimon et al., 2015;Hirschfeld et al., 2016;Jongsma et al., 2016). Although these Bd-host infection data vary in spatial and temporal resolution, they suggest a pattern of Bd emergence that may represent a threat to African amphibian species, warranting a comprehensive analysis encompassing the broadest possible temporal and geographic scales.
We set out to address an important prediction regarding continent-wide infection dynamics of Bd: that the emergence of Bd in Africa is not limited to isolated outbreaks but instead represents a broad pattern of synchronized emergence across the entire continent. We address this prediction by comparing the historical frequency of Bd occurrence (spanning over a century) to the recent frequency of Bd occurrence (roughly the last 20 years). These comparisons of the frequency of Bd occurrence (revealing epizootic dynamics) are made both within African countries and within broader geographic regions. In addition, we include comparisons of changes in overall Bd infection intensity (zoospore counts) for specific countries in central Africa (Cameroon) and eastern Africa (Kenya) that are known to harbor widespread Bd in recent years (Hirschfeld et al., 2016;Kielgast et al., 2010). For our analyses, we have gathered 16,920 original and previously published museum-and field-collected records between 1852 and 2017 from a total of 36 countries. Importantly, our survey of Africa is unique in that it includes both positive and negative Bd samples as opposed to a more restricted reporting of only the number of Bd positives over time (Zimkus et al., 2020).
Additionally, if there is a continent-wide synchronous emergence of Bd, there could be regional differences in the magnitude of Bd prevalence. Such differences are likely associated with the distributions of the two described Bd lineages detected on the African continent (Bd-CAPE and Bd-GPL), depending on the presence of one or a combination of both in certain areas. Bd-CAPE has been detected in Cameroon and was found to be widely present in South Africa, while the Bd-GPL lineage has invaded almost every region of the continent (Rosenblum et al., 2013;O'Hanlon et al., 2018;Byrne et al., 2019). While it is speculated that the hypervirulent Bd-GPL is the more virulent of these two lineages based on ex situ experiments, less is known about the impacts of Bd-CAPE in nature, and there is evidence that Bd-CAPE does exhibit high virulence under certain conditions (Doddington et al., 2013;O'Hanlon et al., 2018;Fisher et al., 2021). Competition among invasive Bd strains, and the potential for hybridization has been documented on other continents (Jenkinson et al., 2016) and genotyping has revealed lineages in Africa that are intermediary to the more derived Bd-GPL complex and the more early diverging Bd-CAPE clade (Rosenblum et al., 2013;O'Hanlon et al., 2018;Byrne et al., 2019). These intermediaries are likely to be hybrids of the two described lineages and could be impacting disease dynamics (Farrer et al., 2011;Greenspan et al., 2018).
To understand how the dynamics of Bd lineage competition and/ or hybridization may lead to future epizootics in some areas and not in others, we utilize our large database to document changes in Bd prevalence in regions that have experienced invasion by different lineages. To predict the risk to future populations of amphibians across the continent, we develop a spatially explicit model of Bd risk across regions of Africa, based on our current understanding of the habitat suitability of these Bd lineages and amphibian species richness. Combined, this work provides novel predictions for the regions in Africa that are at high risk of Bd invasion and is complementary to a recent spatially explicit suitability model for the African continent (Zimkus et al., 2020).

Generating and compiling Bd records dataset
We sought to characterize historical and recent Bd prevalence in African amphibians. To do this we first compiled a dataset including (1) original records we generated (referred to as "our data") from amphibian skin swabs assayed for Bd and (2) published Bd records (referred to as "published data") for African amphibians.
For records we generated (our data), we collected skin swabs from 2,972 amphibian museum specimens collected from 1908-2009 in Cameroon, Ethiopia, Kenya, Lesotho, Tanzania and Uganda, and from 1,651 live field animals collected by coauthors and collaborators from Burundi and Equatorial Guinea (in 2011), and from Cameroon and the Democratic Republic of Congo (in 2011 and2013) for a combined total of 4,623 samples. We used a standard Taqman Bd-qPCR assay to determine the frequency of Bd occurrence and Bdinfection intensity (Boyle et al., 2004;Hyatt et al., 2007). This assay was designed to detect Bd-GPL, and therefore a limitation is that divergent lineages with mutations in the Taqman binding site may not amplify (Mutnale et al., 2018). Other retrospective studies utilized historical specimens (collected during a similar time-frame) and using the same Bd-qPCR assay, found evidence for Bd-infection even in their oldest specimens Rodriguez et al., 2014;Talley et al., 2015). This technique was shown to generate a recovery rate of 82% from museum specimens , and has robust power of detection, especially with large datasets. We focused historical sampling on areas with putative Bd-related declines [Lesotho (Minter et al., 2004), Ethiopia (Gower et al., 2012;Gower et al., 2013), Tanzania (Channing et al., 2006;Weldon et al., 2020)], areas with many available samples (Uganda and Kenya), and those with the most extensive recent sampling [Cameroon (Blackburn et al., 2010;Hirschfeld et al., 2016) and Kenya (Kielgast et al., 2010)]. Museum specimens are housed in the permanent collections of the California Academy of Sciences (CAS), the Museum of Vertebrate Zoology (MVZ) at the University of California, Berkeley, and Harvard University's Museum of Comparative Zoology (MCZ). Field sampling locations (live animals) were selected opportunistically based on field trips of coauthors and collaborators. As our goal was to assess Bd prevalence across all amphibians, we sampled amphibians from all taxonomic classes that were available in collections or captured at field sites (these included Anura and Gymnophiona).
Gloves were worn during swab sampling and changed between animals to prevent cross contamination. In addition, archived museum specimens were rinsed with 70% ethanol before swabbing to reduce cross contamination between animals housed in the same jar. We used standard skin swabbing procedures (which include 30 strokes of ventral surfaces) for Bd qPCR assays on live animals and museum specimens using Bd-specific ITS1 primers (Boyle et al., 2004;Hyatt et al., 2007;Cheng et al., 2011). Swabs were dried and stored in 1.5 ml microcentrifuge tubes at 4°C until DNA extraction using Prepman Ultra DNA Extraction Kit following the manufacturer's protocol (using 40 ml Prepman Ultra per sample). Standard real-time PCR quantification techniques were used: one replicate per sample was run with standard controls of 0, 0.1, 1, 10, and 100 zoospore equivalents (ZE) (Boyle et al., 2004;Hyatt et al., 2007;Cheng et al., 2011). Bd standards were created using Bd isolate CJB7 collected in 2009 from host Rana muscosa during an epizootic at Sixty Lake Basin in the southern Sierra Nevada, California , identified as part of the highly virulent Bd Global Panzootic Lineage (Bd-GPL1; James et al., 2015;Piovia-Scott et al., 2014;Rosenblum et al., 2013). We multiplied the qPCR genomic equivalents (GE) by a PCR dilution factor of 80 in order to quantify the number of zoospore equivalents on each swab (Z swab ; Briggs et al., 2010;Vredenburg et al., 2010). Any sample with a Z swab >0 was considered Bd-infected or Bdpositive. Because some animals had very high Z swab values, we report infection intensity on a base-10 logarithmic scale, represented with units of logZE.
For our compilation of published data, we obtained records from Bd surveys conducted in continental Africa (excluding Madagascar but including the islands of the Gulf of Guinea) from scientific literature and from the Bd-maps Legacy Project (Olson et al., 2021; dataset available at: https://geome-db.org/workbench/projectoverview?projectId=300). We included records from Anura, Caudata, and Gymnophiona from all life-stages with location data (i.e., country, latitude and longitude when available), collection year, and Bd infection status (i.e., presence or absence), and excluded food markets and the pet trade. While many published records used methods for swab sampling and Bd qPCR described above, we also included records using different sample materials (including skin scrapes/brush samples, toe clips, skin/mouthpart tissue, Bd isolates) and/or Bd screening methods (including conventional PCR, histology, microscopic visualization, Bd genome sequencing or genotyping). These published data include 12,297 Bd records (collected 1852-2017). When combined with our data (museum and field records we generated; N=4,623), we include a total of 16,920 Bd records from 36 African countries collected between 1852 to 2017 for this study. For more information on the dataset, see the Supplementary Material.

Bd emergence
To characterize the temporal pattern of Bd throughout Africa and estimate the timing of Bd emergence, we examined the frequency of Bd occurrence over time for all of Africa, for individual countries, and for regions (as defined by the United Nations). For frequency of Bd occurrence across all of Africa, we analyzed 16,857 records from our compiled dataset (summarized in Table 1). These included our data (all new records generated for this study; N=4,623) and 12,234 records from published data (we excluded 63 of the 12,297 published records in analyses of frequency of Bd occurrence because the associated studies did not report both Bd-negative and Bd-positive samples for Bd detection in amphibian hosts; these 63 Bd-  positive records were included, however, in our Bd risk analysis described below that only required Bd-positive samples with latitude and longitude data). To examine temporal changes in Bd occurrence throughout Africa, we binned records by decade (excluding an additional 171 records from published data for which collection year did not fall into a single decade; N=16,686; see Supplementary Material; Table S1) and calculated the frequency of Bd occurrence by dividing the number of Bd-positive records by the total number of animals tested. We calculated a 95% binomial confidence interval (95% CI) for each decade using the R function binom.test (R Core Team, 2022;Phillips and Puschendorf, 2013;De Leoń et al., 2017). We assessed our power to detect a prevalence <11%, a conservative estimate of Bd endemism based on a previous study that spanned over a 100 year period (Talley et al., 2015), by calculating the probability of detecting zero positives for each decade (defined as "Pr(no Bd)"). To do so, we used the R function dbinom (R Core Team, 2022;Doraj-Raj, 2014), which is based on the binomial distribution; we used 11% as the "true" probability that each sample was positive and the total number of samples evaluated in a particular time period as the number of trials.
Next, we examined the frequency of Bd occurrence and Bdinfection intensity by decade in countries for which we had the most comprehensive data: Cameroon and Kenya. For these countries, we analyzed a subset of records for which we had data on Bd infection intensity in addition to infection status (this included our data generated for this study for Cameroon and Kenya, data from Hirschfeld et al. (2016) and Zimkus et al. (2020) for Cameroon, and data from Kielgast et al. (2010) for Kenya; see Supplementary Material; Tables S2, S3). In addition to calculating frequency of Bd occurrence with 95% CI and Pr(no Bd) (as described above), mean and maximum Bd-infection intensity for Bd-positive samples were calculated and visualized by decade. We emphasize, however, that estimates of infection intensity from museum specimens are unreliable due to DNA degradation that occurs over time, so our historical infection intensities may be inaccurate estimates.
Finally, we defined a time of emergence based on the increase in frequency of Bd occurrence we observed across decades for all of Africa, Cameroon, and Kenya (year = 2000). We binned Bd records by whether they were pre-emergence ("historical," which we define as "pre-2000") and post-emergence ("recent," defined as "2000- present") and calculated frequency of Bd occurrence with 95% CI and Pr(no Bd) (as described above) by country and region. To assess whether data across Africa support our proposed timing of emergence statistically, we tested whether the frequency of samples with Bdinfection was dependent on time period for countries and regions that had both pre-emergence and post-emergence records (Chi-square tests using the chisq.test function in R; R Core Team, 2022). For cases where we rejected the null (p<0.05), we inferred the direction of change based on the values for frequency of Bd occurrence in historical and recent time periods. Statistical analyses were conducted in R (version 4.0.4) in RStudio (version 1.4.1106) using 'pwr,' 'binom,' and 'Hmisc' packages (R Core Team, 2022;Harrell, 2012;Doraj-Raj, 2014;RStudio Team, 2016;Champely, 2018).

Bd genotyping assay
To determine which lineages of Bd were present, we genotyped Bd from a subset of our swab samples taken from 32 frogs collected in Cameroon in 2013 (N=25) and Burundi in 2011 (N=7) using a custom genotyping assay (Table S4) (Byrne et al., 2017). Briefly, this genotyping assay uses the Fluidigm Access Array platform to perform microfluidic multiplex PCR on 192 target loci of which 188 are on the Bd nuclear genome, 3 are on the Bd mitochondrial genome, and one is designed to target the internal transcribed spacer (ITS) region of the related pathogenic fungus Batrachochytrium salamandrivorans. Each target locus is 150-200 base pairs (bp) long. First, extracted DNA from swab samples was cleaned and concentrated from 30µl to 10µl using an isopropanol precipitation protocol. All samples were preamplified in two separate PCR reactions, each containing 96 primer pairs at a final concentration of 500nM. For each preamplification PCR reaction, we used the FastStart High Fidelity PCR System (Roche) at the following concentrations: 1x FastStart High Fidelity Reaction Buffer with MgCl 2 , 4.5mM MgCl 2 , 5% DMSO, 200µM PCR Grade Nucleotide Mix, 0.1 U/µl FastStart High Fidelity Enzyme Blend. We added 1µl of DNA to each preamplification reaction and used the following thermocycling profile: 95˚C for 10min, 15 cycles of 95˚C for 15sec and 60˚C for 4min. Preamplified products were treated with 4µl ExoSAP-it (Affymetrix Inc.) and incubated for 15min at 30˚C, then 15min at 80˚C. Treated products were diluted 1:5 in PCR-grade water. The diluted products from each of the two preamplification reactions were combined in equal proportions and used for downstream amplification and sequencing.
Each preamplified sample was loaded into the Fluidigm Access Array IFC (Fluidigm, Inc.) for amplification. All samples were barcoded and pooled for sequencing on ¼ of an Illumina MiSeq lane using the 300bp paired-end kits at the University of Idaho IBEST Genomics Resources Core. We pre-processed all sequencing data as described in Byrne et al. (2017) and generated consensus sequences for all variants present for each sample at each locus. We filtered reads by selecting sequence variants that were present in at least five reads and represented at least 5% of the total number of reads for that sample/locus using dbcamplicons (https://github.com/ msettles/dbcAmplicons).
We constructed a phylogeny to explore the relationship of the Bd sequences collected for this study (N=32) to those from previously-published Bd isolates (N=31; Rosenblum et al., 2013;O'Hanlon et al., 2018). The isolates from Rosenblum et al. (2013) (N=28) were sequenced and processed exactly as described for the swab samples. The consensus amplicon sequences for the isolates from O'Hanlon et al. (2018) (N=3) were calculated from raw sequencing reads. Briefly, we downloaded raw reads from NCBI's Sequence Read Archive (accession numbers SRS2757071, SRS2757072, SRS2757205) and cleaned the reads using seqyclean v 1.9.9 (https://github.com/ibest/seqyclean). We aligned the paired reads to the reference genome JEL423 (Broad Institute v. 17-Jan-2007) using BWA MEM (Li, 2013). Aligned reads were sorted using picard v.2.9.0 (http://broadinstitute.github.io/picard) and variants were called using samtools mpileup v.1.6 and bcftools call v.1.6 (Li, 2011). We produced a consensus genome for each isolate using the bcftools consensus. After producing a consensus genome for each of the three isolates, we used BLAST v.2.7.1 to match our reference amplicon sequences to the amplicon sequences in each genome. We then extracted the best-scoring BLAST match for each isolate to produce a single FASTA containing a consensus sequence for each target locus.
The isolates represent three of the major Bd clades: Bd-GPL (Global Panzootic Lineage), Bd-BRAZIL/Bd-ASIA-2, and Bd-CAPE (Farrer et al., 2011;O'Hanlon et al., 2018). First, we trimmed our consensus sequence dataset to eliminate loci that had more than 50% missing data. The remaining dataset contained 186 loci (183 nuclear, 3 mitochondrial). To ensure proper alignment, we filled in missing sequences with a string of Ns (length = maximum sequence length of that locus) using custom R scripts. Next, we individually aligned all loci using the MUSCLE package v.3.18.0 (Edgar, 2004) Team, 2022) and concatenated all aligned loci for each sample for a total length of 25,617bp. We visually checked the concatenated alignment for errors in Geneious v.10.2.3 (Kearse et al., 2012) and used the RAxML plugin (Stamatakis, 2014) in Geneious to search for the best scoring maximum likelihood tree using rapid bootstrapping (1000 replicates).

Bd risk model
We integrated abiotic and biotic factors to predict the areas in Africa with the greatest risk of Bd establishment, following previous studies (Murray et al., 2011;Liu et al., 2013;Yap et al., 2018). We first created a presence-only habitat suitability model (HSM) driven by climate and land use factors using Maxent v.3.4.1 (Phillips et al., 2017). We used 1,914 georeferenced Bd occurrence records for continental Africa and the islands in the Gulf of Guinea (Bioko, São Tome, Prıńcipe); north of the Sahara, recent data were only available from Morocco. Records included all Bd-positives with latitude and longitude from our data and from published data. To reduce sampling bias, we applied a 10 arc-minute spatial filter by randomly choosing one Bd occurrence site from every 10 arc-minute (~12 km 2 ) area (Boria et al., 2014) using R packages 'dismo' (Hijmans et al., 2013) and 'maptools' (Bivand and Lewin-Koh, 2014). This resulted in 202 Bd-positive sites; 195 sites had environmental data and were used for model training. We further minimized sampling bias by restricting background sampling areas to a minimum convex hull polygon around large clusters of occurrence points (Phillips et al., 2009;Rödder et al., 2009;Thorne et al., 2012;Merow et al., 2013;Phillips and Elith, 2013;Mainali et al., 2015;Yap et al., 2018). Similar to Yap et al. (2018), we used 19 bioclimatic variables from the Worldclim database (http://www.worldclim.org/bioclim) (Hijmans et al., 2005) and the global human footprint (HF) (Venter et al., 2016) as predictor variables. To reduce the chances of overfitting the model due to multicollinearity, we calculated Spearman's rank correlations (r) among all the variables to determine which variables were highly correlated with each other (James et al., 2015). We removed one variable from each pair that had an r 2 > 0.7. We used the resulting subset of 11 environmental variables to create the Bd HSM: mean diurnal temperature range (Bio2), isothermality (Bio3), maximum temperature of the warmest month (Bio5), minimum temperature of the coldest month (Bio6), mean temperature of driest quarter (Bio9), precipitation of the wettest month (Bio13), precipitation of the driest month (Bio14), precipitation seasonality (Bio15), precipitation of the warmest quarter (Bio18), precipitation of the coldest quarter (Bio19), and human footprint (HF). These variables are biologically relevant to Bd and performed well in previous Bd HSMs (Ron, 2005;Rödder et al., 2009;Murray et al., 2011;Liu et al., 2013;James et al., 2015;Yap et al., 2018). We ran 20 replicates using cross-validation: occurrence data were divided into 20 equal-sized folds, or groups, and for each replicate 19 folds were used for model training and one fold was used for model testing, and the results were then averaged by Maxent to produce a final occurrence probability (Phillips et al., 2017).
In addition to abiotic factors, biotic factors such as host availability are critical for disease maintenance and spread (Lloyd-Smith et al., 2005). Although host abundance data would be ideal, such data at the community level for a continental assessment are not available for Africa. However, several studies have shown species richness can be an important factor in Bd spread (Rödder et al., 2009;Liu et al., 2013;Boria et al., 2014;Yap et al., 2016;Hydeman et al., 2017;Yap et al., 2018). Thus, in lieu of abundance data, we incorporated species richness to refine our predictions of Bd risk. We obtained amphibian richness data from AmphibiaWeb (amphibiaweb.org; AmphibiaWeb, 2021), which included range maps from the IUCN Red List as well as expert-based range maps supplied by AmphibiaWeb to calculate richness.
We calculated the product of the Bd HSM and amphibian richness to determine Bd risk (i.e., areas of overlap between suitable Bd habitat and potential hosts) (Yap et al., 2015;Yap et al., 2016;Yap et al., 2018). We show Bd risk on a scale of low to high using the Jenks optimization method, which identifies natural clusters, or classes, of data and reduces variance within the classes while maximizing the differences between classes (Jenks, 1967). This method is more informative than the binary threshold options provided by Maxent when predicting species range expansions (Creley et al., 2019).

Bd emergence
Our Bd-assay data for 4,623 amphibian specimens (collected 1908-2013) compiled with 12,234 Bd records from published data (collected 1852-2017) over a 165-year period show a pattern of rapid Bd emergence beginning largely at the turn of the century ( Figure 1A).
From 1852-1999, we found low overall Bd prevalence (3.2%, 2.7-3.7% CI) and limited geographic spread, but after the year 2000 (from 2000-2017), we found a sharp increase in continent-wide prevalence (18.7%, 18.1-19.4% CI) and wider geographic spread. When we examined Bd prevalence by decade, we found that it remained below 5% for every decade until the 2000s when it increased to 17.2% (16.6-18.0% CI; Figure 1B; Table S1). In the second decade of the new century (2010s), overall frequency continued to increase to 21.6% (20.4-22.9% CI). In several countries where we had more complete data, the pattern was more pronounced. For example, in Cameroon, we found only one Bd positive in samples we collected up to 1999, but in the next decade, the 2000s, the frequency of Bd occurrence rose to 11.0% (8.9-13.3% CI), followed in the 2010s, by an increase to 36.2% (33.0-39.5% CI; Figure 1C; Table S2). In Cameroon, the average Bd infection intensity also increased over time and in the most recent samples, maximum values (20002010present, 6.4 logZE) surpassed a threshold infection intensity associated with chytridiomycosis-caused mortality (4.0 logZE; Vredenburg's "10,000 zoospore rule"; Vredenburg et al., 2010;Kinney et al., 2011). In Kenya, frequency of Bd occurrence was zero or very low in the first four decades of data (the highest frequency in these early decades was in the 1970's at 5.2%, 2.5-9.3% CI; Figure 1D; Table S3) and increased approximately six-fold in the 2000s (31.5%, 28.4-34.7% CI; Figure 1D; Table S3). The average Bd infection intensity in Kenya also increased during this same time period, with maximum Bd-infection intensity surpassing the mortality threshold in the 2000s by several orders of magnitude (7.9 logZE; Figure 1D; Table S3).
In our test of whether frequency of samples with Bd-infection was dependent on time period (historical vs recent), we found that the overall frequency of Bd occurrence across Africa in the historical time period (pre-2000; 3.2%, 2.7-3.8% CI), was significantly different from the recent time period (2000-2017; 18.7%, 18.1-19.4% CI; Chi-sq., p<0.001; Table 1). We also found that the frequency of Bd occurrence was significantly different regionally (see regions; Figure 1A) across central, eastern, and southern Africa, between historical (<7%) and recent (>17%) time periods (Chi-sq., p<0.001, p<0.001, p<0.001, respectively; Table 1). One country in central Africa, Equatorial Guinea, showed a significant decrease in the frequency of Bd occurrence between historical and recent time periods (all samples from Bioko Island, Equatorial Guinea; Chi-sq., p=0.02), though all the Bd-positives from the historical time-period were collected close to the year 2000 (in 1998). Burundi, in eastern Africa, had the highest frequency of Bd infections after 2000 (73.7%, 66.3-80.2% CI) with no Bd positives detected historically; however, a Chi-square p-value could not be accurately estimated for this country because of the small historical sample size (N=13). We did not detect a recent surge of Bd infections in the country of South Africa (25.7%, 17.8-34.9% CI, historic vs. 23.3%, 21.1-25.6% CI, recent; Chi-sq, p=0.65), though we did in the region overall (6.6%, 5.0-8.5% CI, historic vs. 24.2%, 22.0-26.4% CI, recent; Chi-sq, p<0.001). There were no historical data for northern Africa, which only included data from Morocco, and the recent frequency of Bd occurrence there was low (2.0%, 0.9-4.0% CI). In western Africa, although the frequency of Bd occurrence also showed an increase from 0% historically to >20% recently, low historical sample size (N=10) prevented accurate estimation of a Chi-square p-value. Bd was almost entirely absent from western Africa (Table 1) except in Nigeria, where the recent frequency of Bd occurrence was 44.5% (40.8-48.2% CI).

Bd lineage genotyping
We found multiple lineages of Bd in Burundi and Cameroon (Figures 2A, B; Table S4; Byrne et al., 2019). In Cameroon, all genotyped samples formed a well-supported clade that includes the previously-published Bd-CAPE isolate CCB1 (Rosenblum et al., 2013), indicating that amphibians in Cameroon are widely infected by the Bd-CAPE lineage ( Figure 2B). In contrast, Bd samples collected in Burundi were more variable, with some samples grouping with previously-published Bd-GPL isolates (CAS250667, CAS250764, CAS250826, CAS250829) while others were not clearly members of either the Bd-GPL or Bd-CAPE lineages (CAS250662, CAS250742, CAS250744) ( Figure 2B). For one of these ambiguous samples, CAS250742, the average number of alleles found at each locus across all amplified loci was 2, which was higher than for all other sequenced samples ( Figure S1). Bootstrap support in this group of samples is low, thus no firm conclusions can be made from data for these three samples. However, the ambiguous Burundi samples are distinct from both from Bd-CAPE samples as well as other previously published Bd-CAPE/Bd-GPL hybrids ( Figure 2B; O'Hanlon et al., 2018).

Bd risk model
The Bd risk model predicted moderate to high Bd vulnerability throughout much of western, central, and eastern Africa, including within the Ethiopian Plateau (Figure 3). Some portions of Mediterranean northern Africa and the southeastern tip of Africa were also predicted to have some Bd risk. Our model predicted the highest Bd risk in the tropical rainforests of western and central Africa and the montane forests surrounding the Great Rift Valley, while less vulnerable areas included the arid Sahara Desert throughout most of northern Africa and the semi-arid Kalahari Desert in the western region of southern Africa.

Discussion
The 165 years of amphibian sample collections show a distinct signal of Bd emergence and geographic spread across the African continent that began in the first decade of the 21st century. In general, in the 20 th century (1852-1999), we found few, geographically western North America, studies found few and scattered Bd-infected specimens historically in areas that decades later experienced rapid Bd expansion, epizootics and host population collapse Cheng et al., 2011;Sette et al., 2015;Yap et al., 2016;Adams et al., 2017;De Leoń et al., 2017;Yap et al., 2018;Vredenburg et al., 2019;Sette et al., 2020). In Africa, we saw only one major exception to this general pattern of emergence: in South Africa, infections were widespread even in our oldest samples (Table 1). Thus, generally we found strong evidence indicating a recent surge (beginning in 2000) in the frequency of Bd occurrence in amphibian specimens collected across a large portion of the continent (Figures 1A, B; Tables 1, S1). Disturbingly, we found this broad pattern of Bd emergence in areas of Africa where many amphibian species are either suspected or known to have collapsed (Channing et al., 2006;Gower et al., 2012;Gower et al., 2013;Hirschfeld et al., 2016;Weldon et al., 2020). Considering that Bd invasion was followed by Bd-epizootics resulting in extinctions of amphibian populations in western North America Vredenburg et al., 2010), Central America (Berger et al., 1998;Lips et al., 2006), South America (Catenazzi et al., 2011;Carvalho et al., 2017), Europe (Walker et al., 2010;Bates et al., 2018), and Australia (Laurance et al., 1996;Berger et al., 1998), amphibian declines in Africa that coincide with widespread expansions of Bd may also be linked. Recent Bd emergence was most pronounced in eastern and central Africa where many species are threatened amid reports of declines and die-offs. In fact, most reports place the declines and extinctions occurring near or after the year 2000 (e.g. in Tanzania, Ethiopia, and Cameroon; Channing et al., 2006;Gower et al., 2012;Gower et al., 2013;Hirschfeld et al., 2016;Weldon et al., 2020), which coincides with our Bd emergence data. The historical and recent Bd occurrence data from Cameroon and Kenya are the most complete. In these countries, we saw a large and significant increase in the number and proportion of Bd infected hosts after 2000 (Figures 1C, D; Tables  S2, S3). Studies have shown that high levels of infection (i.e. Bd infection intensity or Bd load) are associated with severe chytridiomycosis, death, and population collapse (e.g., mortality threshold, Vredenburg's "10,000 zoospore rule"; Vredenburg et al., 2010;Kinney et al., 2011); and, low Bd infection intensity is associated with host survival and coexistence with Bd suggestive of enzootic dynamics Knapp et al., 2016;Catenazzi et al., 2017). In both Cameroon and Kenya, we detected a recent increase in the average Bd-infection intensity of hosts. We note, however, that Bd infection intensities from historical Bd occurrences may be unreliable (see Methods). Regardless, the maximum values for infection intensity that we documented from live animals in the most recent decade for Cameroon and Kenya surpass the mortality threshold (Figures 1C, D;  Tables S2, S3). Thus, we conclude that our evidence of recent and rapid emergence of Bd and presence of highly infected animals signal that Bd epizootics may be occurring in these regions. By contrast, in South Africa, we found no evidence of a recent Bd emergence  Table S4 for sample information associated with the newly genotyped swabs.
( Table 1). There have been a few cases of observed mortality and population bottlenecks associated with Bd in some regions of South Africa, however no widespread declines or extinctions have been linked to Bd despite comprehensive population survey data (Hopkins & Channing, 2003;Tarrant et al., 2013;Griffiths et al., 2018;Doherty-Bone et al., 2020). This suggests that either Bd may be endemic to South Africa, that there may be Bd genetic lineages with lower virulence there, or that sampling remains incomplete. Global studies on Bd genomics suggest that Bd has a complex history with many lineages that appear to have varied effects on amphibian populations (Rosenblum et al., 2013;Fisher and Garner, 2020). Generally, lineages identified as endemic to a single part of the world are thought to have a long evolutionary history of cooccurrence with hosts, and thus have evolved lower virulence than invasive (or novel) lineages that cause epizootics and host die-offs. In some places, endemic and invasive lineages are hypothesized to be cooccurring and hybridizing to form new lineage variants (O'Hanlon et al., 2018), and the effects of these hybrid lineages on hosts is unknown. Bd-epizootics can also be triggered by enigmatic abiotic or biotic factors that have not been deeply explored Fisher et al., 2021). Previous work has reported multiple genetic lineages of Bd in Africa (O'Hanlon et al., 2018;Byrne et al., 2019), suggesting the potential for complex hybridization and epizootic dynamics.
In Cameroon, we found evidence of Bd-CAPE, a lineage generally considered to be endemic to Africa and less virulent than Bd-GPL; however, our emergence data and results of histology (see Supplementary Material) suggest this lineage may in fact be virulent in Africa. All genotyped samples from Cameroon clustered genetically with Bd-CAPE ( Figure 2B), which supports previous work (O'Hanlon et al., 2018). Histological examination of animals from Cameroon show evidence of chytridiomycosis (see Supplementary Material;  Table S5; Figure S2), and we found Bd prevalence in the area has increased and spread geographically. The reported amphibian declines from Cameroon (Hirschfeld et al., 2016) coincide with the recent increase in frequency of Bd occurrence that we report here. We concur with previous work that states that pathogen lineage alone cannot fully predict disease outcomes since lineages that appear less pathogenic in native hosts may cause negative impacts on new or naïve host populations, especially when abiotic factors promote infection and disease (O'Hanlon et al., 2018;Fisher et al., 2021). We suggest this may be the case with Bd-CAPE in parts of Africa. It is possible that Bd-CAPE was present in Cameroon historically and has only recently increased in frequency due to shifts in some environmental, human, or other unknown factor, but it is also possible that the recent increase in prevalence represents a new invasion and spread of Bd-CAPE in Cameroon, perhaps from elsewhere on the continent. Bd-CAPE has been widely detected in South Africa (O'Hanlon et al., 2018;Byrne et al., 2019) and may be endemic there. In South Africa, we found a relatively constant proportional occurrence of Bd in hosts through a long time-period, and while there have been cases of Bd-associated mortality there, these could be associated with different Bd lineages (Bd-GPL is also present) or result from cofactors that are not yet known. Previous work has proposed that Bd-CAPE spread from South Africa to Europe, where it has thus far been associated with host declines (Doddington et al., 2013;O'Hanlon et al., 2018). Thus, we propose that Bd-CAPE may have recently spread within continental Africa and may also be associated with Bd-epizootics in some locations where conditions are suitable. Studies have reported recent evidence of mass amphibian die-offs in Cameroon (Blackburn et al., 2010;Hirschfeld et al., 2016) and Lesotho (Minter et al., 2004), areas where more detailed studies are needed to discern whether Bd-CAPE is associated with Bd-epizootics.
Our genomic data also identified Bd-GPL (Global Panzootic Lineage) in samples from Burundi. The Bd-GPL lineage is considered the most virulent Bd lineage (Farrer et al., 2011;Rosenblum et al., 2013;O'Hanlon et al., 2018;Byrne et al., 2019). Consistent with this, histological examination of specimens from Burundi showed signs of chytridiomycosis (see Supplementary Material; Table S5). In fact, our samples from Burundi had the highest proportion of Bd infected hosts in our recent time frame (73.7%; 66.3-80.2% CI), and this could suggest Bd-epizootics are occurring there. Furthermore, we found no historical Bd positive data in Burundi (Table 1) or in adjacent Tanzania (samples collected 1926-2003). In 2003, the Tanzanian Kihansi spray toad (Nectophrynoides asperginis) went extinct in the wild purportedly due to a Bd epizootic (Weldon et al., 2020). While there was no statistical support that frequency of Bd occurrence increased in the recent timeframe in Burundi and Tanzania (due to low availability of historical samples), we did find that frequency of occurrence was significantly higher statistically after 2000 in eastern Africa overall (Table 1). Therefore, the recent presence of Bd-GPL detected after a dramatic increase in Bd occurrence in eastern Africa along with histological evidence of chytridiomycosis, may indicate that Bd-GPL driven epizootics are occurring there. Predicted Bd disease risk in Africa from our habitat suitability model (HSM); Based on 1,914 georeferenced Bd occurrence records with 195 unique Bd-positive localities in the African mainland and islands in the Gulf of Guinea. The Bd risk model is a product of a Bd HSM model (similar to Yap et al. (2018); here using 10 bioclimatic variables and the global human footprint as predictor variables) and host richness maps (AmphibiaWeb) for the entire African continent. The results are shown on a color scale from low risk (tan color) to high risk (red color).
Hybridization between distinct Bd lineages has been proposed as a way that virulent lineages have emerged, and these hybrid lineages may be more virulent than either parent lineage (Farrer et al., 2011;Greenspan et al., 2018). Three of our genotyped samples from Burundi did not clearly cluster with previously characterized Bd lineages on the phylogenetic tree, but rather were positioned between Bd-GPL and Bd-CAPE ( Figure 2B). We see three possibilities for the identity of Bd lineages from these samples: (i) Bd-GPL and Bd-CAPE are able to hybridize, and while our ambiguous samples did not cluster with known hybrids, they could be novel hybrids of these two lineages ( Figure 2B; O'Hanlon et al., 2018); (ii) these samples could each contain multiple lineages of Bd, as evidenced by a higher average number of alleles found in one of these samples than any other ( Figure S1), and this supports the possibility that more than one lineage was present in that sample; or (iii) Bd detected in these samples could represent novel genetic variation that has not previously been characterized. Further work on ambiguous lineages in Burundi could resolve this issue and help in our general understanding of this complex pathogen. Further studies may identify lethal and sublethal infections (potentially from different Bd lineages) that could help us understand factors that could affect pathogen success (and how it may vary across lineages) such as differences in pathogen dispersal ability, reproductive ability, and vulnerability to predators. Additionally, we emphasize that future studies should prioritize genotyping Bd throughout Africa. The genotyping work presented here is limited, and further Bd genotyping is needed to better characterize the distribution of lineages across the continent.
Habitat suitability models combined with potential host availability provide an opportunity to understand the potential distribution of Bd risk over large geographic areas (Murray et al., 2011;Liu et al., 2013;Yap et al., 2018). The results of our HSM generally align with previous potential distribution models (Penner et al., 2013;Zimkus et al., 2020); however, two key distinctions in our model are that (1) we incorporated potential host availability to determine risk of disease spread (Yap et al., 2015;Yap et al., 2016), and (2) we included the global human footprint (Venter et al., 2016). Although host abundance data at the community and continental scale would be ideal to identify host availability, such data are not available. Therefore, as a proxy, we used species richness, which has been shown to be an important factor in Bd spread (Rödder et al., 2009;Liu et al., 2013;Boria et al., 2014;Yap et al., 2016;Hydeman et al., 2017;Yap et al., 2018). We included the global human footprint (Venter et al., 2016) in our model because the presence of humans or human activities can play a role in disease spread (Liu et al., 2013;James et al., 2015). We also took additional steps to reduce sampling bias and minimize the chances of overfitting due to multicollinearity of the abiotic variables. Thus, our model is complementary to previous potential distribution models (Penner et al., 2013;Zimkus et al., 2020). Our model identifies geographic areas with moderate to high risk of Bd spread as well as further highlighting African regions with high Bd suitability and supporting the assertion that countries in these regions are potential hotspots for the negative impacts of Bd (Figure 3).
In central Africa, the HSM predicted high Bd risk in Angola, Cameroon, Democratic Republic of Congo, Equatorial Guinea, Gabon, and Republic of Congo. The observed Bd dynamics in Cameroon and Democratic Republic of Congo, coupled with population declines in Cameroon (Hirschfeld et al., 2016) suggest a recent emergence, yet more information and monitoring is also needed in Angola, Equatorial Guinea, Gabon, and Republic of Congo to determine their Bd status (Table 1). We found a high proportion of Bd positives from Bioko Island (Equatorial Guinea) from 1998, and Bd-GPL has also been detected there (O'Hanlon et al., 2018;Byrne et al., 2019;Zimkus et al., 2020), but the lack of population survey data across a longer timescale prevents a proper understanding of Bd-host dynamics there.
In eastern Africa, the HSM model predicted high Bd risk in Kenya, Malawi, Mozambique, Tanzania, Zambia, and Zimbabwe. Mozambique and Zimbabwe border South Africa, where Bd-CAPE may be endemic and there is evidence suggesting spread to Lesotho could be associated with amphibian declines (Minter et al., 2004). Tanzania experienced a documented Bd-driven species extinction in the wild in 2003 (Weldon et al., 2020). Additionally, Tanzania borders Kenya, where we found evidence of Bd emergence and high Bd risk, and Burundi, where we report high Bd prevalence and ambiguous Bd lineages (Seimon et al., 2015). Because eastern Africa is identified as having widespread high Bd risk and multiple locations have high recent Bd prevalence, Bd-host dynamics in these countries should be closely monitored and studied.
In western Africa, our HSM predicted high Bd risk in Ghana, Cote d'Ivoire, Liberia, Guinea, and Sierra Leone. However, there are no records of Bd positives there, and it has been hypothesized that the Dahomey Gap (a non-forested area between western Nigeria and eastern Ghana) may act as a natural barrier creating a safe haven from chytridiomycosis for amphibians (Penner et al., 2013). However, it is worrisome that in nearby Nigeria (a country where we identified moderate Bd risk) the frequency of Bd occurrence is very high (44.5%; 40.1-48.2% CI), and our model predicted highly suitable habitat just west of the Dahomey Gap, similar to previous studies (Penner et al., 2013;Zimkus et al., 2020). Also, other studies have shown that once Bd invades an area, it spreads over large distances and into remote areas rapidly (Lips et al., 2006;Berger et al., 2016;Yap et al., 2016;Vredenburg et al., 2019;Fisher and Garner, 2020). We suggest that future research should focus on Bd-host relationships in western Africa, including monitoring disease presence and spread to determine whether species are threatened by future Bd invasion and surveillance of amphibians in Nigeria to determine whether population declines have occurred.
In northern Africa, our HSM predicted moderate risk in Morocco and in northern Algeria, Tunisia, and Egypt. However, this risk may be underestimated as our analysis did not include recent studies showing an expansion of Bd positives in Morocco (El Cadi et al., 2019;Thumsováet al., 2022). Despite these newer records for Morocco, data for northern Africa are limited; in fact, we are not aware of any published records for other countries in the region. Thus, northern Africa represents another understudied region where monitoring of amphibian populations and Bd are crucial.
New studies of Bd-host dynamics and interactions between Bd lineages in African amphibians could greatly advance our knowledge of chytridiomycosis. By bringing together new and published data spanning 165 years, our study demonstrates recent increases in the frequency of Bd occurrence in much of Africa; yet, we also find areas where frequency of Bd occurrence has not changed over time.
In some areas, increases in frequency of Bd occurrence are coincident with amphibian declines (e.g., Cameroon Tanzania, and Ethiopia; Hirschfeld et al., 2016;Channing et al., 2006;Weldon et al., 2020;Gower et al., 2012;Gower et al., 2013), whereas in regions where Bd has not changed in frequency over long periods of time, evidence of widespread declines has not been detected (e.g., South Africa; Hopkins and Channing, 2003;Tarrant et al., 2013). Combined, these results indicate that Africa may harbor regions of Bd endemism (where host species may not be susceptible to epizootics) as well as regions where Bd is invading (posing large risks to host species). We find that multiple genetic lineages of Bd occur in Africa and suggest that (1) lineages generally considered to be low in virulence may actually be responsible for epizootic dynamics in some areas (e.g. Bd-CAPE), and that (2) co-occurring lineages may be hybridizing. Our results provide an opportunity to further investigate virulence of Bd lineages in nature and to determine how Bd hybridization may affect virulence, pathogenicity, and health outcomes in hosts. Finally, our initial data and model predictions help expand our study to enable risk prediction in areas where no Bd data are presently available. Our study highlights the importance of determining the temporal context of host-pathogen dynamics and emphasizes the need to bridge historical and contemporary datasets to better describe host-pathogen dynamics over larger temporal scales and predict future scenarios.

Ethics statement
The animal study was reviewed and approved by San Francisco State University Institutional Animal Care and Use Committee (#A20-05-SFSU).

Author contributions
VV and DB conceived and designed the study. DB, EG, MTK, and DP collected field samples. DB and SG collected museum samples. SG and AB performed laboratory work. DM performed histology. SG, TY, HS, AC-A, SC, KL, AM, EP, MSK, and HR compiled and organized data from the literature. SG, TY, AB, ER, VV, and HS contributed to the analysis of the data. SG, VV, TY, AB, HS, DM, and AZ wrote the manuscript. All authors contributed to the article and approved the submitted version.