Wildfires disturb the natural skin microbiota of terrestrial salamanders

Environmental change can disturb natural associations between wildlife and microbial symbionts, in many cases to the detriment of host health. We used a North American terrestrial salamander system to assess how the skin microbiota of amphibians responds to wildfires. In northern California ’ s red-wood/oak forests, we assessed how recent wildfires affected the skin microbiota of three different salamander species ( Taricha sp., Batrachoseps attenuatus , and Ensatina eschscholtzii ) over two different sampling seasons in 2018 and 2021. We found species-specific responses to wildfire disturbance on the alpha diversity of the skin microbiota of terrestrial salamanders, although burning in general altered the composition of the skin microbiota. The effect of burning on alpha diversities and body condition indices varied by sampling season, suggesting an additional effect of annual climatic conditions on body condition and skin microbiota response. We tested all salamanders for Batrachochytrium dendrobatidis and found four infected individuals in 2018 and none in 2021. Our study documents correlations in the skin microbiota response to an increasing source of disturbance in western North American ecosystems. In addition, our results highlight the need to consider the effects of increased wildfire regimes/intensities and longitudinal effects on wildlife-associated microbiota and animal health.


INTRODUCTION
The future of many wildlife populations is uncertain due to looming threats from climate change, pollution, and habitat loss.With thousands of plant and animal species at risk of extinction, the need for effective conservation strategies is more urgent than ever (Funk et al., 2019;Trevelline et al., 2019).Developing such a strategy requires an understanding of cryptic factors that affect the health of endangered and threatened species.The microbial communities that live inside and outside animals are a prime example of interactions that play an essential role in the health of the host (i.e., microbiota; Barron & Young, 2022).Microorganisms contribute to the fitness and adaptability of animals by regulating processes related to digestion, reproduction, resistance to infection, etc. (Cullen et al., 2020;Hutchins et al., 2019).Consequently, microbiota disturbances can adversely affect overall population health through increased infectious disease prevalence or other illnesses (Redford et al., 2012).Therefore, the stability and diversity of microbiota living on anything from animal intestines to skin could help inform and support conservation efforts (Cullen et al., 2020).Because the composition of an animal's microbiota is closely related to its environment and lifestyle, studying these microbial communities can provide a more holistic insight into the response of wildlife to disturbances such as habitat loss, temperature change, pollution, and chemical contamination.
The impact of increased wildfire regimes and intensities on environmental free-living and host-associated microbial communities remains largely unexplored.The rise in global temperatures associated with climate change is expected to influence hydrological and temperature patterns across the world (Prudhomme et al., 2014).As forested areas become drier and hotter, the frequency and intensity of wildfires are expected to increase (Loehman et al., 2020).In the last few years, the world has witnessed extreme forest fire seasons across the Americas, Europe and Australia (Duane et al., 2021).Wildfires, in addition, can disturb local hydrological patterns through increased erosion and sediment loads in streams (Cole et al., 2020).Changes in land cover and hydrology have the potential to negatively impact both terrestrial and aquatic habitats, alike.Furthermore, interactions between wildfires and other threats, such as infectious diseases or pollutants, among already vulnerable aquatic/semiaquatic wildlife populations are probable (Hossack et al., 2013).Current climatic models predict an increase in fire intensities in these wildfire-prone regions (Goss et al., 2020), which suggests a possible synergistic effect of increased temperatures, decreased moisture and wildfires on these ecosystems.As both temperature, hydrology, and wildfires regimes continue to shift, there is an urgent need to fully understand their combined influence on sensitive wildlife species and the assembly of host-associated microbiota.
Identifying wildfire-impacted host species and determining whether urgent management action is necessary might be key to wildlife management success in response to increased wildfire frequencies (Todd et al., 2016).In the case of infectious diseases, increased fire regimes may alter the environment in ways that either benefit the host (Hossack et al., 2013), transmission vectors (MacDonald et al., 2018) or pathogens (Beh et al., 2012).These changes in host, vector and microbial ecologies merit the need to investigate how wildfires may predispose vulnerable wildlife to infectious diseases.Amphibians serve as excellent models to evaluate hostpathogen responses to environmental change.The amphibian skin harbours diverse microbial communities that provide protection against lethal amphibian pathogens (e.g., Batrachochytrium dendrobatidis -Bd; Harris et al., 2006;Bataille et al., 2016;Jiménez & Sommer, 2017).Because the skin microbiota of amphibians recruits environmental microbes (Walke et al., 2014), environment disturbances can result in notable alterations to the amphibian skin community structure (Hern andez-G omez et al., 2020b;Hern andez-G omez et al., 2020a;Jiménez et al., 2020;Buttimer et al., 2021).Characterizing environmental contributors to variation in the amphibian-associated microbial communities is crucial, as changes in the microbiome (i.e., dysbiosis) can interfere with the defensive phenotype of the skin microbiota resulting in the proliferation of sometimes lethal fungal diseases (Jiménez & Sommer, 2017).
Of particular interest as well is the effect wildfires may have on the prevalence of infectious diseases across the landscape.In some wildfire affected ecosystems, such as the Pacific Coast of the United States, the wildfire season occurs when some amphibians (e.g., terrestrial salamanders) are estivating underground or under water sparing them from the direct effects of burning.Despite this asynchrony, canopy burning can induce moisture/temperature changes in the soil and impact the suitability of microhabitats for both salamanders and Bd, who prefer cool and moist environments (Hossack et al., 2013).Burning may also reduce available microhabitats (i.e., fallen logs/leaves), increasing the density of hosts in remaining habitable environments and increasing opportunities for pathogen transmission among individuals (Albery et al., 2021).Anticipating Bd's response to wildfires in forest floors can be difficult because this pathogen mainly spreads through water, although transmission among terrestrial individuals is possible (Kolby et al., 2015).Furthermore, quantifying the effects of wildfires on a moisturedependent pathogen in terrestrial systems is important, given that direct developing amphibians with 100% terrestrial life histories can be more susceptible to chytridiomycosis (Mesquita et al., 2017).
We evaluated population health, prevalence of Bd, and microbiota structure of wildfire-affected and unaffected populations in the San Francisco Bay Area, California.Several wildfires occurred in the San Francisco Bay Area between 2015 and 2020 (Figure 1).We managed to collect salamanders in areas that experienced wildfires in 2015, 2017, 2019, and 2020 (Figure 1).Our sampling efforts occurred during December 2018-March 2019 (2018 sampling season) and January-February 2021 (2021 sampling season), allowing us to use the forementioned wildfire events as natural experiments to evaluate the response of these populations and their microbial symbionts to wildfire burning.We chose to focus our investigations on terrestrial salamanders as the local terrestrial amphibian fauna contains both direct (Ensatina eschscholtzii and Batrachoseps attenuatus) and larval developers (Taricha torosa, T. granulosa, and T. rivularis), with the former laying eggs that hatch into miniature adults and the latter laying eggs in water bodies where they hatch into aquatic larvae.In the case of Taricha sp., adults can inhabit both aquatic and terrestrial habitats after metamorphosis.We anticipated wildfires to differ in their effect on population health based on species traits.Because of their reliance on wildfire ground fuel for microhabitat, we expected body condition, Bd infection prevalence, and microbiota (i.e., bacterial communities) of direct developers to respond more strongly to wildfire events than those of salamanders with more amphibious life cycles.

