Seasonality influences skin bacterial community structure and anti-Bd function in two anuran species

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).

Table 1. Sample information of anuran species included in this study.

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. 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.

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).

Table 5. The composition and relative abundance of the top five phyla in each groups.

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.

Table 6. The composition and relative abundance of the top five genera in each groups.

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.

Google Scholar

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

PubMed Abstract | Crossref Full Text | Google Scholar

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

Crossref Full Text | Google Scholar

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

PubMed Abstract | Crossref Full Text | Google Scholar

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