- Shaanxi Key Laboratory of Qinling Ecological Security, Shaanxi Institute of Zoology, Xi’an, China
Microorganisms on amphibian skin reduce disease susceptibility and play an important role in pathogen defense. We hypothesized that anuran skin bacterial communities would change in response to seasonal variation and host species. To test this hypothesis, we used 16S rRNA amplicon sequencing to identify cutaneous bacterial communities of two frogs from the Qinling Mountains of China, Pelophylax nigromaculatus and Nanorana quadranus. We matched the amplicon sequence variants (ASVs) of microbes exhibiting protective effects against the pathogen Batrachochytrium dendrobatidis (Bd), using a database containing over 1900 16S rRNA gene sequences from amphibian skin bacteria. The results showed that seasonal variation had a stronger effect than host species on the structure (alpha-diversity, beta-diversity, species composition and abundance, and biomarkers) and anti-Bd function of cutaneous bacterial communities. These communities were highly dynamic but varied similarly between hosts. Their structural similarities were more consistent at the phylum level, but markedly less so at finer taxonomic levels. The highest relative abundance of anti-Bd reads was observed in P. nigromaculatus during summer, but anti-Bd reads were present in both frog species during different seasons. Therefore, the protective function of cutaneous microbial communities appears to be continuous despite between-species differences in anti-Bd ASV abundance. This observation does not directly explain why Bd infections have not been recorded in the region, butprovides important insight on anuran pathogen defense mechanisms. Our findings also suggest that specific seasons can be periods of high infection risk, with major implications for research on amphibian ecology and conservation.
1 Introduction
Bacterial communities on the skin are an important first line of defense against pathogens in amphibians (Brunetti et al., 2021; Jani et al., 2021; Longo et al., 2015). Experimentally removing skin bacteria increases morbidity associated with chytridiomycosis, whereas the addition of certain nascent bacteria reduces morbidity and mortality in some amphibian species (Becker et al., 2009; Harris et al., 2009; Kueneman et al., 2016; Walke et al., 2015). Symbiotic bacterial communities have gained greater recognition in mediating protection against a wide range of pathogens by modulating and contributing to host immunity (Bletz et al., 2017b; Fraune et al., 2015; Khosravi and Mazmanian, 2013; Rosenberg et al., 2007; Woodhams et al., 2007). Moreover, skin bacterial protection may be linked to specific bacterial metabolites and volatile compounds (Becker et al., 2009; Woodhams et al., 2018). For instance, violacein, prodigiosin, and volatile organic compounds produced by amphibian skin bacteria inhibit two species of chytrid fungi, Batrachochytrium dendrobatidis (Bd) and B. salamandrivorans (Bsal) (Brucker et al., 2008a; Woodhams et al., 2018), the criminal ringleader of amphibian chytridiomycosis.
Chytridiomycosis is a fatal skin disease that has significantly contributed to global decline and extinction of amphibian populations (Scheele et al., 2019). Of the two pathogens, Bd exhibits a broader host range, including a large number of Anuran, Caudate, and Gymnophiona species. While Bsal typically infects salamanders, it can also affect anuran species, including Osteopilus septentrionalis (Basanta et al., 2022; Gray et al., 2023; Laking et al., 2017; Towe et al., 2021). Bd infection exhibits temporal variation, and Bsal infection may exhibit similar pattern (Bletz et al., 2017b; Longo et al., 2015; Kinney et al., 2011; Longo et al., 2010; Phillott et al., 2013). This temporal variation is largely driven by temperature fluctuations (Bletz et al., 2017b). Therefore, temporal or seasonal variations in skin microbiota may be an important factor in disease dynamics. Amphibian skin microbial communities may also vary in response to habitat and pathogen presence, as well as host-specific factors such as species and developmental stage (Bletz et al., 2017a; Longo and Zamudio, 2017; McKenzie et al., 2012; Prest et al., 2018; Rebollar et al., 2016; Rebollar et al., 2019).
The Qinling Mountains serve as the natural demarcation line between northern and southern China. These mountains are important biodiversity hotspots but have recently experienced a sharp decline in amphibian populations. The most evident causes include climate change and habitat destruction, but pathogen infection has not been ruled out. Interestingly, no Bd- or Bsal-infected amphibians have been reported to date. However, no studies have investigated the assemblages and structures of amphibian cutaneous bacterial community in the Qinling Mountains; hence, we cannot conclude whether the lack of infection is due to superior anti-pathogen protection on amphibian skin or an absence of pathogens in the region.
To address these knowledge gaps, our study focused on two questions. First, what happens to the composition and diversity of amphibian cutaneous bacteria communities in response to seasonal variation and host species differences? Second, do any resultant fluctuations in microbial communities influence their protective function against pathogens such as Bd? We obtained data from two common and abundant aquatic anuran species in the Qinling Mountains, Pelophylax nigromaculatus and Nanorana quadranus. We sampled the cutaneous bacterial community across seasons, calculated their diversity and structure, as well as assessed whether chytridiomycosis was present and amplicon sequence variants (ASVs) from cutaneous microbes exhibited anti-Bd properties. Our findings should facilitate the elucidation of the cutaneous microbial ecology of frogs from the Qinling Mountains. Additionally, understanding the biological mechanism associated with pathogen susceptibility, even in species that are not currently in decline, can help us to establish how this range of susceptibility relates to the cutaneous skin microbiota of amphibians.
2 Materials and methods
2.1 Sample collection
Adult individuals from two aquatic anurans [the running-water frog, N. quadranus (n = 20), and the quiet-water frog, P. nigromaculatus (n = 12)] were sampled in spring, summer, and autumn of 2023 at Huangbaiyuan, located in the southern foot of the Qinling Mountains. Frogs were not sampled in winter because they hibernate during that season. The two species were selected for their abundance, increasing the likelihood of obtaining adequate sample sizes across seasons. For details on specific sampling times and locations (see Table 1).
Frogs were captured manually by researchers wearing sterile high-density polyethylene gloves (one pair per subject) and then rinsed thrice with purified water to flush away dirt and transient bacteria. Skin on the back, abdomen, and limbs was swabbed 30 times using a sterile skin swab. These swabs were placed in 2 mL of DNA storage solution (consisted of Tris, EDTA-2Na, and NaCl; Shanghai Langfu Industrial) and stored at 4°C. Collected samples were transported on dry ice to Beijing Biomarker Technologies for sequencing.
2.2 DNA extraction and sequencing
Total genomic DNA was extracted from all swabs using a TGuide S96 Magnetic Soil/Stool DNA Kit (Tiangen Biotech, Beijing, China), including lysozyme pretreatment. Extracted DNA was used as templates to amplify the V4 region of the 16S rRNA gene using barcoded primers pairs (515F: 5′-GTGYCAGCMGCCGCGGTAA-3′; 806R: 5′-GGACTACNVGGGTWTCTAAT-3′) in a polymerase chain reaction (PCR) reaction (Bletz et al., 2017b). The PCR products were quantified through agarose gel electrophoresis and purified using a DNA purification kit (Omega, Norcross, GA, USA). Purified amplicons were subjected to paired-end sequencing (2 × 250 bp) on the Illumina Novaseq 6,000 platform.
2.3 Sequence processing
Sequence processing and analyses were performed using BMK Cloud.1 Raw data were primarily filtered based on single-nucleotide quality in Trimmomatic (version 0.33) (Bolger et al., 2014). Primer sequences were identified and removed in Cutadapt (version 1.9.1) (Martin, 2011).
Further quality control of data was performed using filterAndTrim function, setting maxEE to 2 [EE = sum(10^(−Q/10))] and other parameters as default. Model construction was performed using learnErrors function, de-noising was performed with dada2 function, and double-end reads splicing was performed using mergePairs function (setting parameters: minOverlap: 18, maxMismatch: 18*0.2). Chimera removal was performed using the removeBimeraDenovo function (select the consensus method) (Edgar et al., 2011; Edgar, 2013; Callahan et al., 2016).
2.4 Sequence analysis
Feature classification was conducted using clean reads to generate ASVs through dada2, and ASVs with counts <2 across all samples were filtered (Callahan et al., 2016). These ASVs were matched against the SILVA database (release 138.1) in QIIME2 for taxonomic annotation, based on the Naive Bayes classifier with a confidence threshold of 70% (Quast et al., 2013).
Shannon, Simpson, Chao1 and ACE indices were calculated for frog skin samples using QIIME2 and displayed using ggplot2 (version 3.1.1). Between-group (species) differences in alpha diversity were determined with Wilcoxon test. Significance was set at p < 0.05. Beta diversity was calculated with unweighted UniFrac and Binary jaccard and visualized with principal coordinates analysis based on dissimilarity matrices (PCoA) (Bolyen et al., 2019). Between-group differences in beta diversity were tested using PERMANOVA [Adonis function in the R (version3.1.1) Vegan (version 2.3–0) package] (Dixon, 2003; Anderson, 2017). Unweighted pair group method with arithmetic mean (UPGMA) in Python was used to determine clustering patterns across samples. Linear discriminant analysis-effect size (LEfSe) was used to test for significant taxonomic differences between groups (Segata et al., 2011). A logarithmic LDA score of 4.5 was set as the threshold for discriminative features to analyse the effect of season on bacterial communities, while the LDA score of 4.0 was set to analyse the effect of host species on bacterial communities.
Next, to explore the protective effect of cutaneous bacterial communities, a database was utilized which containing over 1900 16S rRNA gene sequences from amphibian skin bacteria that have been tested for activity against the pathogen, Bd (Woodhams et al., 2015). Positive hits were then matched with ASVs present in the samples to calculate the proportion of anti-Bd reads (100%matching rate). Correlations between proportion of inhibitory reads and seasonal or species differences were tested using ANOSIM.2
In the sequence analysis, frog samples were divided into groups based on capture date. Therefore, P. nigromaculatus caught in spring, summer, and autumn were, respectively, labeled as “SpringPn” (SpringPn1–8), “SummerPn” (SummerPn1–6), and “AutumnPn” (AutumnPn1–6). Because no N. quadranus was caught in the spring, these frogs were divided into two groups: “SummerNq” (SummerNq1–6) and “AutumnNq” (AutumnNq1–6) (Table 2).
Table 2. The number of amplicon sequence variants (ASVs) and alpha diversity indices in each sample.
2.5 Detection of chytridiomycosis
The extracted DNA was then subjected to detection through PCR utilizing Bd and Bsal nested primer pairs. The outer nest primer pairs of Bd and Bsal are ITS1f1 (5′- CTT GGT CAT TTA GAG GAA GTAA −3′) and ITS4 (5′-TCC TCC GCT TAT TGA TAT GC-3′). The inner nest primer pairs of Bd are Bd1a (5’-CAGTGTGCCATATGTCACG-3′) and Bd2a (5’-CATGGTTCATATCTGTCCAG-3′), and the inner nest primer pairs of Bsal are STerF (5′TGCTCCATCTCCCCCTCTTCA3′) and STerR (5′TGAACGCACATTGCACTCTAC3′).
The PCR reaction system consisted of 2 × Taq PCR Mix II (catalog# KT211-02, Tiangen Biotech) in a volume of 10 μL, with 1 μL each of upstream and downstream primers (10 μM), 1 μL of template DNA, and 7 μL of deionized water (ddH2O), resulting in a total reaction volume of 20 μL. The PCR amplification conditions were 10 min at 95°C, followed by 30 cycles of 10 s at 95°C, 10 s at 53°C (outer nest primer pairs) /52°C (inner nest primer pairs of Bd) /58°C (inner nest primer pairs of Bsal), and 10 s at 72°C and a final 10 min at 72°C. The resulting target gene fragments of Bd and Bsal were of lengths 296 bp and 161 bp, repectively.
3 Results
3.1 Overview of cutaneous bacterial communities
After sequencing, data filtering, and sequence splicing of 16S rRNA amplicons from 33 frog skin samples, we obtained 1,054,720 sequences. These sequences were further processed and clustered into 1943 ASVs, predominantly from five phyla: Proteobacteria (558), Firmicutes (283), Actinobacteriota (199), Bacteroidetes (276), and Acidobacteriota (157).
Table 2 and Figure 1 show the number of ASVs per sample and relative ASV abundance by phyla per sample, respectively. We found 88 ASVs common to all samples, belonging to Proteobacteria (38), Firmicutes (17), Bacteroidota (15), Actinobacteriota (10), Acidobacteriota (2), Cyanobacteria (1), Acidobacteriota (1), Verrucomicrobiota (1), Desulfobacteriota (1), Fusobacteriota (1), unclassified_Bacteria (1) and Unassigned (1) (Supplementary Figure 1).
Figure 1. Histogram of skin microbial distribution at the phylum level. Different colors indicate different species; stacked columns are the top 10 taxa in relative abundance at each taxonomic level.
3.2 Seasonal variation influenced bacterial communities
The SpringPn, SummerPn, and AutumnPn groups shared 202 ASVs, or 10.66% of total ASVs found in these samples (Supplementary Figure 2). The SummerNq and AutumnNq groups shared 255 ASVs, accounting for 14.74% of the total across both datasets (Supplementary Figure 3).
The Chao1 and ACE indices showed that the microbial species richness of P. nigromaculatus significantly differed across the three seasons, decreasing in the order of AutumnPn > SpringPn > SummerPn (Table 2; Figures 2A,B). Shannon and Simpson indices showed that the AutumnPn group had significantly greater species diversity than the SpringPn and SummerPn groups; however, the latter two groups did not differ (p > 0.05, Table 2; Figures 2C,D). For N. quadranus, species richness and diversity were significantly greater in the AutumnNq group than in the SummerNq group (Figures 2E–H).
Figure 2. Box plot of variation in alpha diversity indices for cutaneous bacteria on Pelophylax nigromaculatus across three seasons (A–D) and on Nanorana quadranus across two seasons (E–H). The horizontal coordinates are the group names, and the vertical coordinates are the values of the corresponding alpha diversity indices. (A,E) Wilcoxon test of Chao1 index; (B,F) Wilcoxon test of ACE index; (C,G) Wilcoxon test of Shannon index; (D,H) Wilcoxon test of Simpson index.
Using PCoA to analyse the unweighted UniFrac and Binary jaccard distance indices, we found seasonal differences in the species composition of cutaneous bacterial communities. On the PC1 axis, SummerPn clustered to the left, AutumnPn to the right, and SpringPn in the center. A large gap in their bacterial communities existed between the three seasons. On the PC2 axis, SpringPn clustered on the upper side, while SummerPn and AutumnPn were separated by a large gap with SpringPn (Figures 3A,B). For N. quadranus, the PC1 axis was the main factor contributing to a large difference between SummerNq and AutumnNq (Figures 3A,B). PERMANOVA (p = 0.001, Treatments 1, 2 in Tables 3, 4) verified the seasonal difference in the bacterial communities between the different seasons for both frog species.
Figure 3. Principal coordinates analysis of beta diversity in cutaneous bacteria on all groups. (A) unweighted UniFrac distance; (B) Binary jaccard distance. Each point represents the skin bacterial community of an individual sample.
Table 3. Summary of PERMANOVA models (unweighted UniFrac distance) of beta diversity for microbial communities on frog skin.
Table 4. Summary of PERMANOVA models (Binary jaccard distance) of beta diversity for microbial communities on frog skin.
The top five phyla composition of skin microbes in the SpringPn, SummerPn and AutumnPn groups possessed high similarity, although their relative abundance varied (Figure 1; Table 5). In the SummerNq and AutumnNq groups, three of the top five phyla were also identical, but with different relative abundance (Figure 1; Table 5).
The UPGMA clustering tree (unweighted UniFrac distance) combined with the species distribution histogram (genus level) revealed that samples collected in the same season were more similar in species composition (Figure 4). Between-season differences in composition and relative abundance were significant (Figures 4A,B and Table 6). The top five genera observed in AutumnNq accounted for only approximately 20.89% of cutaneous microorganisms (Table 6). The relative abundance of Lysobacter (1.97% ± 0.20%) and unclassified_Microscillaceae (1.81% ± 0.12%) in AutumnPn group ranked 11th and 12th, respectively, and did not appear in SummerNq group.
Figure 4. Clustering tree from unweighted pair group method with arithmetic mean (unweighted UniFrac distance) combined with species distribution histogram (genus level) for P. nigromaculatus in three seasons (A) and N. quadranus in two seasons (B). Clustering tree on the left, species distribution histogram on the right.
Next, LEfSe (LDA score of 4.5) of bacterial populations with between-group differences in relative abundance revealed multiple biomarkers (Figure 5). At the genus level, the most abundant biomarker in SpringPn was Gardnerella, while that in summerPn was unclassified_Erwiniaceae (Figure 5A). The most abundant biomarkers in SummerNq was Endozoicomonas (Figure 5B).
Figure 5. Histogram of the distribution of LDA values (LDA score of 4.5), comparing P. nigromaculatus in three seasons (A) and N. quadranus in two seasons (B). The vertical axis represents the taxonomic units exhibiting significant differences between the groups, while the horizontal axis displays bar graphs illustrating the logarithmic scores of LDA analysis for each respective taxonomic unit.
3.3 Effect of host species on cutaneous bacterial communities
The Venn diagram of cutaneous bacterial communities from frog samples revealed that 468 ASVs were shared between SummerPn and SummerNq, accounting for 56.05% of all ASVs in these two groups (Supplementary Figure 4). Additionally, 1,227 ASVs (87.39% of all) were shared between the AutumnPn and AutumnNq groups (Supplementary Figure 5).
Alpha diversity did not differ (p > 0.05) between the two species during summer or autumn (Figures 6A–H). PCoA analysis demonstrated that the elliptical circles formed by each of the samples of two frog species in same season partially overlapped (Figures 3A,B). Therefore, P. nigromaculatus and N. quadranus appear to have similar bacterial communities in summer and autumn. The results from PERMANOVA (p > 0.05, Treatments 3, 4 in Tables 3, 4) confirmed that the beta diversity of the two frog species did not differ during the same season.
Figure 6. Box plot of variation in alpha diversity indices for cutaneous bacteria on P. nigromaculatus and N. quadranus in the same seasons, including summer (A–D) and autumn (E–H). The horizontal coordinates are the group names, and the vertical coordinates are the values of the corresponding alpha diversity indices. (A,E) Wilcoxon test of Chao1 index; (B,F) Wilcoxon test of ACE index; (C,G) Wilcoxon test of Shannon index; (D,H) Wilcoxon test of Simpson index.
Samples from the two frog species differed in microbial phyla and genera composition and abundance during the same season (Figures 1, 7). The UPGMA clustering tree (Figures 7A,B) indicated that skin samples from the two frog species share similar microbial species composition within the same season.
Figure 7. Clustering tree from unweighted pair group method with arithmetic mean (unweighted UniFrac distance) combined with species distribution histogram (genus level) for P. nigromaculatus and N. quadranus in the same seasons, including summer (A) and autumn (B).
When LDA score was set to 4.0, the SummerNq group had one biomarker, Romboutsia, and the SummerPn group had none (Figure 8). Neither frog species exhibited any biomarker in autumn.
Figure 8. Histogram of the distribution of LDA values (LDA score of 4.0), comparing P. nigromaculatus and N. quadranus in the same seasons, including summer. The vertical axis represents the taxonomic units exhibiting significant differences between the groups, while the horizontal axis displays bar graphs illustrating the logarithmic scores of LDA analysis for each respective taxonomic unit.
3.4 Anti-Bd activity was stable across seasonal variation and host species
Of the 1943 ASVs identified, 924 (47.56%) matched the anti-Bd ASV database. Within all samples, anti-Bd ASVs accounted for an average relative abundance of 53.33% and were predominantly from five phyla: Proteobacteria (301), Bacteroidota (164), Firmicutes (161), Actinobacteriota (100), Acidobacteriota (66). Some anti-Bd ASVs appeared only in one frog species, while others appeared in both. We detected 46 anti-Bd bacterial ASVs common to all the groups. These were predominantly Proteobacteria (19), Bacteroidetes (9), Firmicutes (2), Actinobacteriota (5), and Acidobacteriota (6).
The percentage of anti-Bd ASVs varied between groups. The SummerPn group had the highest percentage (327/616, 53.08% anti-Bd ASVs), followed by SummerNq (349/687, 50.80% anti-Bd ASVs), SpringPn (490/1000, 49.00% anti-Bd ASVs), AutumnNq (573/1298, 44.14% anti-Bd ASVs), and AutumnPn (577/1333, 43.29% anti-Bd ASVs). To summarize, in P. nigromaculatus, more anti-Bd ASVs were present in the summer than in spring or autumn, with the fewest present in autumn. For N. quadranus, more anti-Bd ASVs were present in summer than in autumn.
The variation in the relative abundance of anti-Bd reads (percentage of summed reads associated with anti-Bd ASVs, relative to all reads) differed from that in the percentage of anti-Bd ASVs. The relative abundance of anti-Bd reads in SpringPn, SummerPn, AutumnPn, SummerNq, and AutumnNq was 40.20, 72.56, 46.60, 68.47, and 43.21%, respectively. ANOSIM indicated a significant difference across seasons for P. nigromaculatus (Bray-Curtis, R = 0.843, p = 0.001) and N. quadranus (Bray-Curtis, R = 1, p = 0.002). However, anti-Bd reads did not differ across host species in same seasons.
In addition, the PCR results for chytridiomycosis showed all the samples tested negative for Bd and Bsal.
4 Discussion
Seasonal fluctuations and host-specific factors can influence pathogen resistance through their effects on the amphibian skin microbiome (Costa et al., 2016; Ellison et al., 2019a; Ellison et al., 2019b; Longo et al., 2015; Muletz-Wolz et al., 2017). Herein, we used samples from P. nigromaculatus and N. quadranus across different seasons to find empirical support for this hypothesis in anurans from the Qinling Mountains. Our results showed that seasonal variation had a significantly stronger effect on skin bacterial community structure than host species. Stronger effects were observed in changes to microbial alpha-diversity, beta-diversity, species composition and abundance, biomarkers, and anti-Bd function.
Notably, while the phylum level showed higher consistency, the compositions of bacterial communities were much less consistent at the finer taxonomic levels. For example, Lysobacter and unclassified_Microscillaceae appeared on the skin of both frog species in autumn but was not found in summer and was very rare in spring. Lysobacter species are known as “peptide production specialist” and produce peptides that disrupt the cell walls or cell membranes of other microorganisms (Panthee et al., 2016). This means that seasonal variation are associated with peptide production and bacteriostatic action of frog skin. Gardnerella was present only in the SpringPn group, despite a high relative abundance. Because we did not capture any N. quadranus in spring, we could not verify whether Gardnerella spp. would be present on N. quadranus in that season. Gardnerella is associated with the etiology of human bacterial vaginosis (Shvartsman et al., 2023), but its pathogenic effects on frog skin remain unclear. In addition, unclassified_Erwiniaceae appeared in two samples from the SummerPn group. Erwiniaceae can produce antibiotics and inhibit the growth of entomopathogenic fungi (Cambronero-Heinrichs et al., 2023). Therefore, the presence of pathogenic fungi in these two samples from the SummerPn group may have caused an increase in unclassified_Erwiniaceae. We ruled out the possibility of Bd or Bsal infection in these two samples using PCR.
Our observations on P. nigromaculatus and N. quadranus suggest that skin bacterial communities are highly dynamic environments. Seasonality in particular had a strong effect on cutaneous bacterial community structure. Bacterial species richness on P. nigromaculatus decreased from spring to summer and increased from summer to autumn, while bacterial richness on N. quadranus increased from summer to autumn. The consistent increase from summer to autumn may be attributable to new bacterial taxa colonizing the skin. Indeed, we observed that 1,076 ASVs on P. nigromaculatus and 1,043 ASVs on N. quadranus were present in autumn but not in summer (Supplementary Figures 2, 3). Notably, P. nigromaculatus and N. quadranus showed similar variation in community structure despite being different species and from different locations. In both frogs, microbial species richness and diversity were significantly greater in autumn than in summer. Additionally, Actinobacteria was more abundant in autumn than in summer, whereas Proteobacteria was more abundant in summer than in autumn.
The prominent seasonal shift in bacterial community structure in these two aquatic frog species is probably related to changes in water temperature at both locations. Temperature directly influences the growth of bacterial community assemblages, community-member interactions, and antifungal functions (Bletz et al., 2017b; Daskin et al., 2014; Woodhams et al., 2014). Although environmental microorganisms were not sampled and analysed in this study literature shows that environmental temperature can alter amphibian gut microbiota and microbial communities in surrounding habitats (Bletz et al., 2017b; Kohl and Yahn, 2016). For aquatic frogs, bacteria in paddy fields, ponds, rivers, and near-water soils are important microbial reservoirs that colonize the skin; bacterial communities in these habitats also exhibit considerable temporal variation (Crump and Hobbie, 2005; Rasche et al., 2011). Therefore, seasonal (temperature-related) shifts in bacterial communities may dictate host-associated microbiota variations. The absence of environmental controls lead us to the current lack of clarity on the changes in the bacterial community specific to frog skin.
Anti-Bd bacterial ASVs (47.56%) within our all samples were higher than those reported on anuran species from Panama (8.47%), but lower than those reported from India (51.7%) (Mutnale et al., 2021; Varela et al., 2018). The factors contributing to the varying percentages of anti-Bd ASVs included soil pH, precipitation, prevalence of Bd infection, and others (Mutnale et al., 2021; Varela et al., 2018). Notably, the average percentage of anti-Bd ASVs in summer (51.94%) were higher than those in spring and autumn in our study. Therefore, the sampling season had a significant effect on the percentage of anti-Bd bacterial ASVs.
The top five most abundant genera in SummerPn all possessed anti-Bd properties, representing 60.15% of anti-Bd ASV abundance in the group. Overall, P. nigromaculatus skin bacteria in summer had the highest anti-Bd ASV abundance (72.56% average relative abundance), but all the test groups exhibited potential anti-Bd activity. This functional stability resulted from the constant presence of specific ASVs and replacement by different ASVs in some seasons. For example, anti-Bd Pseudomonas ASV (ASV 376) and Acinetobacter ASV (ASV 1770) were consistently present across all the three seasons, whereas anti-Bd Serratia ASV (ASV 162) only appeared in spring for P. nigromaculatus and autumn for N. quadranus. Moreover, Janthinobacterium ASV (ASV 31908) only appeared in spring and summer for P. nigromaculatus. In vitro experiments have demonstrated that Pseudomonas and Acinetobacter, isolated from amphibians, produced metabolites with Bd inhibitory activity (Brucker et al., 2008b; Becker et al., 2015; Woodhams et al., 2015). Additionally, Janthinobacterium lividum and Serratia marcescens can produce anti-Bd (and generally antifungal/antibacterial) metabolites such as prodigiosin and violacein (Woodhams et al., 2018; Zhang et al., 2024). Therefore, our predicted anti-Bd genera align with prior finding from in vitro experiments. Overall, seasonal shifts in host-associated bacterial communities were associated with significant variation in potential anti-Bd functions. Consequently, hosts may temporarily lose important protective bacterial taxa during certain seasons, drastically increasing their susceptibility to this disease.
Anti-Bd bacteria can produce Bd inhibitory metabolites or biofilms on frog skin. The host itself may also possess mechanisms to select anti-Bd metabolite-producing bacteria on the skin (Woodhams et al., 2018; Mutnale et al., 2021; Loudon et al., 2014; Piovia-Scott et al., 2017; Loudon et al., 2016). In this study, we demonstrated that the two tested frog species exhibited abundant anti-Bd bacterial ASVs on the skin. However, it is unclear whether this is directly related to the temporary absence of Bd infections in anuran populations from the Qinling Mountains. Moreover, the anti-Bd ASV database assembled by Woodhams et al. (2015) does not include any sampling from Qinling Mountains of China; hence, it cannot be ruled out that there are certain bacteria in the Qinling Mountains that might have anti-Bd activity but have not been matched with the anti-Bd ASV database. This area may be the focus of our subsequent studies.
5 Conclusion
Our results indicated that seasonal variation exerted a greater effect than host species on the structure and anti-Bd function of cutaneous bacterial communities. These communities exhibited more consistent structural similarities at the phylum level but were more diverse at finer taxonomic levels. Although both frog species hosted bacteria with anti-Bd function, the skin of P. nigromaculatus during summer had the highest anti-Bd bacterial ASV abundance. Our study provides important insights into the seasonality of cutaneous bacterial community structure on amphibians inhabiting the Qinling Mountains. The findings establish foundational knowledge for future studies on host–bacteria interactions and the structure–function relationship within amphibian skin bacterial communities. The findings from this study and other related studies can potentially benefit research on disease etiology and control in other wildlife hosts of cutaneous microbes.
Data availability statement
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://www.ncbi.nlm.nih.gov/, NCBI BioProject, PRJNA1134887.
Ethics statement
Ethical approval was not required for the study involving animals in accordance with the local legislation and institutional requirements because we only swab the skin samples of two anuran species, which are subsequently released.
Author contributions
HaZ: Conceptualization, Data curation, Writing – original draft, Writing – review & editing. HM: Conceptualization, Writing – original draft. JD: Conceptualization, Writing – original draft. HuZ: Writing – review & editing. CF: Data curation, Writing – review & editing. JZ: Supervision, Writing – review & editing. QW: Methodology, Writing – review & editing. HoZ: Methodology, Writing – review & editing. WJ: Methodology, Writing – review & editing. FK: Supervision, Writing – review & editing.
Funding
The author(s) declare that financial support was received for the research, authorship, and/or publication of this article. This study was supported by the Basic Research Program of Shaanxi Academy of Sciences (2023k-18), Basic Research Program of Shaanxi Academy of Sciences (2023k-10), General social development project of Shaanxi Science and Technology Department (2024SF-YBXM-549), Shaanxi Science and Technology Department, Key R&D Program of Shaanxi Province (2022NY-045), and Shaanxi Science and Technology Department, Key R&D Program of Shaanxi Province (2022NY-101).
Acknowledgments
The authors gratefully thank Beijing Biomarker Technologies Co., Ltd. (Beijing, China) for their assistance in amplicon sequencing.
Conflict of interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Publisher’s note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Supplementary material
The Supplementary material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmicb.2024.1463563/full#supplementary-material
Footnotes
References
Anderson, M. J. (2017). “Permutational multivariate analysis of variance (PERMANOVA)” in inWiley StatsRef: Statistics reference online, American cancer society. eds. N. Balakrishnan, T. Colton, B. Everitt, W. Piegorsch, F. Ruggeri, and J. L. Teugels (Hoboken, NJ: John Wiley & Sons, Ltd), 1–15.
Basanta, M. D., Avila-Akerberg, V., Byrne, A. Q., Castellanos-Morales, G., González-Martínez, T. M., Maldonado-López, Y., et al. (2022). The fungal pathogen Batrachochytrium salamandrivorans is not detected in wild and captive amphibians from Mexico. Peer J. 10:e14117. doi: 10.7717/peerj.14117
Becker, M. H., Brucker, R. M., Schwantes, C. R., Harris, R. N., and Minbiole, K. P. C. (2009). The bacterially produced metabolite violacein is associated with survival of amphibians infected with a lethal fungus. Appl. Environ. Microbiol. 75, 6635–6638. doi: 10.1128/AEM.01294-09
Becker, M. H., Walke, J. B., Murrill, L., Woodhams, D. C., Reinert, L. K., Rollins-Smith, L. A., et al. (2015). Phylogenetic distribution of symbiotic bacteria from Panamanian amphibians that inhibit growth of the lethal fungal pathogen Batrachochytrium dendrobatidis. Mol. Ecol. 24, 1628–1641. doi: 10.1111/mec.13135
Bletz, M. C., Archer, H., Harris, R. N., McKenzie, V. J., Rabemananjara, F. C. E., Rakotoarison, A., et al. (2017a). Host ecology rather than host phylogeny drives amphibian skin microbial community structure in the biodiversity hotspot of Madagascar. Front. Microbiol. 8:1530. doi: 10.3389/fmicb.2017.01530
Bletz, M. C., Perl, R. G. B., Bobowski, B. T., Japke, L. M., Tebbe, C. C., Dohrmann, A. B., et al. (2017b). Amphibian skin microbiota exhibits temporal variation in community structure but stability of predicted Bd-inhibitory function. ISME J. 11, 1521–1534. doi: 10.1038/ismej.2017.41
Bolger, A. M., Lohse, M., and Usadel, B. (2014). Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics 30, 2114–2120. doi: 10.1093/bioinformatics/btu170
Bolyen, E., Rideout, J. R., Dillon, M. R., Bokulich, N. A., Abnet, C. C., Al-Ghalith, G. A., et al. (2019). Reproducible, interactive, scalable and extensible microbiome data science using QUME 2.1. Nat. Biotechnol. 37, 852–857. doi: 10.1038/s41587-019-0209-9
Brucker, R. M., Baylor, C. M., Walters, R. L., Lauer, A., Harris, R. N., and Minbiole, K. P. C. (2008a). The identification of 2, 4-diacetylphloroglucinol as an antifungal metabolite produced by cutaneous bacteria of the salamander Plethodon cinereus. J. Chem. Ecol. 34, 39–43. doi: 10.1007/s10886-007-9352-8
Brucker, R. M., Harris, R. N., Schwantes, C. R., Gallaher, T. N., Flaherty, D. C., Lam, B. A., et al. (2008b). Amphibian chemical defense: antifungal metabolites of the microsymbiont Janthinobacterium lividum on the salamander Plethodon cinereus. J. Chem. Ecol. 34, 1422–1429. doi: 10.1007/s10886-008-9555-7
Brunetti, A. E., Bunk, B., Lyra, M. L., Fuzo, C. A., Marani, M. M., Spröer, C., et al. (2021). Molecular basis of a bacterial-amphibian symbiosis revealed by comparative genomics, modeling, and functional testing. ISME J. 16, 788–800. doi: 10.1038/s41396-021-01121-7
Callahan, B. J., McMurdie, P. J., Rosen, M. J., Han, A. W., Johnson, A. J. A., and Holmes, S. P. (2016). DADA2: high-resolution sample inference from Illumina amplicon data. Nat. Methods 13, 581–583. doi: 10.1038/nmeth.3869
Cambronero-Heinrichs, J. C., Battisti, A., Biedermann, P. H. W., Cavaletto, G., Castro-Gutierrez, V., Favaro, L., et al. (2023). Erwiniaceae bacteria play defensive and nutritional roles in two widespread ambrosia beetles. FEMS. Microbiol. Ecol. 99:fiad144. doi: 10.1093/femsec/fiad144
Costa, S., Lopes, I., Proença, D. N., Ribeiro, R., and Morais, P. V. (2016). Diversity of cutaneous microbiome of Pelophylax perezi populations inhabiting different environments. Sci. Total Environ. 572, 995–1004. doi: 10.1016/j.scitotenv.2016.07.230
Crump, B. C., and Hobbie, J. E. (2005). Synchrony and seasonalityin bacterioplankton communities of two temperaterivers. Limnol. Oceanogr. 50, 1718–1729. doi: 10.4319/lo.2005.50.6.1718
Daskin, J. H., Bell, S. C., Schwarzkopf, L., and Alford, R. A. (2014). Cooltemperatures reduce antifungal activity of symbioticbacteria of threatened amphibians – implications fordisease management and patterns of decline. PLoS One 9:e100378. doi: 10.1371/journal.pone.0100378
Dixon, P. (2003). VEGAN, a package of R functions for community ecology. J. Veg. Sci. 14, 927–930. doi: 10.1111/j.1654-1103.2003.tb02228.x
Edgar, R. C. (2013). UPARSE: highly accurate OTU sequences from microbial amplicon reads. Nat. Methods 10, 996–998. doi: 10.1038/nmeth.2604
Edgar, R. C., Haas, B. J., Clemente, J. C., Quince, C., and Knight, R. (2011). UCHIME improves sensitivity and speed of chimera detection. Bioinformatics 27, 2194–2200. doi: 10.1093/bioinformatics/btr381
Ellison, S., Knapp, R. A., Sparagon, W., Swei, A., and Vredenburg, V. T. (2019b). Reduced skin bacterial diversity correlates with increased pathogen infection intensity in an endangered amphibian host. Mol. Ecol. 28, 127–140. doi: 10.1111/mec.14964
Ellison, S., Rovito, S., Parra-Olea, G., Vásquez-Almazán, C., Flechas, S. V., Bi, K., et al. (2019a). The influence of habitat and phylogeny on the skin microbiome of amphibians in Guatemala and Mexico. Microb. Ecol. 78, 257–267. doi: 10.1007/s00248-018-1288-8
Fraune, S., Anton-Erxleben, F., Augustin, R., Franzenburg, S., Knop, M., Schröder, K., et al. (2015). Bacteria-bacteria interactions within the microbiota of the ancestral metazoan Hydra contribute to fungal resistance. ISME J. 9, 1543–1556. doi: 10.1038/ismej.2014.239
Gray, M. J., Carter, E. D., Piovia-Scott, J., Cusaac, J. P. W., Peterson, A. C., Whetstone, R. D., et al. (2023). Broad host susceptibility of north American amphibian species to Batrachochytrium salamandrivorans suggests high invasion potential and biodiversity risk. Nat. Commun. 14:3270. doi: 10.1038/s41467-023-38979-4
Harris, R. N., Brucker, R. M., Walke, J. B., Becker, M. H., Schwantes, C. R., Flaherty, D. C., et al. (2009). Skin microbes on frogs prevent morbidity and mortality caused by a lethal skin fungus. ISME J. 3, 818–824. doi: 10.1038/ismej.2009.27
Jani, A. J., Bushell, J., Arisdakessian, C. G., Belcaid, M., Boiano, D. M., Brown, C., et al. (2021). The amphibian microbiome exhibits poor resilience following pathogen-induced disturbance. ISME J. 15, 1628–1640. doi: 10.1038/s41396-020-00875-w
Khosravi, A., and Mazmanian, S. K. (2013). Disruption of the gut microbiome as a risk factor for microbial infections. Curr. Opin. Microbiol. 16, 221–227. doi: 10.1016/j.mib.2013.03.009
Kinney, V. C., Heemeyer, J. L., Pessier, A. P., and Lannoo, M. J. (2011). Seasonal pattern of Batrachochytrium dendrobatidis infection and mortality in Lithobates areolatus: affirmation of Vredenburg’s ‘10,000 zoospore rule’. PLoS One 6:e16708. doi: 10.1371/journal.pone.0016708
Kohl, K. D., and Yahn, J. (2016). Effects of environmental temperature on the gut microbial communities of tadpoles. Environ. Microbiol. 18, 1561–1565. doi: 10.1111/1462-2920.13255
Kueneman, J. G., Woodhams, D. C., Harris, R., Archer, H. M., Knight, R., and Mckenzie, V. J. (2016). Probiotic treatment restores protection against lethal fungal infection lost during amphibian captivity. Proc. Biol. Sci. 283:20161553. doi: 10.1098/rspb.2016.1553
Laking, A. E., Ngo, H. N., Pasmans, F., Martel, A., and Nguyen, T. T. (2017). Batrachochytrium salamandrivorans is the predominant chytrid fungus in Vietnamese salamanders. Sci. Rep. 7:44443. doi: 10.1038/srep44443
Longo, A. V., Burrowes, P. A., and Joglar, R. L. (2010). Seasonality of Batrachochytrium dendrobatidis infection in directdeveloping frogs suggests a mechanism for persistence. Dis. Aquat. Org. 92, 253–260. doi: 10.3354/dao02054
Longo, A. V., Savage, A. E., Hewson, L., and Zamudio, K. R. (2015). Seasonal and ontogenetic variation of skin microbial communities and relationships to natural disease dynamics in declining amphibians. R. Soc. Open Sci. 2:140377. doi: 10.1098/rsos.140377
Longo, A. V., and Zamudio, K. R. (2017). Environmental fluctuations and host skin bacteria shift survival advantage between frogs and their fungal pathogen. ISME J. 11, 349–361. doi: 10.1038/ismej.2016.138
Loudon, A. H., Holland, J. A., Umile, T. P., Burzynski, E. A., Minbiole, K. P. C., and Harris, R. N. (2014). Interactions between amphibians’ symbiotic bacteria cause the production of emergent anti-fungal metabolites. Front. Microbiol. 5, 1–8. doi: 10.3389/fmicb.2014.00441
Loudon, A. H., Venkataraman, A., Treuren, W. V., Woodhams, D. C., Parfrey, L. W., McKenzie, V. J., et al. (2016). Vertebrate hosts as islands: dynamics of selection, immigration, loss, persistence, and potential function of bacteria on salamander skin. Front. Microbiol. 7:333. doi: 10.3389/fmicb.2016.00333
Martin, M. (2011). Cutadapt removes adapter sequences from high-throughput sequencing reads. EMBnet J. 17, 10–12. doi: 10.14806/ej.17.1.200
McKenzie, V. J., Bowers, R. M., Fierer, N., Knight, R., and Lauber, C. L. (2012). Co-habiting amphibian species harbor unique skin bacterial communities in wild populations. ISME J. 6, 588–596. doi: 10.1038/ismej.2011.129
Muletz-Wolz, C. R., Almario, J. G., Barnett, S. E., DiRenzo, G. V., Martel, A., Pasmans, F., et al. (2017). Inhibition of fungal pathogens across genotypes and temperatures by amphibian skin bacteria. Front. Microbiol. 8:1551. doi: 10.3389/fmicb.2017.01551
Mutnale, M. C., Reddy, G. S., and Vasudevan, K. (2021). Bacterial Community in the Skin Microbiome of frogs in a Coldspot of Chytridiomycosis infection. Microb. Ecol. 82, 554–558. doi: 10.1007/s00248-020-01669-5
Panthee, S., Hamamoto, H., Paudel, A., and Sekimizu, K. (2016). Lysobacter species: a potential source of novel antibiotics. Arch. Microbiol. 198, 839–845. doi: 10.1007/s00203-016-1278-5
Phillott, A. D., Grogan, L. F., Cashins, S. D., McDonald, K. R., Berger, L., and Skerratt, L. F. (2013). Chytridiomycosis and seasonal mortality of tropical stream-associated frogs 15 years after introduction of Batrachochytrium dendrobatidis. Conserv. Biol. 27, 1058–1068. doi: 10.1111/cobi.12073
Piovia-Scott, J., Rejmanek, D., Woodhams, D. C., Worth, S. J., Kenny, H., McKenzie, V., et al. (2017). Greater species richness of bacterial skin symbionts better suppresses the amphibian fungal pathogen Batrachochytrium dendrobatidis. Microb. Ecol. 74, 217–226. doi: 10.1007/s00248-016-0916-4
Prest, T. L., Kimball, A. K., Kueneman, J. G., and McKenzie, V. J. (2018). Host-associated bacterial community succession during amphibian development. Mol. Ecol. 27, 1992–2006. doi: 10.1111/mec.14507
Quast, C., Pruesse, E., Yilmaz, P., Gerken, J., Schweer, T., Yarza, P., et al. (2013). The SILVA ribosomal RNA gene database project: improved data processing and web-based tools. Nucleic Acids Res. 41, D590–D596. doi: 10.1093/nar/gks1219
Rasche, F., Knapp, D., Kaiser, C., Koranda, M., Kitzler, B., Zechmeister-Boltenstern, S., et al. (2011). Seasonality and resource availability control bacterial and archaeal communities in soils of a temperate beech forest. ISME J. 5, 389–402. doi: 10.1038/ismej.2010.138
Rebollar, E. A., Bridges, T., Hughey, M. C., Medina, D., Belden, L. K., and Harris, R. N. (2019). Integrating the role of antifungal bacteria into skin symbiotic communities of three Neotropical frog species. ISME J. 13, 1763–1775. doi: 10.1038/s41396-019-0388-x
Rebollar, E. A., Hughey, M. C., Medina, D., Harris, R. N., Ibáñez, R., and Belden, L. K. (2016). Skin bacterial diversity of Panamanian frogs is associated with host susceptibility and presence of Batrachochytrium dendrobatidis. ISME J. 10, 1682–1695. doi: 10.1038/ismej.2015.234
Rosenberg, E., Koren, O., Reshef, L., Efrony, L., and Zilber-Rosenberg, I. (2007). The role of microorganisms in coral health, disease and evolution. Nat. Rev. Microbiol. 5, 355–362. doi: 10.1038/nrmicro1635
Scheele, B. C., Pasmans, F., Skerratt, L. F., Berger, L., Martel, A., Beukema, W., et al. (2019). Amphibian fungal panzootic causes catastrophic and ongoing loss of biodiversity. Science 363, 1459–1463. doi: 10.1126/science.aav0379
Segata, N., Izard, J., Waldron, L., Gevers, D., Miropolsky, L., Garrett, W. S., et al. (2011). Metagenomic biomarker discovery and explanation. Genome Biol. 12, 1–18. doi: 10.1186/gb-2011-12-6-r60
Shvartsman, E., Hill, J. E., Sandstrom, P., and MacDonald, K. S. (2023). Gardnerella revisited: species heterogeneity, virulence factors, mucosal immune responses, and contributions to bacterial vaginosis. Infect. Immun. 91:e0039022. doi: 10.1128/iai.00390-22
Towe, A. E., Gray, M. J., Carter, E. D., Wilber, M. Q., Ossiboff, R. J., Ash, K., et al. (2021). Batrachochytrium salamandrivorans can devour more than salamanders. J. Wildl. Dis. 57, 942–948. doi: 10.7589/JWD-D-20-00214
Varela, B. J., Lesbarrères, D., Ibáñez, R., and Green, D. M. (2018). Environmental and host effects on skin bacterial community composition in Panamanian frogs. Front. Microbiol. 9, 1–13. doi: 10.3389/fmicb.2018.00298
Walke, J. B., Becker, M. H., Loftus, S. C., House, L. L., Teotonio, T. L., Minbiole, K. P. C., et al. (2015). Community structure and function of amphibian skin microbes: an experiment with bullfrogs exposed to a chytrid fungus. PLoS One 10:e0139848. doi: 10.1371/journal.pone.0139848
Woodhams, D. C., Alford, R. A., Antwis, R. E., Archer, H., Becker, M. H., Belden, L. K., et al. (2015). Antifungal isolates database of amphibian skin-associated bacteria and function against emerging fungal pathogens. Ecology 96:595. doi: 10.1890/14-1837.1
Woodhams, D. C., Brandt, H., Baumgartner, S., Kielgast, J., Küpfer, E., Tobler, U., et al. (2014). Interacting symbionts and immunity in the amphibian skin mucosome predict disease risk and probiotic effectiveness. PLoS One 9:e96375. doi: 10.1371/journal.pone.0096375
Woodhams, D. C., LaBumbard, B. C., Barnhart, K. L., Becker, M. H., Bletz, M. C., Escobar, L. A., et al. (2018). Prodigiosin, violacein, and volatile organic compounds produced by widespread cutaneous bacteria of amphibians can inhibit two Batrachochytrium fungal pathogens. Microb. Ecol. 75, 1049–1062. doi: 10.1007/s00248-017-1095-7
Woodhams, D. C., Vredenburg, V. T., Simon, M. A., Billheimer, D., Shakhtour, B., Shyr, Y., et al. (2007). Symbiotic bacteria contribute to innate immune defenses of the threatened mountain yellow-legged frog Rana muscosa. Biol. Conserv. 138, 390–398. doi: 10.1016/j.biocon.2007.05.004
Keywords: seasonal variation, host species difference, skin, bacterial community, anurans, amphibians
Citation: Zhang H, Ma H, Deng J, Zhao H, Fang C, Zhang J, Wang Q, Zhang H, Jiang W and Kong F (2024) Seasonality influences skin bacterial community structure and anti-Bd function in two anuran species. Front. Microbiol. 15:1463563. doi: 10.3389/fmicb.2024.1463563
Edited by:
Brianna R. Beechler, Oregon State University, United StatesReviewed by:
Hong Mingsheng, China West Normal University, ChinaSunil Banskar, University of Arizona, United States
Alexandra Alexiev, Oregon State University, United States
Paul Snyder, Yale University, United States
Copyright © 2024 Zhang, Ma, Deng, Zhao, Fang, Zhang, Wang, Zhang, Jiang and Kong. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Fei Kong, k.coffee@163.com