Field methods
Several of our salamander sampling locations experienced multiple fires in the past five years.For example, parts of Pepperwood Preserve were burnt during the late summer/fall of 2017 and 2019.Pepperwood Preserve is dominated by savannah and oak habitat, and is located east of the California Coastal Range.Another savannah/oak dominant site located east of the California Coastal Range is the UC Stebbins Cold Canyon Preserve, which experienced major wildfires in 2015 and 2020.Our two most western sampling sites occurred within Point Reyes National Seashore.Habitat on these sites included a mix of the previously described savannah/oak interspersed with patches of mixed redwood/deciduous forest habitat.The northern site at Point Reyes (Bear Valley Trail) experienced a wildfire in 2020, and we could not find any recent burn records for the southern sampling site in the national park (Five Brooks Ranch).To ensure the safety of field researchers, our sampling was restricted to transects along the edges of the fire perimeters where trees were more structurally sound.Across every site, the habitat sampled in burnt areas possessed heterogenous patches of surviving forest floor possessing ashes, burnt logs, and live herbaceous vegetation given that sampling occurred 2-3 months after the last fire and during the rainy season.
We collected salamanders using log flip surveys along transects that begin within $100 m of each site's location coordinates (Table 1).We usually collected salamanders for 1 h or until we obtained more than 20 individuals for a species to limit the time salamanders are handled.Field researchers donned new unsterile gloves between the handling of each salamander.Upon capture, we stored each individual salamander in a new unsterile plastic bag.We marked each log where salamanders were collected from with an identification number in order to return individuals back to their original logs after microbiota sampling.Before collecting skin swabs, we rinsed every salamander with 250 mL of sterile water in order to remove any transient bacteria (Culp et al., 2007).Following rinsing, we swabbed the dorsum of each salamander 30 times then stored the swab on dry ice in a 1.5-mL sterile microcentrifuge tube.We restricted swabbing to the dorsum in order to match the field collection methods of Buttimer et al. (2021) whose data we are also using later in the study.We then measured each salamander's tail and snout-vent length to the nearest millimetre using callipers and mass to the nearest tenth of a gram on a 500-g digital scale.We processed and returned all F I G U R E 1 Salamander sampling sites and forest fire historical perimeters across the San Francisco Bay Area, California.The site marked with a circle did not have a history of wildfires or evidence of recent wildfires on-site (e.g., charred bark, ashes, etc.).Sites marked with stars have both a history of wildfires and evidence of recent wildfires on-site.Wildfire and prescribed fire perimeters downloaded from CalFire are presented, and wildfire season perimeters relevant to this study are coloured separately from all historical fires.salamanders to their marked logs within 1 h of capture.Skin swabs were transferred to a À 20 C freezer within 24 h upon transfer to the laboratory.

