Weathering the Storm: Legacies of Extreme Meteorological Events and Daily Weather Variability Shape the Skin Microbiota of the Endangered Golden Alpine Salamander Salamandra atra aurorae (Trevisan, 1982)
Emily L. Pascoe, Federico Polli, Matteo Girardi, Michele Dalponte, Matteo Marcantonio, Antonio Romano, Luca Roner, Giulio Galla, Lucia Zanovello, Paolo Pedrini, Heidi C. Hauffe
TL;DR
Extreme weather events and daily weather changes affect the skin microbes of endangered golden Alpine salamanders, potentially influencing their disease resistance.
Contribution
First characterization of the skin microbiota of Salamandra atra aurorae and its response to a major windstorm and weather variability.
Findings
Salamanders from windstorm-impacted plots had less variable bacterial communities and higher Aeromonas abundance.
Rainfall reduced bacterial diversity but increased fungal beta diversity, showing contrasting moisture preferences.
Weather conditions during sampling explained up to 9% of bacterial and 6% of fungal variation in skin microbiota.
Abstract
Ecosystems worldwide are undergoing unprecedented changes, and as a result amphibians are experiencing devastating population declines driven by subsequent habitat loss and emerging pathogens. The skin microbiota is an important first line of defence for amphibians against pathogens. Here, for the first time, we characterised the bacteria and fungi comprising the skin microbiota of 56 individual golden Alpine salamanders ( Salamandra atra aurorae , Trevisan, 1982), a highly endemic and endangered amphibian subspecies. In addition, we investigated the impact of the 2018 Vaia windstorm on skin microbiota of salamanders in plots classified as impacted or non‐impacted based on windthrows. Salamander sex, weather during sampling, and dominant tree species in plots were also investigated as influencers of microbiota. Beta diversity estimates revealed greater variation in bacterial microbiota…
Genes, proteins, chemicals, diseases, species, mutations and cell lines named across the full text — each resolved to its canonical identifier and authoritative record.
Click any figure to enlarge with its caption.
FIGURE 1
FIGURE 2
FIGURE 3
FIGURE 4
FIGURE 5| Target gene/region | Taxa | Relative abundance | LEfSe | ||
|---|---|---|---|---|---|
| Non‐impacted ( | Impacted ( | ||||
| 16S rRNA | Genus | ||||
|
| 0.01% (±0.09) | < | 0.03% (±0.04) | LDA = 2.12; | |
|
| 0.25% (±1.36) | < | 5.36% (±6.26) | LDA = 4.43; | |
|
| < 0.01% (±0.02) | < | 0.02% (±0.04) | LDA = 2.26; | |
|
| < 0.01% (± < 0.01) | < | 0.02% (±0.04) | LDA = 2.50; | |
|
| 0.47% (±1.53) | < | 4.22% (±4.27) | LDA = 4.27; | |
| ITS1 | Genus | ||||
|
| 0.83% (±1.12) | > | 0.28% (±0.51) | LDA = 3.62; | |
|
| 0.04% (±0.17) | < | 0.65% (±1.48) | LDA = 3.67; | |
|
| 0.49% (±0.65) | > | < 0.01% (±0.01) | LDA = 3.49; | |
|
| 1.02% (±1.74) | > | 0.15% (±0.53) | LDA = 3.68; | |
|
| 0.14% (±0.44) | < | 2.68% (±3.50) | LDA = 4.16; | |
|
| 0.32% (±0.79) | > | < 0.01% (± < 0.01) | LDA = 3.51; | |
|
| 0.90% (±1.41) | > | 0.03% (±0.12) | LDA = 3.74; | |
|
| < 0.01% (± < 0.01) | < | 2.21% (±7.51) | LDA = 4.12; | |
- —HORIZON EUROPE Marie Sklodowska‐Curie Actions10.13039/100018694
Peer Reviews
No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.
Videos
No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.
Taxonomy
TopicsAmphibian and Reptile Biology · Dermatological diseases and infestations · Turtle Biology and Conservation
Introduction
1
Ecosystems worldwide are experiencing unprecedented disturbances that are directly (e.g., land‐use change, introduction of non‐native species, pollution) and indirectly (e.g., extreme weather events linked to anthropogenically induced climate change) associated with human activities (Ummenhofer and Meehl 2017). As the frequency and intensity of these disturbances continue to increase (Gardiner et al. 2013; Ummenhofer and Meehl 2017), it is imperative that we understand the knock‐on effects for the organisms and their microbiota within impacted habitats, particularly endangered taxa, to inform effective conservation and restoration strategies.
Amphibians in particular are extremely sensitive to changes in climate and habitat, in part because they are ectothermic, are sensitive to ultraviolet radiation, have highly permeable skin, and require high air humidity for efficient cutaneous respiration (Alton and Franklin 2017; Blaustein and Kiesecker 2002; Hayden Bofill and Blom 2024). These characteristics also make amphibians susceptible to several emerging pathogens: both Ranaviruses and chytrid fungi (Batrachochytrium spp.) can be transmitted directly between individuals, or through contact with contaminated substrates (e.g., water, soil), and have been implicated in mass mortalities and major amphibian declines worldwide (Gray et al. 2009; Fisher et al. 2021). The skin is an important first line of defence for amphibians against environmental risks, including pathogens and pollution (e.g., Colombo et al. 2015; Cordeiro et al. 2024). More specifically, skin microbiota is crucial in the modulation of the amphibian immune system (Colombo et al. 2015) and skin permeability (Cordeiro et al. 2024). In addition, the anti‐fungal properties of certain bacteria in the skin microbiota of some amphibians provide protection against chytrid fungi (Rebollar et al. 2020), while a skin microbiota with high species richness may be protective against Ranavirus mortality (Harrison et al. 2019). At the same time, habitat disturbances related to agricultural activities (Jiménez et al. 2020), as well as (semi‐)natural phenomena such as wildfires (Mulla and Hernández‐Gómez 2023) and chronic weather anomalies such as drought (Buttimer et al. 2024) may influence the skin microbiota of amphibians. However, it is currently unknown whether habitat disturbances caused by acute extreme weather events also impact amphibian skin microbiota and, by extension, their health.
Just over half of all recorded damage to forests in Europe has been attributed to wind events (Schelhaas 2008; Gardiner et al. 2013). In 2018 the ‘Vaia’ windstorm caused severe windthrow damage to more than 42,000 ha of forest, including across the entire distribution range of an endangered subspecies of Alpine salamander: the golden Alpine salamander Salamandra atra aurorae (Trevisan, 1982; Giannetti et al. 2021). This subspecies occupies patches of fragmented habitat covering no more than 26 km^2^, distributed across a total area spanning 70 km^2^ at 1200–1800 m a.s.l. in the Venetian Prealps (Romanazzi and Bonato 2014). Such a localised distribution, combined with a viviparous reproductive strategy that includes a long gestation period of 2–3 years for the production of just one or two offspring, makes populations slow to recover from sudden and severe changes in habitat and climate (Guirguis et al. 2023). Due to its limited distribution, population fragmentation, and habitat deterioration, the IUCN assessment classifies this salamander as Endangered (EN) under criteria B1ab(iii) + B2ab(iii) (Rondinini et al. 2022). Here, for the first time, we characterise the skin microbiota of S. a. aurorae, using the Vaia windstorm as a unique opportunity to investigate whether severe habitat disturbance impacts the diversity and composition of amphibian skin microbiota and, consequently, resilience to environmental changes (Bahrndorff et al. 2016). At the same time, we investigate if daily weather fluctuations may also impact skin microbiota, an understudied but potentially important and chronic influencer of microbial communities (Varela et al. 2018; Estrada et al. 2019; Ruthsatz et al. 2020). To capture a greater breadth of this microenvironment, we extend beyond the bacterial focus of most skin microbiota studies and also characterise the fungal community, targeting the 16S rRNA gene and ITS1 region, respectively.
Materials and Methods
2
The Extreme Weather Event: The Vaia Windstorm
2.1
Between 27 and 30 October 2018, ‘Vaia’, an extreme meteorological event with wind gusts of almost 200 km/h and more than 600 mm of precipitation in 72 h at some localities, devastated the forests of northeastern Italy (Trenti 2018; Lucianetti et al. 2019; Udali et al. 2021; Vecchiato et al. 2023; Gazzea et al. 2024). Over 42,000 ha of forest are estimated to have been damaged, primarily affecting coniferous forests dominated by Norway spruce ( Picea abies ), silver fir ( Abies alba ), and pine (Pinus spp.), with also some damage in stands of European beech ( Fagus sylvatica ) (Chirici et al. 2019; Giannetti et al. 2021; Udali et al. 2021).
Study Area
2.2
The study area is located at 1450 m a.s.l. on the Vezzena plateau, in the Autonomous Province of Trento (PAT), Italy, within forests affected by Vaia, as described above (Figure 1). The climate is mesothermic‐to‐microthermic, with almost constant humidity (perhumid climate) (Romano et al. 2018). Sampling took place within 7.8 ha of forest, across 33 plots of 20 × 20 m, as part of an ongoing long‐term monitoring program for the subspecies (Romano et al. 2018; Centomo et al. 2023). The size of each plot exceeded that of a typical salamander home range to minimise the likelihood of sampling the same individual in multiple plots (Bonato and Fracasso 2003; Romano et al. 2018).
Overview of area (northern Italy) where golden Alpine salamanders ( Salamandra atra aurorae ) were sampled for skin microbiota analyses. Windthrows caused by the Vaia windstorm occurred in patches throughout the forest. (a) Salamander sampling plots are indicated, with squares representing plots where a skin microbiota sample was collected and triangles where the plot was surveyed but no salamanders were found. Windthrow surveys were conducted at each sampling plot as well as within larger areas surrounding the plots (hatched polygons) to augment environmental data. Plots and polygons in red represent where trees had fallen as a result of the storm (‘impacted’), while those in blue indicate where no trees had fallen (‘non‐impacted’). Representative photos of (b) an impacted and, (c) a non‐impacted plot are provided. Photographs courtesy of A. Romano and L. Roner.
Environmental Variables
2.3
To assess whether tree falls (windthrows) had an impact on the environment, we tested for differences in air temperature, normalised difference vegetation indices (NDVI), and normalised difference water indices (NDWI) between areas impacted and non‐impacted by Vaia. NDVI provides a proxy of land surface ‘greenness’ and, therefore, vegetation cover (Tucker 1979), while NDWI is used to detect changes in surface water content (McFeeters 2013). Since many fallen trees were quickly removed by the Forest Service, we calculated the difference in the percentage of standing trees in each 20 × 20 m salamander sampling plot before (2017) and after (2022) the 2018 storm to provide an assessment of tree loss (Romano et al. 2018). Plots were considered ‘non‐impacted’ if no trees had fallen, while plots with windthrows were classified as ‘impacted’ (Figure 1). To augment the quantity of data we could use to compare temperature, NDVI and NDWI among impacted and non‐impacted forest we also used data from areas nearby to salamander monitoring plots that were surveyed for fallen trees separately, by the PAT Forestry Service (https://siatservices.provincia.tn.it/idt/vector/p_TN_85e14889‐bb85‐4841‐bd8c‐692ce812f554.zip; accessed 27 November 2023); one area we determined manually using PAT classification methods to be non‐impacted (Figure 1). For these additional areas, an impacted polygon was defined as one in which at least 10% of the tree volume (m^3^) was classified as ‘fallen’. The classification of the percentage of fallen trees was obtained from a semi‐automatic classification of images from the Pleiades, Superview, and Geoeye satellites, supported by field surveys carried out by forestry personnel. The total timber volume of a polygon used to estimate the volume of fallen trees was obtained by calculating the average of the harvested timber per hectare of the forest cadastral unit in which the fallen tree was located. Both salamander sampling plots and PAT survey areas were considered impacted by Vaia if fallen trees were recorded after Vaia, and non‐impacted if there were no fallen trees.
Daily average air temperatures (°C) were retrieved from 108 MeteoTrentino meteorological stations distributed throughout PAT for 1 January 2016–29 July 2022. The mean daily temperature data for the centroid of each sampling plot and each windthrow survey underwent spatial interpolation through a linear kriging method, which factored in both location and digital terrain model values. Root mean square error (RMSE) of the daily predictions was computed using the leave‐one‐out cross‐validation strategy over the 108 weather stations, obtaining an RMSE of 1.01°C for the mean daily temperatures. To calculate NDVI and NDWI, satellite imagery encompassing the study area was downloaded from Planet API using functions available in the R package planetR (https://github.com/bevingtona/planetR). All imagery using the PS2, PS2.SD, and PSB.SB sensors within the analytic bundles PlanetScope Ortho Analytic 4B SR (orthorectified, surface reflectance 4‐band, 3 m resolution) and PlanetScope Ortho Analytic 8B SR (orthorectified, surface reflectance 8‐band, 3 m resolution), with < 2% cloud cover, were downloaded for the 2 years preceding Vaia (29 September 2016–25 October 2018) and for the post‐Vaia period up until salamander sampling (15 November 2018–29 July 2022). It was necessary to use imagery from two bundles to cover the time period of interest. All downloaded images were imported in GRASS GIS v.8.3 (Neteler et al. 2012). NDVI and NDWI were calculated using the r.mapcalc GRASS module from the NIR‐red band and NIR‐green band, respectively, using standard formulas (Tucker 1979; McFeeters 2013).
Sampling of
Salamandra atra aurorae Individuals
2.4
Salamanders were captured by gloved hand during the night/morning of 26, 28, and 29 July 2022. Photographs of the unique dorsal patterns of each salamander (Bonato and Fracasso 2003) were taken to avoid sampling an individual multiple times, and to measure total length (as no individuals exhibited tail damage that could compromise length measurements, total length was recorded as the distance from the tip of the snout to the tip of the tail) using ImageJ v.1.53e. Sex was determined following Klewen (1988): individuals were considered male if they had a swollen cloaca, with a rounded outline when seen from the side, and female if they had a flat cloaca. Females were considered pregnant if the posterior part of the trunk appeared enlarged. Salamanders were considered juvenile when total length ≤ 90 mm and adult when the total length > 90 mm (Klewen 1988). Each individual was rinsed with sterile water to eliminate transient microorganisms, and gently swabbed across the body 30 times (following Prado‐Irwin et al. 2017) with a GenoTube Livestock swab (Prionics, Lelystad, NL). Swabs were stored at −20°C until further analyses. All animals were returned to the plot of capture within 6 h of capture and no mortality was observed.
Molecular Methods and Sequencing
2.5
Sterile scissors were used to cut approximately 5 mm from the tip of each swab, and whole DNA was extracted using the DNeasy Blood and Tissue Kit (Qiagen, Valencia, CA, U.S.A.) following the manufacturer's protocol for extraction of DNA from tissue. Two separate PCRs were performed on the resulting DNA using primers which included Illumina adaptors: to identify bacteria, the V3–V4 region of the 16S rRNA gene was amplified using the 341 F and 805 R mod primers (Herlemann et al. 2011), and to identify fungi, the ITS1 region was amplified using the ITS5 F and ITS2 R primers (White et al. 1990). Details regarding primer sequences, and PCR reagents and conditions are available in Table S1. Each amplification product was purified using the MinElute PCR Purification Kit (Qiagen) following the manufacturer's protocol, after which 20 μL was sequenced using paired‐end sequencing (2 × 300 b.p.) on an Illumina MiSeq PE300 (Illumina, San Diego, CA, USA) with a coverage of 100,000 reads per sample.
Bioinformatics
2.6
Quality filtering and processing of sequences were performed using DADA2 v.1.16.0 (Callahan et al. 2016) in R v.4.3.2 (R Core Team 2023). Cutadapt v.3.5 (Martin 2011) was used to remove primer sequences from each read. Sample inference used the ‘pseudo’ pool method and the Divisive Amplicon Denoising Algorithm (DADA; Callahan et al. 2016). An amplicon sequence variant (ASV) table was constructed, and taxonomy was assigned using SILVA v.138.1 (Quast et al. 2013) for bacteria, and UNITE General Fasta v.9.0 for fungi (Abarenkov et al. 2022, 2024). Multiple sequence alignments (MSAs) and phylogenetic tree construction were performed in MICCA v.1.7.2 (Albanese et al. 2015). For 16S rRNA gene data, we performed a template‐based MSA using the Nearest Alignment Space Termination (NAST) algorithm (DeSantis et al. 2006) and the SILVA reference template, while for the ITS1 region we used the MUltiple Sequence Comparison by Log‐Expectation (MUSCLE) algorithm with five iterations. Decontamination of libraries was accomplished using decontam v.1.24.0 (Davis et al. 2017) and negative controls from both DNA extractions and PCRs were used to identify contaminants: 42 ASVs in the 16S rRNA gene dataset, and 19 ASVs in the ITS1 region dataset. Non‐target sequences were filtered out: sequences originating from mitochondria or chloroplasts, and any sequences not taxonomically assigned at the kingdom level to archaea or bacteria were removed from the 16S rRNA gene dataset, and for the ITS1 dataset any sequences not taxonomically assigned at the kingdom level to fungi were removed. Normalisation was performed using the scaling with ranked subsampling (SRS) method, which minimises the likelihood of discarding rare ASVs (Beule and Karlovsky 2020; see Table S2 for more details regarding bioinformatic parameters).
Statistical Analyses
2.7
Time series data for air temperature, NDVI, and NDWI were de‐trended to remove temporal trends that may influence or bias statistical models using the forecast v.8.23.0 package in R (Hyndman and Khandakar 2008). We fitted an exponential smoothing state‐space model (TBATS) to each time series, then obtained the model residuals by subtracting model fitted values from observed values. All de‐trended variables were averaged spatially per 20 × 20 m plot. Linear mixed models (LMMs), were used to test for differences between impacted and non‐impacted plots, where the response variable was the mean of the de‐trended time series of the plot for either air temperature, NDVI, or NDWI (i.e., three LMMs, a separate LMM for each variable). Fixed effects included plot type (non‐impacted or impacted), and Vaia time‐point (before or after the Vaia storm), plus the interaction between these two variables. Random effects included plot ID and Julian day to integrate potential spatial and temporal differences between plots unaccounted for by the covariates and de‐trending. In addition, LMMs were used to test for differences among plots where salamanders were present or absent (thus excluding environmental data derived from PAT windthrow survey areas, as these were not checked for salamanders) for air temperature, NDVI, or NDWI. Salamander presence (present or absent) was the fixed effect and random effects included plot ID and Julian day. Tukey post hoc tests were performed to determine statistical outputs for each level and combination thereof of the interacting factors.
A generalised linear mixed model (GLMM) was used to compare salamander counts, including zeros from plots with no salamanders, between impacted and non‐impacted plots. The model used the Poisson error distribution and log link function, included plot ID and sampling date as a random factor, accounted for zero inflation, and was implemented using the glmmADMB package (v0.8.3.3; Fournier et al. 2012). Fisher's exact tests were used to test for an association between plot type (impacted vs. non‐impacted) and sex ratio (males vs. females), as well as the number of gravid compared to non‐gravid females. LMMs were used to test for associations between plot type, (i) total salamander length, and (ii) salamander mass, and included sex, plot ID, and sampling date as random factors. Variance components were estimated using the Restricted Maximum Likelihood (REML) method in the lme4 package (v1.1.35.5; Bates et al. 2015), and p‐values for fixed effects were obtained using the lmerTest package (v3.1.3; Kuznetsova et al. 2017).
The linear discriminant analysis (LDA) effect size (LefSe; Segata et al. 2011) method was used to test for differential relative abundances of bacterial and fungal genera in microbiota profiles from impacted and non‐impacted plots. Preliminary analyses indicated a significant difference in the relative abundance of a potential amphibian pathogen (Aeromonas spp.) among plot types. To improve taxonomic resolution for this genus we reassigned taxonomy using the species training set of Greengenes2 v.2024.09. We then performed a Mann–Whitney U test on the relative abundances of Aeromonas hydrophila (the known pathogenic species within this genus), including samples where relative abundance of A. hydrophila was zero, to test for differences between impacted and non‐impacted plots.
For skin microbiota alpha diversity, the Chao1, Shannon entropy, inverse Simpson, and Faith's phylogenetic diversity indices were calculated for both the 16S rRNA gene and ITS1 datasets using the microeco package v1.8.0 in R (Liu et al. 2021). The exponential of Shannon entropy estimates were used in analyses to facilitate interpretations in terms of diversity (Jost 2006). For visualisation purposes, all four alpha diversity indices were compared between non‐impacted and impacted plots using box plots and Mann–Whitney U tests, for both the 16S rRNA gene and ITS1 sequences, and *p‐*values were adjusted using the Benjamini‐Hochberg correction to account for comparisons of multiple diversity indices.
We used a single GLMM for each alpha diversity estimate for both the 16S rRNA gene and ITS1 region to test for associations with plot type (impacted or non‐impacted), weather during sampling (sunny, sunny post‐rain, or strong‐to‐light rain), dominant tree species within the plot ( A. alba , F. sylvatica , or P. abies ), and sex (male or female; a single juvenile was sampled but was excluded from this analysis). GLMMs used the Gamma error distribution with inverse (Chao1, Faith's phylogenetic diversity) or log (Shannon entropy, inverse Simpson) function. For each categorical variable included in a GLMM with more than two levels, the emmeans package in R was used to perform pairwise comparisons, adjusting for multiple comparisons, to estimate the marginal means and test for differences between groups (Lenth 2024).
Skin microbiota beta diversity was estimated as Bray–Curtis dissimilarity, Jaccard, and both unweighted and weighted UniFrac indices using microeco, for both the 16S rRNA gene and ITS1 region datasets. Permutational multivariate analyses of variance (PERMANOVA) were performed to test for differences in microbiota by including plot type (impacted or non‐impacted), weather during sampling (sunny, sunny post‐rain, or strong‐to‐light rain), dominant tree species within the plot ( A. alba , F. sylvatica , or P. abies ), and sex (male or female, the single juvenile was excluded from this analysis) on each beta diversity estimate, followed by Permutational Analysis of Multivariate Dispersions (PERMDISP) to test the homogeneity of dispersion for each of these variables. Mantel tests (using Pearson's correlation coefficients and 9999 permutations) and distance‐based redundancy analyses (db‐RDA) were used to explore spatial associations in microbiota composition (i.e., distance between sampled individuals and microbiota composition) based on Bray–Curtis dissimilarity estimates.
Results
3
Environmental Variables
3.1
There was no significant difference in mean daily air temperature associated with plot type (across all years), between the periods before and after Vaia (for all plots), nor when plot type and Vaia time period were considered together as interacting factors. NDVI was significantly higher (t = 9.96; S.E. ≤ 0.01; p < 0.01) before Vaia (2016–2018; NDVI = 0.75 ± 0.08 S.D.), compared with after (2018–2022; NDVI = 0.66 ± 0.20 S.D.) for all plots, regardless of impact type. However, when plot type was taken into account, no differences in NDVI were detected, neither before nor after Vaia, nor for both periods combined. When all plots were considered together, NDWI was significantly lower (t = −10.88; d.f. = 12,530; p < 0.01) after Vaia (NDWI = −0.68 ± 0.19 S.D.) compared to before (NDWI = −0.82 ± 0.08 S.D.), but plot type had no effect: impacted and non‐impacted plots did not differ before Vaia, after Vaia, or when considered across both periods. Details regarding the mean values of environmental variables and the results from these GLMMs are available in Table S3.
Mean air temperature and NDVI did not differ among plots with or without salamanders. However, NDWI was lower (indicating a drier environment) than expected based on the overall trend at plots where salamanders were present (salamander presence = −0.71 ± 0.18, absence = −0.68 ± 0.19; t = −1.99; S.E. ≤ 0.01; p = 0.05). Details regarding the results from these GLMMs are available in Table S3.
Salamandra atra aurorae
Samples
3.2
A total of 56 golden Alpine salamanders were sampled from 14 out of 33 surveyed plots: 44 individuals were from 12 non‐impacted plots, and 12 were from two plots impacted by Vaia (with 17.6% and 52.1% windthrows, respectively). No salamanders were found at 19/33 plots, 10 of which were non‐impacted, and nine impacted, where windthrows ranged from 24.2% to 74.4%. However, the number of salamanders (including plots with zero salamanders) found was not associated with plot type (Z = −0.25; S.E. = 2.8; p = 0.80). Of the captured salamanders, 35 were male (non‐impacted = 63.6%; 28/44; impacted = 58.3%; 7/12), and 20 were female. Visibly gravid females (non‐impacted = 18.2%; 8/44; impacted = 25.0%; 3/12) were more numerous than those showing no signs of being gravid (non‐impacted = 15.9%; 7/44; impacted = 16.7%; 2/12). One juvenile was captured for which sex could not be determined (Table S4). The salamander sex ratio did not differ between impacted and non‐impacted plots (O.R. = 0.8; 95% C.I. = 0.2, 3.6; p = 0.74; Table S4), nor did the proportion of gravid compared to non‐gravid females (O.R. = 1.3; 95% C.I. = 0.1, 19.9; p = 1; Table S4). Accounting for variation associated with sex, there was no difference in body mass between salamanders from non‐impacted plots (mean = 9.0 g; range = 4.2–15.0 g) and impacted plots (mean = 9.1 g; range = 6.0–12.9 g; t = 0.07; S.E. = 0.63; p = 0.95). Similarly, there was no difference in total body length between individuals from non‐impacted (mean = 123.3 mm; range = 94.0–141.5 mm) and impacted plots (mean = 123.9 mm; range = 107.8–134.8 mm) (t = 0.12; S.E. = 3.3; p = 0.91).
Skin Microbiota Composition
3.3
Following sequence trimming, filtering, and quality control 6,596,108 reads were obtained, of which 5,616,479 (reads per sample: median = 104,561; ±29,082.7 S.D.) were bacterial, and 979,629 (median = 7,022; ±21,829.6 S.D.) were fungal. Following removal of contaminants, normalisation, and classification, 12,397 bacterial ASVs were identified from 56 samples, and 2149 fungal ASVs were detected in 48 samples (some samples were removed from the analysed dataset due to low number of sequences).
Proteobacteria was by far the most abundant bacterial phylum (70.0%), followed by Bacteroidota (10.50%) and Actinobacteriota (9.0%) (Figures 2 and S1). The remaining 39 identified phyla each constituted < 5% of the bacterial microbiota, while < 0.01% of reads could not be classified at phylum level. A total of 759 genera were identified, but 11.9% of reads could not be classified to this taxonomic level; 19.3% of reads were Acinetobacter, 17.1% were Pseudomonas, and 5.3% were Pedobacter (Figure 2). All other genera were present in < 5% reads each. LEfSe analyses identified significant differences in the abundances of bacteria in the genera Acidicapsa, Aeromonas, Alloprevotella, Oerskovia, and Pseudarthrobacter between plot types, all of which were higher in relative abundance in the skin microbiota of salamander from impacted compared to non‐impacted plots (Table 1). Following reassignment of Aeromonas spp. taxonomy to species level using Greengenes2, we found that median relative abundances of A. hydrophila were higher in the microbiota of salamanders from impacted (0.01 ± 0.02 S.D.) compared with non‐impacted (0 ± < 0.01 S.D.) plots (W = 109; Hodges–Lehmann median difference = < −0.01; 95% C.I. = −0.02, < −0.01; p < 0.01).
Relative abundance of (a) 16S rRNA gene, and (b) ITS1 amplicon sequence variants at the taxonomic level of (1) phylum, and (2) genus in the skin microbiota of golden Alpine salamanders ( Salamandra atra aurorae ) sampled at plots that were either impacted or non‐impacted by the Vaia storm.
TABLE 1: Overview of bacterial genera with significantly different relative abundances (mean with standard deviation in brackets) in the skin microbiota of golden Alpine salamanders ( Salamandra atra aurorae ) sampled at plots impacted by the Vaia storm, compared with non‐impacted plots. Results of the linear discriminant analysis effect size (LEfSe) method are provided.
Regarding fungi, 14 phyla were identified using the ITS1 gene; 64.4% of reads were Ascomycota and 22.5% were Basidiomycota (Figures 2 and S1). Approximately 10% of reads could not be classified at this taxonomic level (Figure 2). We identified 397 genera of fungi, but more than a quarter (26.5%) of reads could not be classified to genus (Figure 2). Overall, the most common fungal genera were Lophodermium (11.1%) and Cladosporium (6.4%), although there was considerable inter‐individual variation (Figure S2). No fungi in the genus Batrachochytrium (chytrid fungi) were detected in any individual, but 0.35% of reads were identified as the Order Rhizophydiales (the same Order to which Batrachochytrium spp. belong) in a salamander from a non‐impacted plot; however, no clinical signs of disease were noted for this individual. Three fungal genera were present in higher relative abundances in the skin microbiota of salamanders from impacted compared to non‐impacted plots (Herpotrichia, Lachnum, and Rhizidium), while five genera were in lower relative abundances (Exobasidium, Hormonema, Lachnellula, Parafenestella, and Piskurozyma; Table 1).
Alpha Diversity of Skin Microbiota
3.4
Alpha diversity of fungal skin microbiota tended to be higher for salamanders from non‐impacted plots, though this difference was not statistically significant (Figure 3). We did not find any associations between bacterial alpha diversity and plot type, dominant tree species, or sex. In Tukey post hoc tests Chao1 (z ratio = 2.59; S.E. > 0.01; p = 0.03) and Faith's phylogenetic diversity metrics (z ratio = 2.46; S.E. = 0.01; p = 0.04) indicated that bacterial richness and phylogenetic diversity were higher during strong‐to‐light rain as compared to sunny sampling conditions. Mann–Whitney U tests revealed a near‐significant association between plot type and variation in fungal community evenness, as measured by the inverse Simpson index (p = 0.05, Figure 3, Table S5). However, plot type was no longer significant (t = −1.56; S.E. = 0.26; p = 0.12) after accounting for weather, dominant tree species, and sex.
Box plots of alpha diversity metrics for (a) 16S rRNA gene, and (b) ITS1 region amplified from the skin microbiota of golden Alpine salamanders ( Salamandra atra aurorae ). Sampling plots were categorised as non‐impacted (blue) or impacted (red) by the Vaia storm, based on windthrows. The p‐values (black, and in bold where ≤ 0.05) and Benjamini–Hochberg (false discovery rate) adjusted p‐values (grey) from Mann–Whitney U tests comparing the median alpha diversity values for salamanders from non‐impacted and impacted plots, are reported for each metric.
Beta Diversity
3.5
Estimates of all four beta diversity indices from skin microbiota were significantly higher for salamanders from non‐impacted compared with impacted plots, even after adjusting p‐values for multiple comparisons (Figure 4). Weather impacted beta diversity estimates for both bacteria and fungi: bacterial beta diversity was higher when conditions were sunny (including post‐rain), whereas we observed the opposite trend for fungi, for which sunny conditions were associated with lower beta diversity estimates (Figure S3). Beta diversity was typically higher in females compared to males (Figure S4). For both Bray‐Curtis and Jaccard indices, bacterial beta diversity was lower in plots dominated by A. alba compared with F. sylvatica or P. abies , but fungal beta diversity was higher for plots where A. alba was dominant (Figure S5).
Box plots of beta diversity estimates for (a) 16S rRNA, and (b) ITS1 genes amplified from the skin microbiota of golden Alpine salamanders ( Salamandra atra aurorae ). Sampling plots were categorised as non‐impacted (blue) or impacted (red) by the Vaia storm, based on windthrows. The p‐values (black, and in bold where ≤ 0.05) and Benjamini–Hochberg (false discovery rate) adjusted p‐values (grey) from Mann–Whitney U tests comparing the median beta diversity values for salamanders from non‐impacted and impacted plots, are reported for each estimate.
PERMANOVA analyses revealed that weather had the strongest and most consistent association with skin microbiota beta diversity, explaining up to 9% of bacterial variation and 6% of fungal variation, as well as being significantly associated with all beta diversity estimates, except for unweighted UniFrac estimates for bacteria (Table S6). In contrast, the impact of the Vaia storm (plot type: impacted vs. non‐impacted plots) was significantly associated with the bacterial community only (p adjusted = 0.01–0.04), explaining 3% of the variation in Bray‐Curtis, Jaccard, and Unweighted UniFrac distances (but not weighted UniFrac:p adjusted = 0.46). Dominant tree species had a marginal effect on both bacteria and fungi, but following corrections for multiple index comparisons, significance was only observed in Bray‐Curtis and Jaccard distances for fungi (p = 0.03). Sex did not significantly explain variation in skin microbiota composition.
Analyses of multivariate dispersion (PERMDISP) showed significant effects of weather conditions on microbiota dispersion for bacteria and fungi (p adjusted < 0.01–0.03), for all estimates except weighted UniFrac of bacteria (p adjusted = 0.08; Figure 5; Table S7). Dispersion also differed between impacted and non‐impacted plots, particularly for bacteria with respect to Bray‐Curtis (F = 4.54, p adjusted = 0.05), Jaccard (F = 6.25; p adjusted = 0.05), and Unweighted UniFrac distances (F = 4.96; p adjusted = 0.05). For fungi, dispersion differed between plot types only for Jaccard estimates (F = 10.17; p adjusted = 0.01; Figure 5). Dominant tree species influenced dispersion of fungal (p adjusted < 0.01, except for weighted UniFrac = 0.84) but not bacterial (p adjusted ≥ 0.90) beta diversity. Sex had no significant effect across any estimate for bacteria or fungi (p adjusted ≥ 0.62).
Principal coordinate analyses of four beta diversity estimates calculated for (a) 16S rRNA gene, and (b) ITS1 genes amplified from skin microbiota of golden Alpine salamanders ( Salamandra atra aurorae ). These salamanders were sampled at plots, categorised as non‐impacted (blue) or impacted (red) by the Vaia storm based on windthrows.
Across the whole study area (regardless of plot type) there was a significant positive correlation between plot proximity and skin microbiota similarity (based on Bray‐Curtis distances) for bacteria (Mantel test: r = 0.31; p < 0.01, 999 permutations), but not for fungi (r = 0.03, p = 0.33, 999 permutations). Distance‐based redundancy analyses indicated that latitude and longitude together explained 4.8% of the variation in bacterial skin microbiota, representing a small but significant spatial effect (F = 2.40, p < 0.01, 999 permutations; adjusted R ^2^ = 0.05). Both latitude (F = 2.31, p < 0.01) and longitude (F = 2.48, p < 0.01) independently contributed similarly to this pattern. Instead, there was no geographical structuring of the fungal community of the skin; latitude and longitude together accounted for less than 0.5% of the variation (F = 1.13, p = 0.08, 999 permutations; adjusted R ^2^ = 0.01), and neither latitude (F = 1.13, p = 0.13) nor longitude (F = 1.13, p = 0.13) had a significant individual effect.
Discussion
4
The golden Alpine salamander is endemic to just 26 km^2^ of patchy habitat on the Vezzena and Asiago mountain plateaus (Romanazzi and Bonato 2014). Long term monitoring efforts have highlighted that the number of salamanders in the study area has dramatically decreased since the Vaia windstorm (by up to 96% in 2021), such that there were too few individuals to be able to estimate population size (Roner et al. 2022). Since microbiota is known to be essential for health, here for the first time, we provide a detailed analysis of the skin microbial community of this endangered salamander, which we found to vary considerably among individuals. Extreme meteorological events as well as daily weather conditions impacted salamander skin microbiota, including the relative abundance of potentially pathogenic bacteria. These findings underscore the need to integrate microbiota research into conservation strategies, particularly for species with restricted distributions and slow reproductive rates.
The most abundant phyla detected in the golden Alpine salamander skin microbiota (Proteobacteria, Bacteroidota and Actinobacteriota) were comparable to those of the closely related fire salamander ( S. salamandra ; Hill et al. 2010), although relative dominance of these phyla varied. We detected significantly higher abundances of five bacterial and three fungal genera in the skin microbiota of golden Alpine salamanders from impacted compared to non‐impacted forest plots. Notably, relative abundances of A. hydrophila were significantly higher at impacted plots: this is extremely concerning as this bacteria is the causative agent of red‐leg syndrome in amphibians, an infection that can cause widespread mortality (Hill et al. 2010; Pastorino et al. 2023). This result suggests that habitat disturbances may facilitate opportunistic bacterial colonisation, potentially increasing disease susceptibility. Similar trends have been observed in other amphibian populations, where environmental stressors disrupted skin microbial homeostasis, leading to increased pathogen load (e.g., the fire salamander, Salamandra salamandra ; Hill et al. 2010), and green‐eyed frog, Lithobates vibicarius (Jiménez et al. 2020). Other microorganisms that showed significant differences in relative abundances at impacted plots may be linked to changes in forest structure caused by Vaia. For example, the bacterial genus Oerskovia and the fungal genera Hormonema and Lachnellula are associated with decaying plant material, which may have increased in availability following windthrows (Becker et al. 2017; Grasselli et al. 2019; Dondero et al. 2023). Although the relevance and implications of the other genera remain unclear, these results indicate that Vaia may have influenced the skin microbiota of this salamander subspecies, with potentially adverse health effects. Interestingly, while impacted plots showed significant shifts in bacterial diversity, the fungal microbiota exhibited a different pattern. Certain fungal genera, such as Lachnum and Rhizidium, were more abundant in impacted plots, whereas others, including Exobasidium and Piskurozyma, were less so. These shifts may reflect changes in microhabitat conditions, such as increased moisture retention in decomposing logs or altered competition dynamics between microbial species. The absence of Batrachochytrium spp., as previously reported in another study from the same area focussing on B. salamandrivorans (Bsal; Preuss et al. 2020), is reassuring; however, the presence of other fungal taxa within the order Rhizophydiales warrants further investigation. The susceptibility of the golden Alpine salamander to Bsal infection is currently unknown, but is assumed to be high as closely related species, such as the Alpine salamander ( S. atra ) and fire salamander, are highly susceptible to this pathogen (Turner et al. 2003; Schuck et al. 2024). Although several studies have screened for Bsal in urodeles in Italy, including in the golden Alpine salamander (Grasselli et al. 2019; Dondero et al. 2023; Romano et al. 2024), as yet no dedicated, systematic surveillance initiatives have been established, an oversight given the potential severity of the disease and its possible impact on populations of endangered species.
Our data indicated a trend towards lower fungal alpha diversity at impacted plots, suggesting that (together with beta diversity metrics) these conditions may have favoured a few dominant fungal taxa, leading to reduced community evenness. However, differences in alpha diversity were not significant, possibly as few salamanders were present in impacted plots. Lower bacterial alpha diversity in skin microbiota associated with habitat disturbances, e.g., agricultural and other human‐driven activities, has been reported in amphibians (the lesser treefrog, Dendropsophus minutus ; Becker et al. 2017; Preuss et al. 2020; green‐eyed frog; Jiménez et al. 2020). In contrast, natural disturbances like wildfires, even when simulated, appear to have less impact on alpha diversity due to their temporary and infrequent nature (rather than being chronic and continuous, such as in agricultural impacts), leading to differing selection pressures (Turner et al. 2003; Schuck et al. 2024). While we studied the potential impact of Vaia at the microhabitat level (microbiota), studies on other taxa at the mesohabitat scale (invertebrate communities) found higher species richness in Vaia‐damaged plots (Gazzea et al. 2024), likely due to new species colonising empty niches exposed by habitat changes (Seidl et al. 2022). Collectively, this highlights that impacts of disturbance may vary across scales, with factors such as disturbance type, severity, and timing influencing outcomes. Further studies are needed to understand ecological trends at the microhabitat level, including microbial communities (Hobbs and Huenneke 1992; Seidl et al. 2022).
Amphibian skin microbiota is likely acquired, at least in part, from the environment (Loudon et al. 2014; Costa et al. 2016). In contrast to previous studies (Hill et al. 2023), we detected fine geographical scale correlations in skin microbiota; specifically, plot proximity and similarity of the bacterial skin microbiota were positively correlated. In addition, our results demonstrate that both transient environmental fluctuations and stable, long‐term habitat features may influence salamander skin microbiota. Daily weather conditions in particular were consistently associated with differences in skin microbiota: bacterial diversity was higher during sunny compared with rainy conditions for nearly all PERMANOVA metrics, indicating changes in both composition but also relative abundances of bacterial taxa, thus microbiota is likely highly responsive to humidity and precipitation (Woodhams et al. 2016). Fungal beta diversity metrics, in contrast, were higher during rainy conditions, possibly due to the displacement of spores or changes in the condition of the salamander's skin. This is of no surprise as fungi require water to complete their lifecycle, with high relative humidity critical for spore germination (Herman and Bleichrodt 2022). In contrast, bacteria adhesion on drier surfaces tends to be increased (Oh et al. 2018).
Habitat structure and the impact of the Vaia windstorm were also associated with salamander skin microbiota composition, albeit to a lesser extent than weather, with low‐abundance bacterial taxa mainly impacted and little effect on the relative abundance of the dominant taxa (weighted UniFrac). More specifically, bacterial beta diversity at impacted plots was generally lower than at non‐impacted plots. This result was surprising as stressors are generally considered to be associated with greater variability in microbiota composition (Zaneveld et al. 2017). The environmental bacterial community to which the salamander is exposed may be affected by habitat alterations caused by tree falls, such as changes in microclimate, resource availability, or subsequent alterations in host behaviour. However, since bacterial microbiota are known to be closely associated or even selected for by the host (Walke et al. 2014; Liu et al. 2023), it may be that changes in bacterial microbiota are indicative of adaptability and resilience to changing conditions. Fungal communities were less affected by the impacts of Vaia, suggesting that they might be more strongly structured by longer‐term habitat characteristics, such as dominant tree species (Estrada et al. 2019). Plots dominated by European beech ( F. sylvatica ) typically supported higher bacterial but lower fungal beta diversity. Differences in the bacteria community were observed in the most dominant taxa only (weighted UniFrac, but p adjusted > 0.05), while for fungi there was a difference in the presence/absence of closely related taxa (Bray‐Curtis, Jaccard). This is likely a result of differences in the composition of the soil, as well as other substrates such as leaves and bark, that salamanders would come into contact with. For example, Norway spruce ( Picea abies ) causes soil acidification and can alter other soil parameters, such as chemical composition (Augusto et al. 1998), potentially impacting the microorganisms able to survive within such conditions (Sullivan et al. 2017).
Despite windthrows of more than 70% of trees within the study area, we did not detect any difference in environmental variables among impacted and non‐impacted plots. We may not have detected climatic differences among plot types as air temperature values were modelled from data acquired from weather stations, while NDVI and NDWI values were calculated using satellite imagery. Such data has limitations due to the interpolation process and the coarse spatial scale, which may hinder their capacity to detect fine‐scale microhabitat variation. In addition, windthrows can dramatically alter environmental conditions that were not measured in the current study, primarily through canopy opening. For example, shifts in exposure to ultraviolet radiation, vapour‐pressure deficit, and soil composition and chemistry resulting from changes in woody debris and leaf litter availability are consequences of windthrows or reduced canopy cover (Simon et al. 2011; Zhang et al. 2024; Waring et al. 2025). Such effects can persist for years and may in turn impact amphibian microbiota and host‐pathogen dynamics (Bernardo‐Cravo et al. 2020). Nevertheless, both NDVI and NDWI were significantly higher before the storm compared to after, suggesting that drastic changes to forest structure may influence habitat and mesoclimate conditions across a wide‐ranging area, potentially affecting salamander home ranges, even when windthrows did not occur directly within them. In fact, Roner et al. (2022) demonstrate that golden Alpine salamander numbers have been notably lower since the Vaia storm (96% reduction in the number of salamanders in 2021, compared to pre‐Vaia numbers), with salamanders consistently absent from plots that were occupied before the storm. This is of no surprise, as this species tends to occupy inner parts of the forest, is rarely found in more open areas, and is sensitive to ground perturbations and compaction (Romanazzi and Bonato 2014). Instead the golden Alpine salamander relies on leaf litter, coarse woody debris, and other microhabitats that can provide shelter and humidity, nearly all of which were removed or degraded by large machinery removing windthrows after the storm. Importantly, we did not find any differences in salamander population metrics (number per plot, sex ratio, mass and body length) among plot types, indicating that those individuals remaining were likely healthy and able to adapt to habitat changes.
Our findings highlight that both short‐term weather fluctuations and past extreme meteorological events may shape the skin microbiota of amphibians, but their impacts differ in magnitude and across microbial taxa. Daily fluctuations appear to exert immediate and dynamic influences on microbial communities, whereas the legacy effects of disturbances such as storms may contribute to long‐term changes, particularly for bacterial microbiota. Given the decreasing numbers, localised distribution and slow reproductive strategy of S. a. aurorae, any perturbations to its habitat pose significant conservation risks, as well as to ecosystem functions, especially given its generalist predatory feeding strategy (Centomo et al. 2023). Our study emphasises the need to monitor not only population trends but also microbial health as an indicator of environmental stability. Conservation strategies should consider the potential effects of extreme weather events and microclimatic conditions on salamander‐associated microbiota. Mitigating habitat disturbances through controlled forest management which maintains stable microhabitats with sufficient moisture levels could help to support beneficial microbial communities. Furthermore, future studies should investigate whether shifts in microbiota composition translate to changes in disease susceptibility or overall fitness. Longitudinal monitoring of the skin microbiota, coupled with pathogen screening, could provide early warning signs of emerging threats to salamander populations.
Author Contributions
Emily L. Pascoe: conceptualization (lead), data curation (equal), formal analysis (lead), funding acquisition (equal), investigation (equal), methodology (equal), project administration (lead), software (lead), validation (equal), visualization (lead), writing – original draft (lead), writing – review and editing (lead). Federico Polli: conceptualization (supporting), data curation (equal), investigation (equal), validation (equal), writing – review and editing (equal). Matteo Girardi: conceptualization (supporting), data curation (supporting), investigation (supporting), methodology (supporting), resources (supporting), supervision (supporting), validation (supporting). Michele Dalponte: data curation (equal), formal analysis (equal), investigation (supporting), methodology (supporting), software (equal), validation (equal), writing – review and editing (equal). Matteo Marcantonio: data curation (equal), formal analysis (equal), investigation (supporting), methodology (supporting), software (equal), validation (equal), writing – review and editing (equal). Antonio Romano: conceptualization (supporting), data curation (equal), funding acquisition (supporting), investigation (equal), methodology (supporting), validation (equal), visualization (supporting), writing – review and editing (equal). Luca Roner: conceptualization (supporting), data curation (equal), investigation (equal), methodology (supporting), validation (equal), writing – review and editing (equal). Giulio Galla: software (supporting), validation (supporting), writing – review and editing (equal). Lucia Zanovello: investigation (supporting), writing – review and editing (equal). Paolo Pedrini: funding acquisition (supporting), resources (supporting), supervision (supporting). Heidi C. Hauffe: conceptualization (supporting), funding acquisition (equal), project administration (supporting), resources (supporting), supervision (equal), writing – review and editing (equal).
Funding
Data analyses were funded by the European Union's Horizon Europe research and innovation programme under the Marie Skłodowska‐Curie grant agreement no. 101067351 (Project NIPMAP).
Ethics Statement
Permits for capture, handling and swabbing of salamanders were issued by the Italian Ministry of Ecological Transition—authorization MiTE‐0014200 of 7 February 2022.
Consent
The authors have nothing to report.
Conflicts of Interest
The authors declare no conflicts of interest.
Supporting information
Data S1: ece372981‐sup‐0001‐DataS1.doc.
The reference list from the paper itself. Each links out to its DOI / PubMed record.
- 1Abarenkov, K. , R. H. Nilsson , K.‐H. Larsson , et al. 2024. “The UNITE Database for Molecular Identification and Taxonomic Communication of Fungi and Other Eukaryotes: Sequences, Taxa and Classifications Reconsidered.” Nucleic Acids Research 52: D 791–D 797. 10.1093/nar/gkad 1039.37953409 PMC 10767974 · doi ↗ · pubmed ↗
- 2Abarenkov, K. , A. Zirk , T. Piirmann , et al. 2022. “UNITE General FASTA Release for Fungi.” 10.15156/BIO/2483911. · doi ↗
- 3Albanese, D. , P. Fontana , C. D. Filippo , D. Cavalieri , and C. Donati . 2015. “MICCA: A Complete and Accurate Software for Taxonomic Profiling of Metagenomic Data.” Scientific Reports 5: 9743. 10.1038/srep 09743.25988396 PMC 4649890 · doi ↗ · pubmed ↗
- 4Alton, L. A. , and C. E. Franklin . 2017. “Drivers of Amphibian Declines: Effects of Ultraviolet Radiation and Interactions With Other Environmental Factors.” Climatic Change Responses 4: 6. 10.1186/s 40665-017-0034-7. · doi ↗
- 5Augusto, L. , P. Bonnaud , and J. Ranger . 1998. “Impact of Tree Species on Forest Soil Acidification.” Forest Ecology and Management 105: 67–78. 10.1016/S 0378-1127(97)00270-3. · doi ↗
- 6Bahrndorff, S. , T. Alemu , T. Alemneh , and J. Lund Nielsen . 2016. “The Microbiome of Animals: Implications for Conservation Biology.” International Journal of Genomics 2016: 5304028. 10.1155/2016/5304028.27195280 PMC 4852354 · doi ↗ · pubmed ↗
- 7Bates, D. , M. Mächler , B. Bolker , and S. Walker . 2015. “Fitting Linear Mixed‐Effects Models Using lme 4.” Journal of Statistical Software 67: i 01. 10.18637/jss.v 067.i 01. · doi ↗
- 8Becker, C. G. , A. V. Longo , C. F. B. Haddad , and K. R. Zamudio . 2017. “Land Cover and Forest Connectivity Alter the Interactions Among Host, Pathogen and Skin Microbiome.” Proceedings of the Royal Society B: Biological Sciences 284: 20170582. 10.1098/rspb.2017.0582.PMC 557747328835551 · doi ↗ · pubmed ↗