The Fungicide Chlorothalonil Changes the Amphibian Skin Microbiome: A Potential Factor Disrupting a Host Disease-Protective Trait

: The skin microbiome is an important part of amphibian immune defenses and protects against pathogens such as the chytrid fungus Batrachochytrium dendrobatidis (Bd), which causes the skin disease chytridiomycosis. Alteration of the microbiome by anthropogenic factors, like pesticides, can impact this protective trait, disrupting its functionality. Chlorothalonil is a widely used fungicide that has been recognized as having an impact on amphibians, but so far, no studies have investigated its effects on amphibian microbial communities. In the present study, we used the amphibian Lithobates vibicarius from the montane forest of Costa Rica, which now appears to persist despite ongoing Bd-exposure, as an experimental model organism. We used 16S rRNA amplicon sequencing to investigate the effect of chlorothalonil on tadpoles’ skin microbiome. We found that exposure to chlorothalonil changes bacterial community composition, with more signiﬁcant changes at a higher concentration. We also found that a larger number of bacteria were reduced on tadpoles’ skin when exposed to the higher concentration of chlorothalonil. We detected four presumed Bd-inhibitory bacteria being suppressed on tadpoles exposed to the fungicide. Our results suggest that exposure to a widely used fungicide could be impacting host-associated bacterial communities, potentially disrupting an amphibian protective trait against pathogens.


Introduction
Amphibians around the world are increasingly threatened by diseases caused by fungi, viruses, bacteria, and parasites [1,2]. Particularly, the infectious skin disease, chytridiomycosis, is one of the main diseases impacting amphibian health [3,4]. This disease caused by the chytrid fungal pathogens Batrachochytrium dendrobatidis (Bd) and B. salamandrivorans (Bsal) can cause mass die-offs in amphibian species [5,6]. The skin microbiome is considered one of the first lines of defense against pathogenic infections and can mediate disease susceptibility [7][8][9], suggesting it is an essential part of the amphibian's innate immune system. In amphibians, protection against pathogens has been linked to distinct characteristics of the skin bacterial communities, such as bacterial species richness, microbial community assemblage, and the presence and abundance of members in the bacterial communities capable of producing metabolites that suppress pathogen infections [10][11][12][13]. In addition, formerly common throughout the mountain ranges of Tilarán, Cordillera Central, and Talamanca in Costa Rica and western Panama [37]. This species suffered population declines and disappearances across its entire range in the late 1990s and was considered possibly extinct. Six years after disappearing, this species has been re-encountered at different sites in isolated populations in Costa Rica and is reproducing in high numbers [38]. The disease chytridiomycosis was one of the main drivers of declines for L. vibicarius, possibly in combination with habitat disturbance, pesticides, and climate change [38]. The skin bacterial community may provide protection against Bd in this species [39].
We hypothesized that exposure to chlorothalonil would change the skin bacterial communities of tadpoles. We predicted that tadpoles exposed to chlorothalonil would have different bacterial communities than those not exposed to chlorothalonil and that this change would be more emphasized at higher fungicide concentrations. Thus, we investigated changes in the relative abundance of members in the bacterial communities when exposed to varying fungicide concentrations. We identified the presence of putative Bd-inhibitory bacteria with differential abundance between treatments. The present study serves as a baseline to understand the effect of a toxic fungicide on the microbially mediated immune defenses of a montane tropical amphibian, as well as the potential risk faced by remaining populations of L. vibicarius in landscapes where chlorothalonil can potentially be present.