Laboratory methods
We isolated DNA from skin swabs collected partially during the 2018 season and during the 2021 sampling season using the DNeasy DNA Isolation Kit (Qiagen N.V., Hilden, Germany), then diluted each sample at a factor of 1:10 using sterile water.We amplified the 16S rRNA v4 locus using PCR with the forward and reverse primers F515/R806.Each reaction consisted of 5.0 μL of template DNA, 7.5 μL of HotStart Master Mix (New England Biolabs, Ipswich, MA), 1.0 μL of 10 nM forward and reverse primers (combined), and 1.5 μL of sterile water for a total of 15 μL per reaction.Our first PCR conditions consisted of 94 C for 3 min, 30 cycles of 94 C for 45 s, 50 C for 60 s, and 72 C for 90 s, followed by 72 C for 10 min.Our microbiome analysis pipeline included PCR/DNA dilution controls to evaluate for contamination.We pooled amplicon triplicates and cleaned the products using the GeneJet PCR Clean-up Kit (ThermoFisher Scientific, Waltman, MA).We performed the second PCR to attach Illumina adaptors and dual index barcodes.The PCR consisted of 5.0 μL of clean amplicons, 7.5 μL 2X MyTaq Master Mix, 1.0 μL of 10 nM forward and reverse barcode primers (Hern andez-G omez et al., 2017), and 1.5 μL of water for a total of 15 μL reactions.Our second PCR conditions consisted of 94 C for 3 min, 5 cycles of 94 C for 45 s, 65 C for 60 s, and 72 C for 90 s, followed by 72 C for 10 min.We quantified the PCR products using a Qubit Fluorometer HS Kit (ThermoFisher Scientific, Waltman, MA), pooled them in equimolar amounts, and cleaned the pool one more time using the GeneJet PCR Clean-up Kit.We shipped the sample pool in dry ice via the United States Postal Service overnight to the Cornell University Institute for Biotechnology for sequencing on an Illumina MiSeq machine (Illumina Inc., San Diego, CA) using a 250 bp pairedend sequencing kit.
To detect Bd infections, we implemented quantitative PCR (qPCR) reactions in triplicate following the protocol of Boyle et al. (2004) with some modifications.We analysed the DNA extracted in 2021 and the DNA samples from Buttimer et al. (2021) that were stored in À80 C for 3 years.Each qPCR reaction was made using 5 uL of 1:10 diluted template DNA, 12.5 uL of 2X TaqMan Fast Advanced Master Mix (ThermoFisher Scientific, Waltman, MA), 1.5 uL of 10 nM forward and reverse primers, 1 uL of minor groove binder probe MGB2 (Boyle et al., 2004), and 3.5 uL of moleculargrade water for a total volume of 25 uL.We included 3 replicates of the amphibian disease gBlock (Standish et al., 2018) in serial dilutions of 1,000,000 to 1 copies/ uL or zoospore equivalents (ZEs).We ran and analysed these reactions on an Applied Biosystems Ste-pOnePlus RealTime PCR (ThermoFisher Scientific, Waltman, MA) system using the thermoprofile program outlined in Kerby et al. (2013).We classified as negative samples those whose average qPCR quantification fell below 1 ZE after correcting for dilution effects.1).Several materials were discontinued between the time the laboratory work for Buttimer et al. (Buttimer et al., 2021) and the current study were performed, so there were some differences in the brand/manufacturer of the kits we used to prepare both DNA sequencing libraries.However, we made sure to use products (e.g., Qiagen Dneasy Kit) that have been identified to produce similar 16S rRNA results (Mattei et al., 2019;Velasquez-Mejia et al., 2018) to the materials used in Buttimer et al. (2021;e.g., Qiagen PowerSoil Kit).We combined the raw DNA sequence files for both sequencing runs, then quality filtered raw reads using Trimmomatic (Bolger et al., 2014) to remove adapter sequences, bases below phred-20 from both ends of reads, and reads under 30 bp length.We paired reads that passed initial quality control using PANDAseq (Masella et al., 2012), and only paired read files advanced to subsequent analyses.
We analyzed the paired reads using the Quantitative Insights Into Microbial Ecology version 2.2020.2(QIIME2) pipeline (Bolyen et al., 2019) to characterize the microbiota.We processed reads using the DADA2 plugin to filter, dereplicate, remove chimeras, and denoise reads using the default settings then remove reads below/trim reads below 240 bp length (Callahan et al., 2016).We chose to also remove rare and private ASVs (i.e., those found in only one individual) according to Bokulich et al. (2013).We aligned ASV representative sequences using MAFFT and generated a phylogenetic tree in FastTree2 for later use in R (Katoh & Standley, 2013;Price et al., 2010).We applied a pre-trained Naïve Bayes classifier on the Greengenes 13_8 database to assign taxonomy to each ASV at the genus level (DeSantis et al., 2006).In addition, we filtered out any ASV's whose taxonomy matched chloroplast or mitochondria as these were not the target of our amplification protocol.To standardize sequencing depth throughout all samples, we rarefied the filtered ASV table to 2103 sequences per sample.

