Enzootic frog pathogen Batrachochytrium dendrobatidis in Asian tropics reveals high ITS haplotype diversity and low prevalence

Emerging Infectious Diseases (EIDs) are a major threat to wildlife and a key player in the declining amphibian populations worldwide. One such EID is chytridiomycosis caused by Batrachochytrium dendrobatidis (Bd), a fungal pathogen. Aetiology of Bd infection is poorly known from tropical frogs in Asian biodiversity hotspots. Surveys were carried out in four biodiversity hotspots to ascertain the status of Bd fungus. We collected a total of 1870 swab samples from frogs representing 32 genera and 111 species. Nested PCRs revealed low prevalence (8.4%) and high Bd haplotype richness was revealed after sequencing. We document 57 Bd Internal Transcribed Spacer region (ITS) haplotypes, of which 46 were unique to the global database. Bd ITS region showed indels at the Taqman binding site and qPCR reverse primer binding site, suggesting qPCR is unsuitable for diagnosis in Asian Bd coldspots. Our median-joining network and Bayesian tree analyses reveal that the Asian haplotypes, with the exception of Korea, formed a separate clade along with pandemic BdGPL (Bd Global Panzootic Lineage) haplotype. We hypothesise that the frog populations in Asian tropics might harbour several endemic strains of Bd, and the high levels of diversity and uniqueness of Bd haplotypes in the region, probably resulted from historical host-pathogen co-evolution.