Tadpole Collection and Maintenance
In October 2017, we collected 80 tadpoles of L. vibicarius with similar developmental stage (Gosner stage [27][28][29], body size (mean ± standard deviation (SD): 41.5 ± 4.4 mm). and weight (0.71 ± 0.19 g) from a permanent lagoon in the Juan Castro Blanco National Park, Alajuela, Costa Rica. Due to the high abundance of tadpoles we observed in the lagoons of the Juan Castro Blanco National Park during our monitoring program, we concluded that the number of collected animals did not have an impact on this population. The study and ethical procedures were approved by National Commission for the Biodiversity Management of Costa Rica (R-057-2019-OT-CONAGEBIO) and the Ministry of Environment and Energy of Costa Rica-National System of Conservation Areas (SINAC-ACAHN-PI-R-010-2017).
We captured the animals with nets and placed them in sterile plastic trays containing pond water. Animals were transported to a laboratory at the University of Costa Rica to carry out the chlorothalonil exposure experiment. All tadpoles were placed in an aquarium with filtered water and acclimatized to laboratory conditions for 8 days before starting the experiment. We consider an 8-day acclimation period an adequate time to acclimatize both host and microbiome to experimental conditions previous to any manipulations [40,41]. We know that the amphibian skin microbiome can change under captive conditions [42]; therefore, all tadpoles were set up under the same conditions to maintain the same initial microbiome baseline between treatments and reduce any potential effect in the results.

Chlorothalonil Exposure and Sampling
Following the acclimation period, we conducted an 8-day exposure experiment to investigate the effect of chlorothalonil on the skin microbiome of tadpoles. We established four treatments, which consisted of a negative control (filtered water), a solvent control (SC; methanol), and two concentrations of chlorothalonil; low concentration (1 µg/L) and high concentration (5 µg/L) (nominal concentrations). In Costa Rica, chlorothalonil has been detected in the environment (e.g., soil, air, and water), and concentrations above 11 µg/L have been reported in water bodies [18,26,43]. In addition, the concentrations of chlorothalonil used in this experiment are similar to demonstrably nonlethal concentrations used in another study using a species of the same genus, Lithobates taylori (E. Ballestero, unpublished data). Therefore, the exposure levels were chosen to reflect conditions that the species may experience in the wild. We prepared chlorothalonil exposure solutions by adding an aliquot of a stock solution to the exposure medium (filtered water). Stock solution (1061 µg/mL) was prepared from 97.5% pure chlorothalonil standard (Dr. Ehrenstorfer, Germany) dissolved in HPLC-grade 99.97% methanol (J.T. Baker, Phillipsburg, Unites States) and kept at 4 • C. Aliquots were taken using a microvolume syringe (SGE Analytical science, Australia). The quantitative analyses of chlorothalonil were performed with solid-phase extraction (SPE) and gas chromatography-mass spectrometry (GC-MS) at the Instituto Regional de Estudios en Sustancias Tóxicas (IRET), Universidad Nacional, Costa Rica. The actual concentrations for the low-concentration and high-concentration treatments at the beginning of the experiment were 0.9 µg/L and 5.3 µg/L, respectively. We did not detect chlorothalonil in water samples at the end of the experiment (<0.1 µg/L).
In water, the half-life of chlorothalonil ranges from 0.18 to 8 days [44]. The amount of methanol added in the solvent control was the same as that one used in the highest concentration of chlorothalonil. We measured temperature ( • C), pH, and dissolved oxygen (mg/L) during the experiment (Supplemental Data Table S1).
Our experimental unit was one randomly chosen tadpole in a 1 L glass jar containing 800 mL of filtered water. We established 20 replicates per treatment (Table S2). Tadpoles were randomly assigned to one of the experimental treatments. We fed animals ad libitum with organic Spirulina on Days 0, 3, and 6.
On Day 8, we collected skin bacterial samples (skin swabs) from each animal. Swabbing consisted of moving a sterile rayon-tipped swab (Peel Pouch DryswabTM Fine Tip) across the animal skin. The swabbing protocol consisted of 12 strokes on each side (along body and tail), 12 strokes on the dorsal surface of the body, and 12 strokes on the mouth. We placed swabs in sterile vials with 300 µL of DNA/RNA Shield (Zymo Research). The tubes were transported to Ulm University, Germany, and stored at −20 • C until DNA extractions and sequencing. We used tricaine methanesulfonate (MS222) to euthanize all tadpoles at the end of the experiment.

DNA Extraction and 16S rRNA Gene Amplicon Sequencing
We extracted bacterial genomic DNA from swabs using the NucleoSpin Soil kit (Macherey-Nagel, Düren, Germany) following the manufacturer's protocol. We amplified the hypervariable V4 region of the 16S rRNA gene using the primers 515F (5 -GTGCCAGCMGCCGCGGTAA-3 ) and 806R (5 -GGACTACHVGGGTWTCTAAT-3 ). We followed the Fluidigm scheme (Access Array System for Illumina Sequencing Systems, Fluidigm Corporation), in which PCR and barcoding occur simultaneously. The PCR and barcoding (15 µL volume) were performed as described in Jiménez et al. [39]. Barcoded samples were purified using NucleoMag NGS Beads (Macherey-Nagel, Düren, Germany) and quantified with picogreen on Tecan F200. Then, we pooled all samples to an equal amount of 12 ng of DNA and diluted the pool down to 6 nM. Finally, the pooled sample library was paired-end sequenced in a single run on an Illumina MiSeq platform at the Institute of Evolutionary Ecology and Conservation Genomics, Ulm University, Germany. Raw sequence data were deposited into NCBI Repository, BioProject ID PRJNA703661.

Bioinformatics
The initial processing of the sequences was performed using QIIME 2 (version 2019.1) as described in Jiménez et al. [39]. For dada2 analysis, we trimmed the first bases of each read to remove primers (-p-trim-left-f 23, -p-trim-left-r 20) and truncated forward and reverse reads to 200 bp due to decreasing average quality scores of the sequences at the end. We collapsed reads into amplicon sequence variants (ASVs) and assigned bacterial taxonomy using the Greengenes database (version 13_8) as reference (http://greengenes.lbl.gov; accessed on 30 November 2020). We removed sequences classified as chloroplast, mitochondria, archaea, eukaryota, and unclassified phylum. We built a phylogenetic tree of the bacterial ASVs for further diversity analyses using MAFFT [45] and Fast Tree 2 [46]. Then, we imported our data into the R environment version 3.6.3 (https://www.r-project.org/; accessed on 30 November 2020) for further processing of the sequences using the R pack-age "phyloseq" [47]. We removed ASVs with less than 20 reads in the entire dataset and excluded samples with fewer than 9000 sequences. The resulting mean library size across individuals was 18,597 reads (range 9415-43,854). For alpha and beta diversity analyses, we rarefied the ASV table according to the sample with the lowest number of reads.

Statistical Analysis in R Environment
Because we used methanol as a carrier solvent, we first compared the negative control to the solvent control to detect statistical differences in alpha and beta diversity. We did not detect significant differences between negative and solvent controls considering the three alpha diversity measures (ASV richness: p = 0.94, Shannon diversity index: p = 0.97 and PD: p = 0.95; Figure S1). However, we observed that beta diversity differed between negative and solvent controls based on the ASV presence-absence composition (unweighted UniFrac: R 2 = 0.05; p = 0.040) and the ASV abundance-weighted composition (Bray-Curtis: R 2 = 0.12; p = 0.001) ( Figure S2). Thus, the solvent control was used as the basis of comparison in further analysis.
To investigate the effect of chlorothalonil treatments on the skin bacterial alpha diversity measures (ASV richness, Shannon diversity index, and phylogenetic diversity (PD)), we used Generalized Linear Models (GLMs) with Gaussian distribution. We log-transformed the alpha diversity measures prior to model fitting.
To examine the effect of chlorothalonil treatments on skin bacterial beta diversity, we calculated the unweighted UniFrac (based on ASV absence/presence) and Bray-Curtis dissimilarity (based on abundance pattern of ASV) metrics using the R package "phyloseq". We fitted permutational multivariate analyses of variance (PERMANOVAs) using the adonis function of the R package "vegan" [48] to statistically test the effect of fungicide treatments on both beta diversity metrics. We performed permutational pairwise post-hoc tests with a Bonferroni correction to evaluate statistical differences of bacterial beta diversity between chlorothalonil treatments and solvent control. Additionally, we quantified the extent of the difference between group centroids (treatments) as a measure of effect size by calculating Cohen's d and 95% CI with the R package "compute.es" [49]. We performed principal coordinate analysis (PCoA) to visualize the beta diversity distances between treatments.
Then, to identify ASVs that were significantly suppressed or overabundant in the two fungicide treatments, we used a negative binomial model-based approach (exact binomial test generalized for overdispersed counts) using the R package "edgeR" [50]. We present only ASVs that differed significantly between fungicide treatments and solvent control (FDR-corrected p-values at p < 0.001). We used the unrarefied dataset and the Trimmed Mean of M-values (TMM) method for the normalization of samples. Additionally, we explored the presence of putative Bd-inhibitory bacteria being suppressed or overabundant from the results of the previous analysis. To identify the putative Bd-inhibitory ASVs, we queried our ASV sequences against a database of culturable anti-Bd bacteria identified from different amphibian species (Antifungal Isolates Database; [51]). We retained ASVs with a ≥99% sequence identity match to those in the mentioned database following the methods outlined by Muletz-Wolz et al. [52] with the software Geneious version 20.1.2.

Chlorothalonil Disturbs the Skin Microbiome Beta Diversity
We did not detect a significant effect of chlorothalonil treatments on the three alpha diversity measures (ASV richness: p = 0.90, Shannon diversity index: p = 0.90, and PD: p = 0.33; Figure S1). The PERMANOVA models revealed a significant effect of treatments on the ASV presence-absence composition (unweighted UniFrac: R 2 = 0.09, p = 0.02, Figure 1a) and the ASV abundance-weighted composition (Bray-Curtis dissimilarity: R 2 = 0.09, p = 0.004, Figure 1b). Pairwise PERMANOVA tests indicated significant differences in the bacterial community composition between the tadpoles kept in the solvent control and the treatment with a high concentration of the fungicide, whereas those between solvent control and the treatment with a low concentration were similar (Table 1).
Pairwise PERMANOVA tests and Cohen's d effect sizes showed that the R 2 and effect sizes between tadpoles of the solvent control and the high-concentration treatment were higher than between tadpoles of the solvent control and the low-concentration treatment, indicating a greater difference in the higher concentration of fungicide (Table 1). Furthermore, skin bacterial communities of tadpoles exposed to a high concentration of fungicide clustered more distantly from solvent control treatment on the PCoA than those exposed to a low concentration of fungicide (Figure 1a,b). 0.33; Figure S1). The PERMANOVA models revealed a significant effect of treatments on the ASV presence-absence composition (unweighted UniFrac: R 2 = 0.09, p = 0.02, Figure  1a) and the ASV abundance-weighted composition (Bray-Curtis dissimilarity: R 2 = 0.09, p = 0.004, Figure 1b). Pairwise PERMANOVA tests indicated significant differences in the bacterial community composition between the tadpoles kept in the solvent control and the treatment with a high concentration of the fungicide, whereas those between solvent control and the treatment with a low concentration were similar (Table 1). Pairwise PER-MANOVA tests and Cohen's d effect sizes showed that the R 2 and effect sizes between tadpoles of the solvent control and the high-concentration treatment were higher than between tadpoles of the solvent control and the low-concentration treatment, indicating a greater difference in the higher concentration of fungicide (Table 1). Furthermore, skin bacterial communities of tadpoles exposed to a high concentration of fungicide clustered more distantly from solvent control treatment on the PCoA than those exposed to a low concentration of fungicide (Figure 1a,b).

Chlorothalonil Shifts Relative Abundance of Bacterial Strains
We detected that ASVs differed significantly between fungicide treatments and SC (Figure 2a,b). We found three ASVs suppressed when tadpoles were exposed to a low concentration of chlorothalonil, and 14 ASVs showed an increased abundance (Figure 2a). Tadpoles exposed to a high concentration of chlorothalonil showed 13 suppressed ASVs and seven overrepresented (Figure 2b). Tadpoles exposed to a high concentration had a higher number of bacterial taxa with reduced abundance than those exposed to a low concentration of fungicide. In both treatments, a low and high concentration of chlorothalonil, tadpoles showed a significant decrease in abundance for ASVs of the genera Sulfuricurvum and Janthinobacterium, whereas an abundance increment was observed for ASVs of the genera Nevskia, Flavobacterium, and Runella. We identified four putative Bd-inhibitory ASVs with significantly lower abundance in tadpoles exposed to chlorothalonil, one in tadpoles exposed to a low concentration and four in those exposed to a high concentration (Figure 2a,b). These putative Bd-inhibitory ASVs were assigned to the family Comamonadaceae and the genus Janthinobacterium, Acinetobacter, and Novosphingobium. and seven overrepresented (Figure 2b). Tadpoles exposed to a high concentration had a higher number of bacterial taxa with reduced abundance than those exposed to a low concentration of fungicide. In both treatments, a low and high concentration of chlorothalonil, tadpoles showed a significant decrease in abundance for ASVs of the genera Sulfuricurvum and Janthinobacterium, whereas an abundance increment was observed for ASVs of the genera Nevskia, Flavobacterium, and Runella. We identified four putative Bd-inhibitory ASVs with significantly lower abundance in tadpoles exposed to chlorothalonil, one in tadpoles exposed to a low concentration and four in those exposed to a high concentration (Figure 2a,b). These putative Bd-inhibitory ASVs were assigned to the family Comamonadaceae and the genus Janthinobacterium, Acinetobacter, and Novosphingobium.

Discussion
The present study highlights the impact of chlorothalonil, a widely used fungicide, on the immunologically important skin microbial community of a threatened frog that persists despite ongoing exposure to Bd. We provide evidence that exposure to chlorothalonil changes the skin bacterial community of tadpoles of L. vibicarius and, even more importantly, that chlorothalonil can suppress putative Bd-inhibitory bacterial strains at high concentrations. These results raise new concerns and hypotheses that need to be addressed to have a broader understanding of the impact of fungicides on the protective relationship between skin microbiomes and amphibian hosts.

Discussion
The present study highlights the impact of chlorothalonil, a widely used fungicide, on the immunologically important skin microbial community of a threatened frog that persists despite ongoing exposure to Bd. We provide evidence that exposure to chlorothalonil changes the skin bacterial community of tadpoles of L. vibicarius and, even more importantly, that chlorothalonil can suppress putative Bd-inhibitory bacterial strains at high concentrations. These results raise new concerns and hypotheses that need to be addressed to have a broader understanding of the impact of fungicides on the protective relationship between skin microbiomes and amphibian hosts.
Our results show that the skin bacterial community differed when tadpoles were exposed to higher concentrations of chlorothalonil. The effect sizes on beta diversity indicated that differences in bacterial communities of tadpoles increased as animals were exposed to higher concentrations of the fungicide. These findings indicate a change, and potential disruption, of skin bacterial communities as the concentration of chlorothalonil increases. Given that the microbiome can play a role in host health maintenance and that disruption of the natural range of microbial communities may lead to increased incidence of diseases in their hosts [15,53,54], our results suggest that exposure to chlorothalonil may increase susceptibility to diseases. Further work is needed to corroborate these hypotheses. A previous study found an alteration of the skin bacterial communities of tadpoles of Blanchard's cricket frog (Acris blanchardi, family Hylidae) when exposed to an herbicide [55]. Together, these studies suggest that changes in bacterial communities occur with a distinct type of pesticides and highlight an interest in investigating how pesticide mixtures (for example, a mixture of fungicides and herbicides) could be interacting to impact the function of microbial communities of amphibians.
Previous studies indicate that microbial communities are among the first taxa to respond to chemicals exposure [56][57][58]. Chlorothalonil could be directly and/or indirectly interacting with the skin bacteria of L. vibicarius tadpoles by altering their cutaneous bacterial communities. Microorganisms are functionally or nutritionally connected to each other, and changes in one component of a microbial community (e.g., the fungal community (mycobiome)) can influence the structure of the entire community [59][60][61]. Therefore, the fungicide chlorothalonil could be changing the skin mycobiome and thereby influencing the observed changes in the bacterial communities. In this study, we did not evaluate the skin mycobiome, so further research investigating interactions between fungal and bacterial communities after chlorothalonil exposure will provide broader insight into the impact of this fungicide on host-associated microbial communities. In addition, chlorothalonil could be affecting tadpoles' physiology, endocrine, and immune systems, as previously observed in different amphibian species [34], indirectly altering the host bacterial communities through different mechanisms. Here, we have not evaluated the host physiology mechanisms that could be changing these communities. However, we suspect that chlorothalonil could be contributing to an increment of stress on exposed individuals, as observed in tadpoles of the Cuban tree frog Osteopilus septentrionalis (family Hylidae) [34], thus affecting the tadpoles' skin bacterial community and resulting in a potential detriment of the host immune defenses. It is also possible that chlorothalonil may have interfered with the host skin peptide secretions that act as a selective force that controls which microbes can grow on each host's skin. The properties of the skin peptides have been shown to be affected by exposure to the insecticide carbaryl in yellow-footed frogs (Rana boylii, family Ranidae) [62]. These potential alterations to the skin properties by chlorothalonil may have disrupted the appropriate conditions for some bacteria to grow, consequently suppressing their abundance and altering host microbial communities. Further research is needed to investigate these potential mechanisms and will provide a better understanding of the impact of chlorothalonil on an amphibian immune defense trait.
We found different patterns in the relative abundances of bacterial taxa across chlorotha lonil treatments. We also observed the suppression of some putative Bd-inhibitory bacterial strains (for example, strains of genus Janthinobacterium and Acinetobacter) when exposed to chlorothalonil, particularly in the high-concentration treatment. These protective bacteria are known to be capable of producing metabolites that suppress Bd infections [10,11,63]. Together, this suggests that changes in bacterial abundances from chlorothalonil exposure could be disrupting the adequate production of defensive bacterial metabolites that facilitate disease resistance. Further, the suppression of some Bd-inhibitory bacteria provides evidence that chlorothalonil can interfere with these protective taxa, highlighting a potential risk in the disruption of host susceptibility to chytrid infections in early and/or later life stages. This knowledge is relevant as a reduction in bacterial abundances could represent the loss of key bacterial species and functions linked to host health. Previous evidence suggests that amphibian larval stages that have been exposed to chlorothalonil have higher Bd intensity and greater Bd-induced mortality when challenged with Bd after metamorphosis [36]. Therefore, exposure to chlorothalonil in early life might alter the normal bacterial community that establishes a healthy and Bd-protective skin microbiome after metamorphosis, making exposed animals vulnerable to future Bd infections. Based on this information, it would be interesting to investigate if this early-life disruption of bacterial abundances has lasting impacts on host Bd resistance later in life. It is also possible that opportunistic microbes and/or parasites with tolerance to chlorothalonil increase their abundance altering the natural microbial structure. Addressing this gap in our knowledge will allow a better understanding of the development of the immune system and will provide information that will help prevent early disruption of host microbiomes to confer better protection against diseases, such as chytridiomycosis.
We observed some individuals in the high-concentration treatment showing reddish skin, suggesting some type of external dermatitis probably attributed to the fungicide exposure. This is an interesting observation to highlight from our study because skin irritation has been linked to chlorothalonil exposure in animals [64]. This is merely an observation during our experiment but brings out the need for further investigation because chlorothalonil exposure has been associated with skin irritation and contact dermatitis in humans [65][66][67]. Atopic dermatitis can also allow for the colonization of certain types of bacteria that trigger immune response such as inflammation that can worsen symptoms and jeopardize host health [68].
Our understanding of how pesticides influence the amphibian microbiome is still in its infancy [16]. Our study reveals the effect of the exposure to environmentally relevant concentrations of the fungicide chlorothalonil on the skin microbiomes of amphibians in the early-life stage, which may, in turn, impact the stability of host-microbe interactions and microbiome-fitness correlations. Further studies are desperately needed, so we can fully understand the interaction between pesticides, the disease-causing organisms, and how these effects scale up to play a role in amphibian disease dynamics.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/10 .3390/applmicrobiol1010004/s1, Figure S1: Alpha diversity metrics of skin bacterial communities of tadpoles exposed to chlorothalonil, Figure S2: Principal coordinate analysis (PCoA) of skin bacterial communities of tadpoles kept in water and solvent controls, Table S1: Physicochemical parameters measured in water and solvent controls, and chlorothalonil treatments.

Informed Consent Statement: Not applicable.
Data Availability Statement: All raw sequence data were deposited into NCBI Repository, BioProject ID PRJNA703661.