Statistical analysis
To quantify how wildfires influence terrestrial salamanders across the San Francisco Bay Area, we first grouped our samples into categorical variables based on their habitat's documented wildfire history (i.e., burnt or unburnt) and sampling season year (i.e., 2018 or 2021): salamanders collected prior to the summer of 2019 from sites with unknown burnt history and no visible signs of a recent wildfire (sampling year 2018; unburnt habitat), salamanders from sites with documented wildfire histories collected prior to the summer of 2019 (sampling year 2018; burnt habitat), salamanders from sites with zero documented wildfire at the end of the study (sampling year 2021; unburnt habitat), and salamanders from sites that experienced a wildfire event during the span of the study (sampling year 2021; burnt habitat).
We implemented linear mixed models to assess the effects of wildfire status; salamander genus identity; sampling year; the interaction between wildfire status and sampling season year; and the interaction between wildfire status and salamander genus identity on alpha diversity metrics (log-transformed community richness and Shannon Diversity Indices).We used salamander genus identity as an independent variable instead of species identity in these analyses and beyond, because we did not find the same Taricha species across all our sites and in some cases sympatric T. torosa and T. granulosa are indistinguishable.To account for repeated sampling across years, we included location as a random variable.We then used a backwards selection approach to remove variables and test models using AICc criterion (package lmerTest Kuznetsova et al., 2020).We applied Tukey HSD posthoc tests to evaluate alpha diversity variation within the categories of significant variables.We also assessed the effect of genus identity, wildfire status, sampling season year and location on unweighted UniFrac, and Bray-Curtis distances using PERMANOVA tests under the package vegan (Oksanen et al., 2022).We used both phylogenetic-based (unweighted UniFrac) and abundance-based (Bray Curtis) beta diversity metrics in order to cover different aspects of microbiota compositional change among samples.We generated NMDS plots using the packages phyloseq (McMurdie & Holmes, 2013) and ggplot2 (Wickham, 2016) in order to visualize clustering using the unweighted UniFrac and Bray-Curtis distance matrices.We repeated these analyses within Batrachoseps and Ensatina independently in order to evaluate intra-species effects of burning.We excluded Taricha from intra-species analyses given their limited sample size across seasons (Table 1).To assess the effect of burning on alpha diversity within each salamander species, wildfire status was entered as a dependent variable and location as a random variable into linear mixed models.We applied Tukey HSD posthoc tests to evaluate alpha diversity variation within the categories of significant variables.We also used PER-MANOVA tests to quantify and test the effect of wildfire status on salamander skin microbiota composition (unweighted UniFrac and Bray-Curtis).We used the package indicspecies (De C aceres & Legendre, 2009) to identify the ASVs that were significantly more abundant among different salamander genera, and between burnt and unburnt sites.
We found significant differences in baseline body condition indices across the three genera ( p < .001);therefore, we evaluated body condition indices separately among Batrachoseps and Ensatina only and included data for salamanders with and without microbiome data (Table 1).We again excluded Taricha from intraspecies analyses due to their reduced sample size.We used linear models to assess wildfire status, sampling season year, and their interaction as independent fixed variables, location as a random variable, and body condition indices as the dependent variables.We applied Tukey HSD post-hoc tests to evaluate body condition index variation within the categories of significant variables.

