Representing one of the five major clades of terrestrial vertebrates on the planet, amphibians (anurans, salamanders, and caecilians) display an astonishing diversity of life history traits (Vitt and Caldwell, 2013; Crump, 2015; Bardua et al., 2021). Their reproductive diversity (e.g., biphasic tadpole metamorphosis, direct development, etc.) and specialized ecomorphological adaptations (e.g., webbed toes, enlarged toe pads, etc.) have enabled amphibians to diversify across a near complete spectrum of environments on the planet, including fossorial, terrestrial, arboreal, aquatic, and riparian habitats (Afonso and Eterovick, 2007; Gómez and Lires, 2019; Engelkes et al., 2020; Bardua et al., 2021). To persist in such an expansive array of environments, amphibians have evolved specialized skin that serves as an essential respiratory organ and performs numerous physiological functions, including ion transport and osmoregulation, while also acting as a defensive barrier against harmful environmental and pathogenic agents (Duellman and Trueb, 1994; Lillywhite, 2006; Campbell et al., 2012; Kueneman et al., 2014; Wake and Koo, 2018; Varga et al., 2019). However, this same unique morphological feature also makes amphibians highly sensitive to instability in local environments, such as habitat modification and degradation, which has led to their recognition as bioindicators of ecosystem health (Wake, 1991; Blaustein et al., 1994; Houlahan et al., 2000; Ethier et al., 2021) and has contributed to alarming patterns of population declines on a global scale (Ford et al., 2020).

Currently, a staggering 41% of IUCN-assessed amphibian species are threatened with extinction (IUCN, 2023), with global population declines linked to a variety of factors, from pollution and invasive species to climate change and human-mediated habitat modification and destruction (Allentoft and O’Brien, 2010; Campbell Grant et al., 2020; Ford et al., 2020; Green et al., 2020). Furthermore, global declines have been exacerbated by growing threats posed by amphibian infectious diseases, including fungal and viral pathogens (Carey et al., 1999; Daszak et al., 1999; Greenberg and Palen, 2019; Scheele et al., 2019). In particular, Batrachochytrium dendrobatidis (Bd), a fungal pathogen that causes chytridiomycosis (often referenced as chytrid), has been especially devastating to anurans (frogs and toads), the most diverse of the three extant orders of amphibians with more than 7,500 species making up 88% of all amphibian diversity (AmphibiaWeb, 2023). To date, Bd infection has contributed to large-scale population declines of numerous species, as well as the confirmed or presumed extinction of more than 80 lineages to date (Scheele et al., 2019). Given that frogs and toads are essential members of both aquatic and terrestrial ecosystems (Whiles et al., 2006; Hocking and Babbitt, 2014), an improved understanding of the factors impacting anuran health, persistence, and immunity is needed for effective identification and mitigation of threats posed by rapid environmental change, habitat modification, and emerging pathogens.

Recently, a growing focus on the symbiotic microbial communities of host organisms, referred to as microbiomes, has revealed their fundamental roles in numerous processes associated with vertebrate host health, including digestion, nutrient acquisition, metabolism, immunity, development, and behavior (Ley et al., 2008; McFall-Ngai et al., 2013; McFall-Ngai, 2014; Alberdi et al., 2016; Colston and Jackson, 2016; Arizza et al., 2019; Rollins-Smith, 2020; Sehnal et al., 2021). Certain gut microbiota been shown to improve host metabolic activity in response to harsh environmental conditions (Li et al., 2019), with others often altering their composition and gene-expression patterns in response to physiological changes imposed by the host organism, underscoring the important role microbiomes may play in promoting host adaptation (Alberdi et al., 2016). Anuran skin microbiomes also serve a critical role as a protective barrier against pathogens (Harris et al., 2006; Harris et al., 2009; Becker et al., 2011; Kueneman et al., 2016; Ross et al., 2019; Rebollar et al., 2020; Jani et al., 2021). For example, metabolites produced by members of the frog skin microbiome have been shown to inhibit Bd zoospore development (Brucker et al., 2008; Loudon et al., 2014a; Walke and Belden, 2016). Furthermore, previous research has found that the presence of anti-fungal microbiota, such as Janthinobacterium lividum, on the skin of uninfected frogs reduced morbidity and mortality when the organisms were exposed to Bd (Harris et al., 2009; Walke and Belden, 2016). Unfortunately, despite the recognition that skin microbiome composition is influenced by many host-associated and environmental factors (Kueneman et al., 2014; Bletz et al., 2017; Carda-Diéguez et al., 2017; Prado-Irwin et al., 2017; Chiarello et al., 2018; Sehnal et al., 2021), little remains known about the patterns and the degree with which skin microbial communities change in response to environmental variation, particularly among diverse anuran hosts.