collected in 1902. New hypotheses are in contention after Bd was detected in a historical collection of frogs (Acris blanchardi and Rana sphenocephala) from 1888 in the USA 25 , and (Hypsiboas pulchellus) in Brazil from 1894 19 . A recent global study conducted using the next generation sequencing and phylogenetic analysis of 234 isolates revealed the Bd origin to be Asia 26 . Tip dating analysis predicts BdGPL might have originated 120-50 years ago 26 . Clearly, Bd infections are new to frog populations globally, and the factors that lead to their susceptibility to infections are still poorly understood.
The status of chytridiomycosis in Asia is poorly understood, with the exception of a large scale survey involving 15 countries in Asia, that revealed low Bd prevalence, and predicted the widespread occurrence of Bd in several parts of Asia 9 . The first report of a Bd outbreak came from Japan in 2006 27 . Following this, intensive surveys revealed several haplotypes of Bd in Japan 8 , China 28 , and Korea 18 . Such areas are referred to as coldspots of Bd infection and are typically characterised by a high number of haplotypes and low prevalence 15 . From the global Bd phylogeny, it is evident that, Korea harbors two enzootic Bd lineages namely BdASIA-1 and BdASIA-2/BdBRZIL 26 .
As reports on Bd from Asia are growing, it has widened the scope of investigations on the role of the pathogenic strain in the region, and the measures to be taken to secure several narrowly endemic frog populations from lethal infections. India lies at the intersection of four globally recognised megabiodiversity hotspots 29 with high levels of frog species richness and endemism. Bd is expected to be widespread in these hotspots 9,[30][31][32] . India has over 375 frog species, and the list is growing (see, http://amphibiaweb.org/). With such high frog species richness, the stakes are high for understanding the role of Bd on frog populations in the region. Apart from a few scattered efforts in the Western Ghats of India 10,[33][34][35] , the status of Bd in the hotspots is unknown. Since chytridiomycosis is a notifiable aquatic disease, and country-wide reports are far from complete, there is a compelling case for countries like India to report the status of Bd. With these considerations, we measured the prevalence of Bd in frog populations from four biodiversity hotspots in India, using nested Polymerase Chain Reaction (PCR). We found that Bd is prevalent at a low level in all the hotspots in India. We retrieved several Bd haplotypes by targeting the ribosomal Internal Transcribed Spacer region (ITS) and several of them were unique and therefore, we predict that there are enzootic lineages yet to be discovered from the region. Our haplotype network revealed that enzootic haplotypes clustered around BdGPL. The haplotype network further lends support to the hypothesis that Bd might have originated from Asia. We investigated the efficiency of conventional assays used for detection of Bd, and emphasise the need for efficient assays for detecting Bd in Asian coldspots. These findings have important implications for the ongoing global effort to understand impacts of the disease.

Results
Bd detection and prevalence. A total of 1870 swab samples (147 locations) from frogs that belonged to 32 genera and 111 species were analysed. Of these, 158 samples showed amplification in the nested PCR method for Bd (Table 1). We initially analysed 1050 samples using both qPCR and Nested PCR methods. Both methods varied in their sensitivity; qPCR and nested PCR showed 33 and 119 positive samples, respectively. Both qPCR and nested PCR methods had only 19 positives in common. Since nested PCR showed greater sensitivity in detecting Bd, it was used as our primary assay for measuring prevalence.
Prevalence was the highest in Micrixalidae with 17.4% (95% CI 7-37%) ( Table 2); all frogs in this family are endemic to the Western Ghats. Overall, the prevalence of chytridiomycosis in frog populations that were sampled in different regions of India was 8.44% (95% CI 7-97%). Seasonal variation in Bd prevalence coincided with the onset of the monsoons. Bd prevalence gradually increased from June, peaked in July (23.4%), and gradually dropped in September. From this point the prevalence progressed to 14.2% in October and subsequently decreased from November onwards (Table 3).
There was no significant difference in Bd prevalence among frog populations in various geographic regions (F value = 0.51, df = 5, p = 0.766). Rarefaction revealed that the Western Himalaya region had the maximum number of haplotypes of Bd (Fig. 1). Large variations in both prevalence and haplotypes richness were observed within the different geographic regions (Fig. 1).

ITS haplotype diversity and network. From 158
Bd positive samples, 57 haplotypes were retrieved.
We found single haplotypes in all samples, except for 31 samples, wherein multiple haplotypes were detected. From the positive samples, 49 sequences consisting of 34 haplotypes were resolved through cloning. Among them, 46 haplotypes were unique to India, and 11 haplotypes matched with those reported from China, Japan, South Africa, USA, and Italy. Out of the 57 haplotypes, 33 were exclusively from mainland India, and 19 were from Andaman and Nicobar Islands. Only five haplotypes were shared between the mainland and the Islands. Haplotypes varied in their amplicon size ranging from 247 to 263 bp. Two haplotypes (IN02 and IN10) contributed to 65% of the Bd positive samples. Haplotype IN02 was the most common one, and it was represented in 51% of the Bd positive samples. Haplotype IN10 was present in 20 samples, from all regions that had Bd, except for the Western Himalayas and the Nicobar Islands (Fig. 2). We estimate that there were 147 haplotypes (SE ± 0.44) involved in our samples after accounting for the undetected ones.
Four haplotypes, IN05, IN17, IN31, and IN55 had mutations at the Taqman probe binding site. They contributed to 3.1% (N = 158) of positives in our sample. Two out of the four haplotypes showed an addition of two adenine residues at the Taqman binding site and the other haplotypes showed transition and transversion of a single nucleotide base ( Supplementary Fig. S1).
In the Bayesian tree, haplotypes from Korea, Japan, and Brazil formed a distinct clade whereasthe haplotypes from India grouped with haplotypes from China, Japan, Italy, South Africa, and Texas. Both the clades were supported by posterior probability values. As expected, the Taqman binding mutation haplotypes, IN05, IN55, China CN 30, Japan JP02, JP09, and JP10 formed a separate clade ( Supplementary Fig. S2). Median-joining Network has three major clusters, the first one, representing a majority of the haplotypes from India, along with China, Japan, Italy, South Africa, and Texas. The second cluster consists of haplotypes mainly from Korea, Japan, and Brazil. The third cluster consists of haplotypes with mutations at the Taqman binding site from India, China, and Japan (Fig. 3). Some haplotypes show no specific association, such as those from Japan, as they are represented in all three clusters. The most common haplotype (IN02) found in India clustered with haplotypes from Italy, South Africa, Texas, China, and Japan (cluster I in Fig. 3). Bd JEL 423 haplotype clustered with the haplotypes, IN14 and CN13, from India and China, respectively.

Discussion
Using samples collected over a period of six years, representing approximately 25% of India's total frog species, we document for the first time that Bd is found in frog populations across India, in the Western Ghats, Eastern Ghats, Himalayas, North East Hills, and Andaman and Nicobar Islands. These geographic regions were predicted as areas that might harbour Bd, albeit, with low prevalence 9 . We demonstrate that nested PCR is a sensitive assay to quantify Bd prevalence in large samples from Bd coldspots that characteristically show low prevalence and a high number of enzootic haplotypes. Nested PCR is less expensive than qPCR 36 , and it has the ability to detect Bd haplotypes that might go undetected in qPCR. It provides nested PCR with an unassailable advantage for long-term surveillance in Bd coldspots. There is a low prevalence of the Bd infection in the Western Ghats, as documented previously 10,34,35 , but this was a common pattern, and it is repeated in other regions, with no significant differences between them. We attribute the overall low prevalence to detection bias, and patchy occurrence ( Table 1). The factors leading to patchiness in Bd occurrence in relatively homogeneous habitats needs further investigation. Nested PCR fared better than qPCR in retrieving Bd from samples, but several rare haplotypes in the samples represented by very few zoospores would have certainly remained undetected. Therefore, the Bd prevalence measured by this method would be an underestimation. The magnitude of the bias is unknown and therefore, correction factors elude us at this point. We emphasise the need for sensitive and universal assays for detecting all known Bd haplotypes. This would be an important pre-requisite for understanding the dynamics of the disease in Bd coldspots. The family with the highest number of Bd positive samples, Micrixalidae, is known for species that inhabit fast moving streams (torrent frogs), which is not surprising as the Bd zoospores are known to be aquatic. Seasonal variation in abiotic factors such as air temperature and rainfall are known to influence Bd prevalence [37][38][39][40] . The regions sampled receive two monsoons, viz: the southwest, and the northeast. A sharp rise in prevalence was seen in July and October, and these coincide with the monsoons. Such patterns in annual variation in Bd prevalence are known from other parts of the world 37,38,41 .
The number of haplotypes of Bd involved in infecting a frog population seems to fluctuate with the prevalence in that population (Fig. 1). This pattern could emerge when a pool of Bd haplotypes with varying abundance  The occurrence of coldspots of Bd, within global amphibian biodiversity hotspots, such as the case presented here, should be understood in the context of theoretical expectations of host-pathogen evolution 15 . Long-term monitoring of Bd prevalence in these areas would be vital in answering these questions. ITS region is the fungal DNA-barcode region used for species level identification 42 . Using the nested PCR assay, Bai et al. 28 reported 30 new haplotypes from mainland China, with two widespread haplotypes, CN18 and CN02. The haplotype CN18 was present in 78% of the total number of positive samples from this study. Similarly, we identified a single dominant haplotype, IN02, in 51% of the Bd positive samples in India. IN02 was also identical to haplotypes CN18, JP01, SA08 from China, Japan, and South Africa, respectively. We observed similarity in the reports from China and our study, wherein, low Bd prevalence was also associated with multiple haplotypes 8,28 . Based on this evidence, coldspots should be investigated for the role of dominant haplotypes involved in infections.
We now know that Bd occurs in-situ with several enzootic haplotypes in India. In this study, we could sequence only a little over a third of the Bd haplotypes that were probably present in the samples. A single strain of Bd could have multiple Bd ITS regions, resulting in many haplotypes 17 ; several endemic strains in Asia could be involved in causing infections.
Global trade in amphibians is held responsible for the spread of Bd across the world [43][44][45][46] . Historically, India was a major exporter of three large aquatic frogs: Hoplobatrachus tigerinus, Euphlyctis hexadactylus, and Hoplobatrachus crassus until early 1980's 47,48 . Trade of frogs in India ceased after a ban imposed on the collection and export of frogs 49 . There is no known record of the import of frogs into India. Migratory birds, crayfish, and movement of researchers are some of the possible sources 50-52 that might have seeded Bd fungus into India. The preponderance of enzootic haplotypes, and their close relationship to BdGPL, prompts us to hypothesise that frog populations in India might harbour several endemic strains of Bd, and the high levels of diversity and uniqueness of Bd haplotypes in the region, is probably resulted due to historical host-pathogen co-evolution.
Worldwide, Taqman based quantitative PCR is used as the primary diagnostic tool for detecting Bd 53 . The efficacy of this assay was questioned because mutations were identified at the Taqman binding site 8 . We found 30 haplotypes in the database that have mutations at the Taqman binding site. Three haplotypes from China, 13 from Japan, 9 from Korea, and 5 from Brazil show mutations at the Taqman binding site either as indels or as transition-transversions (Supplementary Table S1 and Fig. S1). We additionally report 15 haplotypes with mutations at the reverse primer binding site (Supplementary Table S2 and Fig. S3). The Taqman probe Chytr MGB was designed using 29 Bd sequences from Australia, the USA, Panama, and Ecuador 54 . Twenty six out of the 29 sequences were identical, and the remaining three sequences showed indels 54 . Presently, we have more than 250 ITS1 sequence accessions in NCBI GenBank, obtained either from simple PCRs or Nested PCRs 8,18,28,35,55 . If Bd coldspots are indeed associated with high Bd haplotype richness, and Taqman based qPCR is used as the primary assay, the probability of false negatives would increase, resulting in the failure to report Bd, when it occurs with low prevalence. Multiple sequence alignment of Bd ITS1 revealed that there were no conserved sites, making it an unsuitable target for new Taqman probes. Therefore, we speculate that Chytr MGB Taqman probe might not serve as an effective probe for Bd detection in Asian coldspots.
Bd fungus is present in all the frog hotspots in India and it needs to be matched with a heightened level of surveillance in the region. Our data highlight the need for sensitive and specific assays to detect them as conventional assays show limitations in detecting different haplotypes. Globally, Asia has the highest number of Bd haplotypes and the least number of Bd related mortalities. This presents an opportunity for investigating the strains involved in infections in the region, and the persistence of enzootic strains in Asian tropical frog populations.

Methods
Sample collection. Surveys were carried out in 147 locations representing seven distinct biogeographic regions in India 56 between November 2012 and March 2017 (Fig. 4). For the purpose of the surveys, all anurans were sampled and they are referred to, henceforth as frogs. During the surveys we selected a stream or a pond, and all frogs were captured ad libitum from 18:00 h to 24:00 h. Surveys targeted the two monsoons: southwest monsoon (between June and September), and northeast monsoon (between October and November). Sites that exceeded 600 m above sea level were selected in the Indian mainland, as the sites had optimal abiotic conditions for survival and growth of Bd. Each frog was hand caught using a pair of fresh non-powdered gloves and held in afresh plastic bag to eliminate cross-contamination. We captured at least 10 different frogs, without any preference for each species at any given site. Swab samples were collected from frogs by swabbing the skin using sterile cotton swabs (HIMEDIA ® PW003). Each swab was stroked 70 times on the frog as follows:10 strokes each onthe dorsal surface, lateral sides, ventral surface and undersides of the thighs and five outward strokes on the undersides of each foot 57 . After swabbing each individual was examined for clinical symptoms such as loss of righting reflex, abnormal body posture, skin sloughing, skin lesions, and other abnormalities, after which, the frog was released at the site of capture. Geo-coordinates and altitude were measured using the handheld GPS, Garmin eTrex 10. Used gloves and plastic bags were incinerated at the end of the sampling session. Instruments and boots were cleaned with bleach and wiped with ethanol before proceeding to the next sampling location. All methods were carried out in accordance with Florida International University's Institutional Animal Care and use Committee Guidelines (IACUC: 12 8 ] along with three tungsten carbide beads (3 mm diameter) and were homogenised in Qiagen Tissue Lyser II for 2 min followed by centrifugation at 8,000 rpm for 2 min. The supernatant containing DNA was transferred into fresh 1.5 ml tubes and incubated initially at 50 °C for 2 h, followed by incubation at 95°C for 20 min. After incubation, these tubes were centrifuged at 13,000 rpm for 10 min and the supernatant was transferred into fresh 1.5 ml tubes and stored at −20 °C. This DNA extract was diluted with nuclease-free water in the proportions of 1:10 to be used as a template for PCR analysis. Bd detection and haplotyping. Initially, we tested the efficiency of two assays for detecting Bd from our samples. For this, we used both Quantitative PCR (qPCR) and Nested PCR on roughly 50% of our samples. qPCR was carried out using the protocol from Kriger et al. 58 with a reaction volume of 10 µl. containing 1X Takara master mix (premix Ex Taq), 900 nM of forward, and reverse primer (ITS1-3 Chytr, 5.8S Chytr), 250 nM Taqman MGB probe, 1X ROX (TAKARA Bio.Inc.), 2 µl of diluted DNA, and the remaining volume made up by nuclease-free water. Cycle conditions for qPCR were as follows: initial denaturation at 95 °C for 10 min, followed by 50 cycles of denaturation at 95 °C for 15 sec; annealing, and extension at 60 °C, for 1 min. We used qPCR standards at 100, 10, 1 and 0.1 Genome Equivalents (GE). All reactions were carried out in duplicates with positive (Bd JEL 423 DNA) and negative control (nuclease-free water). For nested PCRs, the ITS1-5.8S-ITS2 region (~300 bp) was used as the marker. We followed Goka et al. 8 nested PCR protocol and performed the first PCR in a total 10 µl reaction containing 1X of Emerald Takara master mix, 1 picomole each of primers Bd18SF, and Bd28SR1 8 , 1 µl of 1:10 diluted template DNA and the remaining volume was made up by adding molecular grade water. The PCR conditions were as follows: initial denaturation at 95 °C for 10 min, followed by 40 cycles of denaturation at 95 °C for 30 sec, primer annealing at 50 °C for 30 sec, extension at 72 °C for 30 sec, and a final extension at 72 °C for 10 min. We included a positive control (Bd JEL 423 DNA), and a negative control (nuclease-free water) in every run. For the second PCR, we used Bd 1a and Bd 2b primers 8,59 , and 1 µl of first PCR product as DNA template for second PCR. PCR conditions as follows: initial denaturation at 95 °C for 10 min, followed by 30 cycles of denaturation at 95 °C for 30 sec, primer annealing at 65 °C for 30 sec, extension at 72 °C for 30 sec, and a final extension at 72 °C for 10 min. PCR products were visualised on 2% agarose gels, 300 bp bands were observed, and samples positive for Bd were re-amplified. Sequencing was done twice using Applied Biosystems (model 3730).
Few samples could not be reliably sequenced due to multiple Bd ITS DNA strands in a single sample, so we created clones in Escherichia coli to address this problem. For this, primers Bda (5′-CAGTGTGCCA TATGTCACG-3′) and Bdb (5′-CATGGTTCATATCTGTCCAG-3′) were modified by inserting the restriction sites for enzymes EcoRI and PstI at 5′-end of the primers. The resultant primers BdaP (5′-CACTGCAGT GTGCCATATGTCACG-3′) and BdbE (5′-GGAATTCCCATGGTTCATATCTGTCCAG-3′) were used for amplification of the 300 bp fragment. The PCR conditions used for amplification were: initial denaturation at 95 °C for 5 min, 30 cycles of amplification with each cycle containing 95 °C for 15 sec, 65 °C for 30 sec, 72 °C for 30 sec, and a final extension at 72 °C for 5 min. The amplicons were gel purified, digested with the enzymes EcoRI, and PstI, and after digestion they were re-purified by gel elution. The purified amplicons were cloned by ligating into a pBS-SK vector that was linearized with EcoRI and PstI. Ligated products were transformed into E. coli DH5α bacterial cells that were made competent. White transformants were selected on LB agar plates containing ampicillin (60 µg/ml) and X-gal (5-bromo-4-chloro-3-indolyl-β-D-galactoside; 20 µg/ml). Colony PCR was performed with the transformants using the M13 forward, and reverse primers. The amplicons obtained were then sequenced with the M13 primer.

Data analyses.
For estimating prevalence, we used 'binconf ' 60-62 function with Wilson score interval in the Hmisc package in R 63 (version 3.4.1) to estimate a 95% confidence interval (CI). We used Estimate S (version 9.1) to: (i) obtain second-order Jack-knife extrapolation estimate 64 to arrive at the number Bd haplotypes that could be present in the samples, as there is expected to be an incomplete detection of rare haplotypes in the samples 65 ; (ii) to check for differences in haplotype richness in different geographic regions by using rarefaction, for the smallest sample (N = 50) in any particular region. Samples from the North-East hills were excluded from this analysis, as there were no positives. One-way analysis of variance (ANOVA) was performed to test for differences in Bd prevalence in the geographic regions, using R (version 3.4.1). A Bayesian tree was prepared using Bd ITS1-5.8S-ITS2 sequences available at genebank, and the sequences generated in this study. The species, Terramyces subangulosm, Boothiomyces sp., Boothiomyces macroporosum, and Kappamyces laurelensis were used as outgroups (Supplementary Table S3). We used, Clustal W 66 in MEGA 67 (version 6.06) to align all the sequences and prepared a nexus file. JModel test was performed to obtainthe best evolutionary model using JModel 68 (version 2.1.10) using Bayesian Information Criteria (BIC). The HKY + G evolutionary model was identified as the best model. We used FastGap 69 (version 1.2) to create a gap matrix. A Bayesian tree was constructed using the new matrix obtained in MrBayes 70,71 (version 3.2.6) with HKY + G evolutionary model. Fifteen million generations were chosen, and sampling was done after every 500 generations. To ensure effective sampling size, Tracer (version 1.6) was used. After the initial 25% of trees were discarded as burn-ins, the resultant tree was visualised in FigTree (version 1.4.2). All Bd ITS sequences, except the outgroup, were aligned using ClustalW 66 in MEGA 67 (version6.06) with some manual editing. A FASTA file was generated through MEGA 67 (version 6.06), and it was used in DnaSP 72 (version 5) to count the number of haplotypes and prepare Roehl data file. Network(version 5.0; http://www.fluxus-engineering.com) was used to build haplotype network using median-joining method 73 .