RESULTS
We collected 36 Batrachoseps attenuatus, 27 Ensatina eschscholtzii, and 23 Taricha sp.across all four sites during the 2018 sampling season (Table 1).Following the three wildfire events at Bear Valley Trail -Point Reyes National Seashore, Pepperwood Preserve, and UC Stebbins Cold Canyon Reserve during the 2019/2020 wildfire seasons, we collected 53 Batrachoseps attenuatus, 15 Ensatina eschscholtzii, and 3 Taricha sp.individuals across all four sites in January/ February of 2021 (Table 1).After DNA sequencing, DNA sequence quality filtering, and PCR control sample removal, the final 16S rRNA amplicon dataset contained pre-and post-wildfire season data for 60 Batrachoseps sp., 16 Ensatina sp., and 21 Taricha sp.Only four salamanders tested positive for Bd across our datasets, and all these positive samples originated from the 2018 sampling season (Table 1).Given the low Bd prevalence, we chose not to perform any statistical analyses using the Bd results.
Model selection for the linear mixed model evaluating the effect of wildfire status, sampling season year, and their interaction on skin microbiota richness among Batrachoseps excluded the interaction effect only (wildfire status: F 1,8.60 = 6.44, p = 0.033; sampling year: F 1,47.18 = 6.62, p = 0.013).Batrachoseps skin microbiota samples collected in 2018 possessed significantly higher richness compared to those collected in 2021 (t = 2.57, df = 47.20,p = 0.0133).Regardless of the effect of sampling year, Batrachoseps skin microbiota richness was significantly higher in burnt areas compared to unburnt areas (t = 2.54, df = 8.6, p = 0.033).A different pattern emerged among Shannon Diversity Indices, where wildfire status was the only retained independent variable during linear mixed model selection (F 1,11.84 = 7.30, p = 0.019).In the case of Shannon Diversity, Batrachoseps collected from burnt sites possessed on average more diverse skin microbial communities than those sampled from unaffected sites (Figure 2A; t = 2.70, df = 11.80,p = 0.019).Model selection for the linear mixed model evaluating the effect of wildfire status, sampling season year, and their interaction on skin microbiota richness and Shannon Indices among Ensatina salamanders excluded all variables (Figures S1A and 2A).
Nine ASVs were associated with Batrachoseps, 119 with Ensatina, and 270 with Taricha salamander skin microbiota.Most of the ASVs associated with each salamander genera were rare (relative abundance <1%); however, we did find a few common ASVs assigned to each salamander genus.Bacteria strongly (effect size >0.6)associated and abundant on Batrachoseps skin included ASVs assigned to the family Alcaligenaceae, genus Halomonas, class Verrucomicobiae, and family Microbacteriaceae.ASVs abundant in Ensatina include several assigned to the family Xanthobacteraceae, the family Alcaliganaceae, and the genus Luteibacter.ASVs strongly associated with Taricha were more numerous and encompassed members of the phyla Actinobacteria, Proteobacteria and Bacteroidetes.Among these ASVs, those associated with the family Micrococcaceae, the genus Sphingomonas, and the genus Pseudomonas were abundant on the skin of newts.
The indicator species analysis performed to evaluate which ASVs differed among the four wildfire status categories (2018 unburnt, 2018 burnt, 2021 unburnt, and 2021 burnt) assigned 347 ASVs to one of the categories (Figure 3).Fifty-two ASVs were associated with 2018 unburnt sites, with the strongest associations belonging to ASVs identified to the genus Aeromonas, the family Alcaligenaceae, and the genus Pseudomonas.Two-hundred and seven ASVs were assigned to 2018 burnt sites, with the most abundant ASVs in this category belonging to the genus Pseudomonas.Fifty-one ASVs were associated with 2021 unburnt sites, with the strongest associated ASVs belonging to the genera Pseudomonas and Sphingomonas.Finally, 37 ASVs were assigned to 2021 burnt sites; however, most ASVs assigned to this category possessed effect sizes <0.6.Of these ASVs, many belonged to the genus Petrimonas and the order Bacteroidales.Burning had a negative effect on Alcaligenaceae and Microbacteriaceae associated with Batrachoseps, and a positive effect on the relative abundance of Micrococcaceae, Sphingomonas, and some Pseudomonas ASVs assigned to the skin of Taricha sp (Figure 3).However, we also found a negative effect on some Pseudomonas ASVs assigned to Taricha sp. as well.We did not observe any crossover effects between Ensatina-associated ASVs and wildfire status.
Batrachoseps collected in 2021 were significantly larger than those collected in 2018 (t = À4.74,df = 92.30,p < 0.001; Figure 4) and the linear model identified sampling season as the only significant variable associated with body condition indices (F 1,92.27= 22.51, p < 0.001).The largest Batrachoseps were collected in 2021 from areas burnt by wildfires, followed by those collected from unburnt areas in 2021, those collected from unburnt areas in 2018, then those collected from burnt areas in 2018 (Figure 4).We found a similar pattern among Ensatina, where sampling season was the only retained variable during linear mixed model selection (F 1,34 = 10.39,p = 0.003).As with Batrachoseps, Ensatina collected in 2021 were larger than those collected in 2018 (t = À3.22,df = 34.00,p = 0.003).On average, the largest Ensatina were collected from burnt areas in 2021, followed by those collected from unburnt areas in 2018/2021, and burnt areas in 2018 (Figure 4).