Studies to date indicate that anuran skin microbiomes are likely shaped by a number of factors such as host evolutionary history (McKenzie et al., 2012; Kueneman et al., 2014; Walke et al., 2014), host ecology (Bletz et al., 2017), and developmental stage (Kueneman et al., 2014; Jiménez et al., 2019). Additionally, local environment (Varela et al., 2018; Ellison et al., 2019; Rebollar and Harris, 2019), geographic location (Kueneman et al., 2014; Belden et al., 2015; Rebollar et al., 2016; Bletz et al., 2017), and the prevalence of emerging infectious disease (i.e. Bd infection; Rebollar et al., 2016; Familiar-López et al., 2017; Bell et al., 2018; Knutie et al., 2018; Rebollar and Harris, 2019; Jani et al., 2021) have been shown to influence anuran skin microbiome diversity. Although previous studies have compared the skin microbiomes of host communities from different sites (Kueneman et al., 2014; Belden et al., 2015; Rebollar et al., 2016; Bletz et al., 2017; Medina et al., 2017), few have systematically sampled frog skin microbiomes across environmental gradients (Bletz et al., 2017; Medina et al., 2017). Yet, climatic variation across ecotones is known to impact anuran community assembly through processes such as environmental filtering (Cortés-Gómez et al., 2013; Díaz-García et al., 2017; Álvarez-Grzybowska et al., 2020), implying that there are major gaps in our understanding of how large-scale climatic and ecological variation may impact anuran skin microbiomes, and in turn, host survival and persistence. Additionally, many anuran skin microbiome studies compare only a few (i.e., 1–3) host species (McKenzie et al., 2012; Kueneman et al., 2014; Walke et al., 2014; Varela et al., 2018). As such, there exists a need for studies that investigate communities of co-distributed anuran species representing a variety of host ecologies across complex environmental landscapes to better understand, and predict, expected future shifts in anuran skin microbiomes as climatic and habitat changes progress.

In this study, we evaluate the skin microbial communities of five widely distributed species of frogs across four major Oklahoma ecoregions to elucidate the roles that geographic location, host evolutionary history, host ecology, and pathogen presence play in anuran skin microbiome assembly. Ecoregions are areas of general ecosystem similarity characterized by abiotic and biotic variables such as climate, geology, soil, and vegetation (Omernik, 1995; Bryce et al., 1999; McMahon et al., 2001). Oklahoma is one of only four states in the United States that possesses more than 10 distinct ecoregions (U.S. Environmental Protection Agency, 2013; Figure 1 ), with 29 native species of frogs distributed across broad regions of this complex landscape (Sievert and Sievert, 2021). Furthermore, recent efforts to determine and monitor the distribution and prevalence of Bd infection among amphibian communities across the state has resulted in robust pathogen datasets that are now publicly available (Marhanka et al., 2017; Watters et al., 2018; Watters et al., 2019; Watters et al., 2021). Therefore, the state represents an ideal spatial framework for investigating the evolutionary and ecological processes involved in anuran skin microbiome assembly dynamics across changing environments (McMahon et al., 2001; Omernik, 2004) and in the presence of growing threats from emerging infectious diseases (Marhanka et al., 2017; Watters et al., 2018; Watters et al., 2019; Watters et al., 2021).

We perform both intra- and interpopulation comparisons of frog skin microbiomes across four ecoregions in Oklahoma for populations of Blanchard’s cricket frog (Acris blanchardi), the American green tree frog (Hyla cinerea), members of the morphologically indistinguishable gray treefrog species complex comprised of Gray’s (H. chrysoscelis) and Cope’s gray (H. versicolor) treefrogs (Jaslow and Vogt, 1977), which we refer to as H. chrysoscelis/versicolor, the American bullfrog (Rana catesbeiana), and the coastal plains leopard frog (R. sphenocephala) ( Figure 1 ; Supplementary Table 1 ). Through these comparisons, our study aims to address the roles that (1) geographic location, (2) host ecology, (3) host evolutionary history, and (4) pathogen presence play in anuran skin microbiome assemblage. If geographic location has a dominant role in shaping frog skin microbial communities, we expect regional populations of each host species to possess distinct skin microbiomes. Conversely, if skin microbiome assemblage is impacted largely by host ecology, we expect frog species with similar life histories (aquatic, semi-aquatic, arboreal) to possess similar microbial community structure and diversity, regardless of their evolutionary history or geographic location. Furthermore, by comparing populations of closely related taxa within two distinct anuran families, we test whether host evolutionary history plays a role in microbiome assemblage by examining whether closely related species harbor greater microbial similarities than can be explained by geographic location or host ecology. Finally, if pathogen presence drives differences in skin microbiome diversity, we predict that anurans infected with Bd will possess distinct skin microbiomes when compared to frogs that tested negative for the fungus. Results from this comparative study contribute to our understanding of the environmental and host-associated processes that drive vertebrate skin microbiome assemblage. Furthermore, we provide the first assessment of variation in skin microbiomes among frog species spanning an ecological gradient in North America, allowing for future studies to investigate the impact environmental stressors, disease, and urbanization will have on population and species persistence.

Among the 179 skin microbiome swabs we collected from our five focal frog species, we obtained 4,050,580 sequences which were clustered into 22,641 OTUs using a 97% sequence similarity threshold against the Silva database (version 138, Quast et al., 2012). After rarefying the 179 samples to a sequencing depth of 1,000 sequences and removing positive and negative control samples, 113 samples and 7,694 OTUs remained and were used in subsequent analyses (A. blanchardi [N = 17], H. chrysoscelis/versicolor [N = 27], H. cinerea [N = 24], R. catesbeiana [N = 30], and R. sphenocephala [N = 15]; Supplementary Table 2 ). Visualization of the data obtained from the two positive control samples indicated that the community compositions of the Extraction-Stage and PCR-Stage ZymoBIOMICS Microbial Community Standards were similar, suggesting that our extraction-to-analysis workflow performed consistently ( Supplementary Figure 2 ).