(A) (B)
F I G U R E 2 Alpha and beta diversities among salamander species and wildfire burnt histories.Alpha diversity measure presented is Shannon Diversity Index (A).Alpha diversity plots include mean and standard error bars for each salamander genus (bold) and each wildfire status category (thin).Tukey HSD Post-hoc multiple comparison groupings are presented in upper-case letters to differentiate groupings by species and in lower-case letters to differentiate groupings within wildfire burnt status/year (except Taricha due to their limited sample size).Beta diversity NMDS plots are presented using Bray-Curtis (B) diversity metrics.

DISCUSSION
Our study characterized the skin microbiota, Bd prevalence, and body condition indices of terrestrial salamanders in areas affected by wildfires across the San Francisco Bay Area of California, USA.We collected three different genera of salamanders before and after the 2020 wildfire season from both affected F I G U R E 3 ASVs (rows) associated with different salamander species and wildfire status.ASV taxonomies are presented.Assignment columns denote the category (salamander species or wildfire status) that each ASV was associated with based on indicator species analyses.Heatmaps denoting the average relative abundance of that ASV across each salamander/wildfire history category are also presented.and unaffected habitats.We only detected Bd in four of the 2018 samples, and the prevalence of the fungal pathogen was below levels indicative of active chytridiomycosis (Carter et al., 2020).In Batrachoseps attenuatus and Ensatina eschscholtzii, we recovered the largest body condition indices from salamanders collected in 2021 compared to 2018.The plethodontids (Batrachoseps attenuatus and Ensatina eschscholtzii) we collected from burnt sites had larger body condition indices than those from unburnt sites.Skin microbial community richness varied among species, and we were able to see host species specific responses to wildfires in both plethodontids.We also observed an interaction of salamander genus identity and wildfire status on beta diversity metrics, with the strongest effect of wildfire status observed on the presence/ absence-based beta diversity metrics (i.e., Unweighted UniFrac).We identified numerous ASVs associated with each genus and four wildfire status categories, and among these only a few ASVs possessed relative abundances above 0.1%.We observed lower relative abundances of bacterial taxa commonly found on Batrachoseps attenuatus and higher relative abundances on bacteria associated with Ensatina eschscholtzii and Taricha sp in wildfire-burnt areas.These results provide evidence for the impact of wildfires on terrestrial salamander skin microbiota and host populations.The high prevalence of significant interaction terms in our models suggests that responses of skin microbiota to wildfire disturbances may vary based on the host species' natural history and skin morphology.
In our observations, Taricha sp. had the richest skin microbial communities across time, followed by Ensatina eschscholtzii and Batarachoseps attenuatus.When we compared Shannon diversity, Ensatina eschscholtzii possessed more diverse communities followed by Taricha sp. and Batrachoseps attenuatus.We also saw some clustering by host genus in our microbiome compositional analyses, and we were able to identify consistent associations between several ASVs and salamander genera across multiple sampling seasons.We anticipate that salamander microbiota assembly mechanisms in these systems result from a combination of physiological and natural history traits of the host.While we made sure to include only terrestrial forms of Taricha, plethodontids and newts possess differences in skin morphology, behaviour, and immunity that could contribute to variation in the microbial microhabitat of the skin (Jani & Briggs, 2018;Wakely et al., 1966;Walke et al., 2014).In addition, home range and microhabitat use differences among species could account for the patterns of variation observed in the skin microbiota's alpha and beta diversities (Hern andez-G omez & Hua, 2023).Taricha sp. can have yearly migrations between aquatic and terrestrial habitats, with a record movement of up to 3.5 miles in some instances (Darrow, 1966).Ensatina eschscholtzii have been documented to move up to 120 meters as well (Staub et al., 1995).In contrast, Batrachoseps attenuatus have significantly smaller bodies and only move around 3 ft, sometimes not even moving from their shelter for long periods of time (Cunningham, 1960).Given that amphibians can sequester bacteria from the environment into their skin microbiota (Loudon et al., 2014;Loudon et al., 2016;Walke et al., 2014), greater dispersion may increase the interaction with different types of environmental bacterial reservoirs resulting in more diverse and variable communities.Because we observed these patterns across multiple seasons, our results offer longitudinal support for the correlation between host eco-morph and patterns of alpha and beta diversity variation among different members of amphibian communities (Bletz et al., 2017a;Bletz et al., 2017b;Kueneman et al., 2019).
We were able to observe some consistent changes in alpha diversity and composition between plethodontid salamanders inhabiting burnt and unburnt habitat; although, the effect of burning sometimes varied by sampling season.We anticipate that the contradictory response of plethodontid microbiota to wildfires across years could result from variation in time since the last wildfire.Our 2018 burnt samples were collected more than 2 years after the last wildfire, while those in 2021 originated months after the last fire event.While our study did not directly measure soil microbiota, our observations correlate with those performed by microbiota surveys of burnt soils.In Villadas et al. (2019), Dove et al. (2022), andPulido-Chavez et al. (2023), a swift reduction in diversity is typically observed immediately after a fire followed by a boost in taxonomic diversity that results from the release of pyrogenic nutrients.In our observations, salamander skin microbiota possessed a greater taxonomic diversity of bacterial phyla in burnt areas immediately after a wildfire (those sampled from burnt areas in 2021) compared to those sampled a few years after the last wildfire (those sampled from burnt areas in 2018).Therefore, our results suggest that salamander microbiota could parallel the ecological succession in soil microbial communities that have been described following a burning event.
Whether burning influences salamander skin microbiota directly or indirectly through colonization from altered environmental microbial communities remains to be explored.To fully understand the effects of wildfires on amphibian microbiota, laboratory-and mesocosm-based studies are needed to appropriately identify the mechanisms behind the changes we observed in our natural experiment.
Among Batrachoseps attenuatus, we observed a consistent negative response of a common ASV identified to the family Alcaligenaceae to burning across sampling seasons.Members of the family Alcaligenaceae have been commonly described on the skin of terrestrial amphibians across the world (Abarca et al., 2018;Belden et al., 2015;Bletz et al., 2017c).Our study is no exception, as we were able to find an abundant ASV categorized as Alcaligenaceae the 100% terrestrial Batrachoseps attenuatus/Ensatina eschscholtzii salamanders.In more than one system, a higher proportion of Alcaligenaceae in the skin microbiota of amphibians has been associated with lower incidence of Bd (Jiménez et al., 2019;Le Sage et al., 2021).Because an increase in diversity following wildfires may negatively influence the relative abundance of common and potentially important skin symbionts like Alcaligenaceae, longitudinal investigations that investigate how this taxonomic rearrangement influences important skin microbiota functions in amphibians such as innate immunity are needed.