Our study sampled the skin microbial communities of anuran hosts across distinct ecoregions in Oklahoma to determine if geographic location, host ecology, pathogen infection, or host evolutionary history correlated with differences in frog skin microbiomes. Results from the CART analysis indicated that host ecology and geographic location were the primary drivers of the skin microbiome diversity differences we observed among the five focal frog species sampled ( Supplementary Figures 4–6 ). Although we observed some compositional differences among skin microbiome samples based on the geographic locations of our hosts, such patterns were not shared across all five host species ( Figure 2 ). Interestingly, the most prominent differences were observed among A. blanchardi skin microbiomes, which differed profoundly between individuals from the two unique ecoregions we sampled ( Figure 2 ). Four phyla, Actinobacteria, Bacteroidetes, Cyanobacteria, and Verrucomicrobia, were present in abundances greater than 1% among samples from Arkansas Valley only, with the exception of two samples from the Crosstimbers ecoregion. This is likely because most of the A. blanchardi samples from the Crosstimbers were dominated largely by Proteobacteria ( Figure 2 ). For the geographic location hypothesis to be fully supported, we would expect to see more ecoregion-specific compositional differences among the skin microbiomes of our other focal anuran species. Additionally, Bd infection status and severity did not emerge as primary drivers of anuran skin microbiome diversity among our focal samples ( Supplementary Figures 4–6 ). However, previous studies have found that Bd infection impacts the composition of anuran skin microbiomes (Rebollar et al., 2016; Familiar-López et al., 2017; Rebollar and Harris, 2019; Jani et al., 2021), and disease data was only available for a subset of samples in our study (N = 71 out of 113 samples). Further, while both Hyla species had associated disease data for 8–13 individuals per ecoregion, R. catesbeiana and R. sphenocephala only had disease information for 2–6 individuals per ecoregion, and no R. catesbeiana individuals from the South Central Plains ecoregion had corresponding Bd disease data ( Figure 1 ; Supplementary Table 1 ). As such, it may be difficult to detect a correlation between Bd disease status/severity and microbiome compositional differences when sample sizes are small. Therefore, increased sampling is needed to further investigate how Bd disease status and severity alters anuran skin microbiome diversity. We did not find evidence of phylogenetic signal among the skin samples from our focal hosts, indicating that anuran skin microbiomes are influenced by host ecology and geographic location, whereas host evolutionary history and pathogen infection did not emerge as primary drivers.

Four dominant bacterial phyla emerged within nearly all samples (N = 111/113), Actinobacteria, Bacteroidetes, Firmicutes, and Proteobacteria, a pattern consistent with findings from other amphibian skin microbiome studies to date (Kueneman et al., 2014; Walke et al., 2014; Walke and Belden, 2016; Bletz et al., 2017). Additionally, previous research has indicated that Chloroflexi (Walke et al., 2014; Ellison et al., 2019), Cyanobacteria (McKenzie et al., 2012; Walke et al., 2014; Ellison et al., 2019), and Planctomycetes (Jiménez and Sommer, 2017; Prado-Irwin et al., 2017; Muletz Wolz et al., 2018; Albecker et al., 2019; Ellison et al., 2019) are also dominant members of the amphibian skin microbiome. Therefore, it is not surprising that we found these phyla in abundances greater than 1% among samples from four of our focal host anurans (excluding H. cinerea). Although Firmicutes was a dominant phylum among skin microbiome samples from all of our focal host species, some R. catesbeiana samples had abundances of Firmicutes as high as 93.54% ( Figure 2 ). This was surprising given that previous research, and results from our study, indicate that Proteobacteria dominates the amphibian skin microbiome (Walke et al., 2014; Walke et al., 2015; Bataille et al., 2016; Medina et al., 2017; Muletz Wolz et al., 2018; Figure 2 ). However, studies have also found that skin microbial communities dominated by Firmicutes are associated with the presence of skin lesions caused by psoriasis (Gao et al., 2008; Yan et al., 2017) and atopic dermatitis (Zheng et al., 2019) in humans, and papillomatous digital dermatitis in cows (Santos et al., 2012). Therefore, it is possible that the R. catesbeiana individuals with high abundances of Firmicutes in their microbiomes could have experienced prior skin infections or were infected at the time of sampling. However, only a small subset (N = 6) of the focal R. catesbeiana hosts had associated Bd infection data. As such, we cannot confirm if Bd presence is correlated with an abundance of Firmicutes in the skin microbiomes of this species at this time, but future studies should investigate if amphibian skin lesions are also colonized by Firmicutes, as observed in the mammalian skin microbiome literature (Santos et al., 2012; Zheng et al., 2019).

We found species-specific compositional differences among the five anuran species sampled. Notably, the Archaea phylum Crenarchaeota was found among aquatic and semi-aquatic species only (A. blanchardi, R. catesbeiana, and R. sphenocephala; Figure 2 ). All cultured Crenarchaeota are extreme thermophiles (Ochsenreiter et al., 2003); however, with increased use of culture-independent sequencing methods, studies have indicated that there are members of the Crenarchaeota phylum that are associated with freshwater sediments and other soil environments (Ochsenreiter et al., 2003; Brochier-Armanet et al., 2008). Acris blanchardi, R. catesbeiana, and R. sphenocephala are commonly encountered within or around water and marshy habitats, which could explain why Crenarchaeota was found on the skin of these species and not the two arboreal host species (H. cinerea and H. chrysoscelis/versicolor). Additionally, we found Myxococcota among skin microbiome samples from four of our host species (excluding H. cinerea). Members of this phylum colonize soil and rotting plant material, including decaying wood and tree bark, but also freshwater environments (Reichenbach, 1999). Hyla cinerea is an arboreal frog which is found most often on green, weedy vegetation on the margin of swamps and ponds (Sievert and Sievert, 2021), possibly limiting their exposure to bacteria that inhabit soil and decaying plant material, such as Myxococcota. In contrast, H. chrysoscelis/versicolor inhabits woody vegetation such as shrubs and trees (Sievert and Sievert, 2021), which could explain why Myxococcota was found on the skin of H. chrysoscelis/versicolor and not H. cinerea.

After a recent splitting of the bacterial class Deltaproteobacteria, Myxococcota and Desulfobacterota emerged as two distinct phyla (Waite et al., 2020; Langwig et al., 2022). Desulfobacterota was detected among a few A. blandchardi and R. catesbeiana samples. This phylum has been found in Oklahoma wetland sediment previously (Murphy et al., 2021), which may be why we detected members of this phylum on the skin of A. blandchardi and R. catesbeiana, two host species that are most often found within or around waterbodies, and emphasizes the importance of collecting environmental samples alongside host microbiome samples. To our knowledge, no studies have reported Crenarchaeota, Desulfobacterota, or Myxococcota as dominant members of the amphibian skin microbiome. However, the recent reclassification of Deltaproteobacteria could explain why Desulfobacterota and Myxococcota have not yet been detected on the skin of amphibians. Deltaproteobacteria was found on the skin of frogs from Madagascar in small relative abundances (1.6%; Bletz et al., 2017), and among the different environmental samples collected by Bletz et al. (2017). Therefore, more studies are needed to determine if Crenarchaeota, Desulfobacterota, and Myxococcota are prominent members of the amphibian skin microbiome or if our findings are an example of host frogs obtaining some of their skin microbiota from their surrounding environment, a concept which has been discussed in numerous other studies (Belden and Harris, 2007; Becker et al., 2014; Loudon et al., 2014b; Bletz et al., 2017; Jiménez and Sommer, 2017; Kueneman et al., 2019; Rebollar and Harris, 2019; Hernández-Gómez et al., 2020; Barnes et al., 2021). For example, soil and water can act as reservoirs for microbiota (Loudon et al., 2014b) which can then colonize the skin of amphibians (Belden and Harris, 2007; Becker et al., 2014; Hernández-Gómez et al., 2020; Barnes et al., 2021). It is thought that these reservoirs support the presence of rare bacterial taxa within the amphibian skin microbiome (Loudon et al., 2014b) and changes to microbial communities in the environment may impact amphibian skin microbiome composition and diversity (Belden and Harris, 2007). Given that the skin microbiome is critical for promoting amphibian immune system function and pathogen defense (Bernardo-Cravo et al., 2020; Rebollar et al., 2020), altering reservoir microbiota with captive breeding or head-start initiatives could impact host health (Belden and Harris, 2007; Becker et al., 2014; Loudon et al., 2014b). As such, future investigations which explore the exchange of rare microbes between the environmental and amphibian skin microbiomes are essential for understanding this complex dynamic and may have critical implications for conservation initiatives. To accomplish this, future studies should consider sampling, and subsequently sequencing, the environmental microbiome present at the exact location where the host organism is collected, a method which has been employed in several other studies to date (Walke et al., 2014; Bletz et al., 2017; Jiménez and Sommer, 2017). Implementing this methodology in future studies may improve our understanding of the mechanisms shaping the composition and diversity of anuran skin microbiomes while also providing a way to differentiate skin microbiome samples that have been potentially contaminated with environmental microbes from those collected from host organisms that have rare microbes within their skin microbiomes. Previous findings indicate that the same dominant microbial phyla comprise both environmental and amphibian skin microbiome samples (e.g., Acidobacteria, Actinobacteria, Bacteroidetes, Firmicutes, Proteobacteria) but in different relative abundances (Walke et al., 2014; Bletz et al., 2017; Hernández-Gómez et al., 2020; Barnes et al., 2021). Therefore, analyzing environmental microbiome samples alongside skin microbiome samples may help identify samples that have been contaminated with environmental microbiota because the compositions of the contaminated samples will be more similar to the environmental samples than to the focal skin microbiome samples. In contrast, we would expect that samples collected from hosts with rare microbial taxa in their skin microbiomes would have compositions that are more similar to the other hosts when compared to the environmental samples but with the addition of the rare microbial taxa.