Terrestrial salamander conservation implications
We did not find Bd positive individuals during the 2021 season.We also observed greater alpha diversities and body condition indices in the later sampling season, and this pattern was consistent between Batrachoseps attenuatus and Ensatina eschscholtzii.Among these plethodontid salamanders, we also saw greater body condition indices in burnt sites compared to their unburnt sites.2020 was a relatively dry year for California (Robinne et al., 2021), and it is possible that the absence of Bd in later samples and larger body conditions/alpha diversities found in 2021 are in response to climatic variation between both seasons.We cannot say the increase in skin microbial diversity/ absence of pathogen infections is related to body conditions as we used a natural experiment and establishing these links would require further experimentation.Still, none or a low Bd infection prevalence has been documented among Batrachoseps, Ensatina, and Taricha in California in previous studies (Bird et al., 2018;Cowgill et al., 2021;Sette et al., 2020), and in our case it is encouraging to observe this same pattern in salamanders collected from unaffected and wildfire burnt areas.
The shift towards larger salamanders following wildfire events is concerning, as it is consistent with previously observed recruitment failures/loss in other plethodontids (e.g., Plethodon shermani) within 6 months of wildfires (Gade et al., 2019).In our case, it is possible that larger salamanders are at a greater advantage at dealing with changes in prey communities, moisture, and temperature following a wildfire explaining their greater frequency.However, long-term studies ($2 years post wildfire) have observed evidence of recruitment following a wildfire event suggesting that terrestrial salamander populations can recover (Cummer & Painter, 2007).Still, in the face of increasing fire intensities and frequencies conservation of terrestrial salamanders in western North America and beyond may benefit from continuous longitudinal research that assesses the effect of burning on chytridiomycosis, host health, and recruitment.Considering how different fire intensities influence the prevalence of diseases and the structure of microbiota in vertebrate hosts remains mostly unexplored as well.In our study we used the edges of wildfire burnt areas, which in recent history tend to be more intense and frequent in the Pacific Coastal Range landscape.Thus, future research into the effect of burning on vertebrate host populations, their pathogens, and microbiota might also benefit from including populations living in areas that have experienced high intensity (i.e., mega wildfires) to low intensity fires (i.e., prescribed fires).
In recent years, microbiota studies have unlocked a wealth of data that provide insight into the nature of host-associated microbial communities (Cullen et al., 2020).With the advances in technology and dramatic cost reductions it has become much easier to study microbiota interactions and influences both within the host and with the external environment (Cullen et al., 2020).Because microbiomes (i.e., the collection of microbial cells and their genes/functionality) are so important to host health, characterizing changes in response to stressors continues to be an important priority.Therefore, identifying factors that are likely to alter natural associations between hosts and microbiomes will enable the development of intervention strategies for wildlife hosts affected by environmental change (Cullen et al., 2020).
T A B L E 1 Salamander sampling frequencies for San Francisco Bay Area, California sites affected by wildfires.Sample sizes for both 2018 and 2021 sampling seasons used in this study are presented separately.Sampling effort (i.e., number of surveyors and survey time) differed between the 2018 and 2021 sampling seasons, as such this table is meant to represent species site occupancy only.Values in brackets represent microbiota 16S rRNA sequencing sample sizes.Values in parenthesis represent individuals who tested positive for Batrachochytrium dendrobatidis.
Salamander body condition indices in wildfire affected/unaffected sites.Data for Batrachoseps attenuatus, Ensatina eschscholtzii, and Taricha sp. are presented and calculated independently.Body condition index represents the residuals for the linear relationship between salamander mass and total body length.Wildfire status (unburnt/burnt) of habitat sampled each season (2018/2021) is noted within the plot.Tukey HSD post-hoc multiple comparison groupings are presented in upper-case letters to differentiate groupings by wildfire burnt status/year for each species except Taricha due to their limited sample size.Means and standard error bars are presented for each wildfire category and salamander combination.