Another consideration for future work is to expand the taxonomic scope of amphibian skin microbiome studies, particularly in a state like Oklahoma which is home to an impressive diversity of amphibian species distributed across a diverse landscape. For the purposes of our study, we aimed to sample co-distributed anuran species that are found commonly in the wild to allow for sufficiently large sample sizes and more robust comparisons of skin microbiomes among taxa. However, continued field surveys across central and eastern Oklahoma would enable expanded comparisons, both taxonomically (across orders, families, genera, and species) and geographically (across more ecoregions) in Oklahoma ( Figure 1 ). For example, Oklahoma is home to several species of true toads (family Bufonidae), narrow-mouthed toads (family Microhylidae), and spadefoot toads (family Scaphiopodidae), along with other members of the focal families Hylidae and Ranidae from this study that lacked sufficient sample sizes to include in our dataset, such as Pseudacris clarkii, P. crucifer, P. fouquettei, and P. streckeri in the family Hylidae and Rana areolata, R. blairi, R. clamitans, and R. palustris in the family Ranidae (Sievert and Sievert, 2021). There is also an incredible diversity of salamanders (order Caudata) distributed across the eastern portion of the state which presents exciting opportunities for future research to compare and contrast the findings of this study with the skin microbiome diversity of this divergent clade of amphibians. Further, given the variable nature in which we sampled frog communities from the four focal ecoregions—frogs from certain ecoregions were sampled across multiple years, other sampled in only two years—future research should investigate to what degree anuran skin microbiomes vary across seasons or years (Douglas et al., 2021).

Symbiotic microbial communities serve fundamental roles in a variety of processes that promote host health. As such, investigating the environmental and evolutionary mechanisms that drive the assemblage of frog skin microbiomes remains essential as we begin to form a more comprehensive understanding of the factors impacting anuran health, persistence, and immunity. Results from this study suggest that multiple factors influence the structure of anuran skin microbial communities; however, host ecology plays the largest role in the assembly of skin microbiome diversity among our sampled Oklahoma frog communities. Comparative studies of skin microbiomes among co-distributed anuran species spanning complex landscapes allow us to better understand, and predict, impending shifts in anuran skin microbial communities in response to continued climate and habitat changes. The investigation of amphibian skin microbiomes remains critical as the vertebrate group continues to face population declines imposed by disease, habitat loss, pollution, and climatic changes.

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number are: NCBI SRA; PRJNA1024480.

The animal study was approved by University of Oklahoma’s Institutional Animal Care and Use Committee. The study was conducted in accordance with the local legislation and institutional requirements.

SS: Conceptualization, Data curation, Formal Analysis, Investigation, Methodology, Project administration, Validation, Visualization, Writing – original draft. JW: Data curation, Funding acquisition, Resources, Writing – review & editing. CS: Conceptualization, Funding acquisition, Resources, Supervision, Validation, Writing – review & editing.