Beyond Composition: Network Topology Reveals Batrachochytrium Dendrobatidis -Linked Changes in Salamander Skin Microbiomes

1. Introduction

The microbiome is increasingly recognized as a key determinant of host health, contributing to immune regulation, pathogen resistance, and metabolic homeostasis in vertebrates [1,2,3,4,5,6,7]. However, translating microbiome knowledge into effective disease mitigation remains challenging because microbial communities are dynamic systems whose composition and function depend on both host traits and environmental context [8,9,10,11]. In amphibians, the skin microbiome is particularly relevant as it constitutes a primary interface with the environment. It can directly modulate susceptibility to skin-associated pathogens through competitive exclusion, niche occupancy, and production of antifungal metabolites [12,13].

Microbiome variation is shaped by multiple factors, including habitat, climate, developmental stage, host species, and sex [8,9,10,11]. Consequently, characterizing “who is there” (community composition) is necessary but often insufficient for understanding microbial community behavior under infection or perturbation. An additional layer of insight is obtained by examining how taxa are organized into association structures, which may reflect shared niches, co-dependence, or antagonistic relationships. Network-based approaches have therefore become helpful in summarizing community organization, identifying modules and highly connected taxa, and evaluating whether perturbations lead to fragmentation, re-wiring, or shifts in putative community hubs [14,15,16]. Under infection or pathogen exposure, microbial networks may lose connections, change their modular structure, or display altered association strengths and topologies; conversely, new associations may emerge as communities reorganize [17,18]. Importantly, such structural changes can occur even when the overall community composition appears to be stable.

A primary amphibian disease motivating microbiome-informed research is chytridiomycosis, caused by the fungal pathogen Batrachochytrium dendrobatidis (Bd), which has contributed to severe amphibian population decline globally [19]. Conservation concerns are heightened because disease impacts frequently interact with stressors, such as climate variability, pollution, and habitat fragmentation, potentially increasing host vulnerability and limiting recovery [20,21,22,23,24,25]. Moreover, the emergence and spread of additional fungal pathogens, including Batrachochytrium salamandrivorans (Bsal), further underscore the need for predictive and mechanistic frameworks to understand host–microbiome–pathogen dynamics [26,27,28]. Despite extensive efforts, there is still no universally effective and safe field-ready intervention for chytridiomycosis; tested strategies include thermal treatments and probiotic inoculations, often with variable outcomes across host taxa and environments [29,30,31,32].

One central challenge is that treatment attempts to target infections embedded in a complex and interacting microbial ecosystem. Cutaneous bacteria can facilitate the colonization of beneficial community members, compete for resources, or inhibit other microorganisms via secondary metabolites, producing context-dependent and potentially non-linear outcomes [33,34,35,36,37,38]. Therefore, understanding community organization and how it changes with Bd detection may help identify candidate taxa or community features associated with resistance, tolerance, or dysbiosis and guide the rational design of microbiome-informed mitigation strategies [39,40]. However, the extent to which Bd is related to changes in the structure and composition of salamander skin bacterial communities remains unclear. Prior work has reported taxa that differ in abundance across lesion types or body regions and highlighted particular antifungal bacteria (e.g., Janthinobacterium lividum) as potential protective symbionts [41,42,43,44,45,46,47,48]. In contrast, evidence for consistent opportunistic bacterial signatures that may facilitate Bd infection or reproducible network-level rearrangements linked to Bd status is still limited.

Here, we investigated whether Bd detection is associated with compositional and structural differences in the skin bacteriome of salamanders. We focus on three species that differ in reported susceptibility or outcomes to Bd and Bsal: Eurycea bislineata, Notophthalmus viridescens, and Desmognathus monticola [49,50]. Using amplicon sequencing of skin bacterial communities, we compared Bd-detectable and Bd-non-detectable individuals within and across species and inferred microbial association networks to evaluate changes in community organization. Specifically, we addressed two core questions: (i) to what extent do community composition and network structure depend on host species, and (ii) are Bd-detectable individuals characterized by consistent genus-level shifts and/or network-level reorganization relative to Bd-non-detectable individuals?

2. Materials and Methods

2.1. Data Source and Study Selection

A literature search was performed in PubMed to identify publicly available amplicon datasets suitable for the comparative analysis of salamander skin bacteriomes under contrasting Bd detection statuses. The eligibility criteria were as follows: (i) 16S rRNA gene amplicon sequencing, (ii) consistent targeted the hypervariable region (V3-V5) of the 16S rRNA gene and sequencing platform MiSeq across samples, and (iii) availability of Bd detection data allowing classification of individuals into Bd-detectable and Bd-non-detectable groups. Raw sequence data were obtained from Muletz-Wolz et al. (2022), Inhibitory Bacterial Diversity and Mucosome Function Differentiate Susceptibility of Appalachian Salamanders to Chytrid Fungal Infection, available under NCBI BioProject PRJNA659464. We analyzed samples from three salamander species (Eurycea bislineata, Notophthalmus viridescens, and Desmognathus monticola). Samples were stratified into two conditions based on Bd detection status (Bd-detectable vs. Bd-non-detectable). The final dataset included 15 Bd-detectable and 17 Bd-non-detectable samples.

2.2. Amplicon Processing and Taxonomic Assignment

Raw reads were processed in R using DADA2 for quality filtering, denoising, chimera removal, and amplicon sequence variant (ASV) inference. Taxonomic assignment was conducted within the phyloseq package using the SILVA v138 reference database and DECIPHER classifier. Downstream analyses were performed on genus-level agglomerated tables when required (e.g., network inference), and diversity analyses were conducted at the ASV level unless otherwise specified.

2.3. Community Composition and Diversity Analyses

Relative abundance profiles were computed to characterize both dominant and low-abundance taxa across the samples. Alpha diversity was quantified using the Shannon, Simpson, and Chao1 indices. Differences in alpha diversity between the Bd-detectable and Bd-non-detectable groups were evaluated using the Wilcoxon rank-sum test (Mann–Whitney U test).

Beta diversity was assessed using Bray–Curtis dissimilarities and visualized using principal coordinate analysis (PCoA). Differences in community composition between groups were tested using permutational multivariate analysis of variance (PERMANOVA) implemented with adonis2 (vegan package; 9,999 permutations). To assess within-group dispersion (i.e., community heterogeneity within conditions), we evaluated variance/dispersion in the multivariate space using appropriate dispersion analyses (e.g., beta dispersion) before interpreting the PERMANOVA results. All statistical analyses were conducted in R.

2.4. Differential Abundance and Biomarker Discovery

Differentially abundant taxa between Bd-detectable and Bd-non-detectable samples were identified using two complementary approaches. First, Linear Discriminant Analysis Effect Size (LEfSe) was applied to detect taxa with consistent differences between conditions, reporting discriminant scores at the family and genus levels. Significance was assessed using the Kruskal–Wallis test as implemented in the MicrobiomeMarker package. Second, ANCOM-BC was used at the family and genus levels to account for compositionality and correct for bias in the differential abundance estimation. For ANCOM-BC, taxa were ranked using the W statistic (or the method’s corresponding test statistic), providing an additional and more conservative validation of candidate biomarkers.

2.5. Microbial Association Network Inference

Microbial association networks were inferred at the genus level to compare community organization across the conditions. Networks were inferred separately for Bd-detectable and Bd-non-detectable samples, and within each salamander species to assess species-specific patterns. Network inference combined mutual information-based inference (ARACNE-like approaches) and sparse graphical model inference, as implemented in R packages such as minet and SPIEC-EASI. Associations were estimated using rank-based correlations (Spearman) and/or model-based sparse inverse covariance estimation, and edges were retained according to stringent significance/selection criteria to reduce the number of spurious links. Networks were visualized using Gephi and Cytoscape software. For each inferred network, we computed standard node centrality measures, including degree, betweenness, and closeness centralities, to identify highly connected and potentially influential taxa (putative hubs).

2.6. Network Topology and Robustness Analyses

To quantify the structural differences between the networks, we evaluated their robustness using the NetSwan package. The program calculates the robustness of each graph after a node is removed by determining three important metrics. (i) SwanEfficiency: measures the change in the sum of distances between all node pairs when a node is excluded. (ii) SwanCloseness: Measures the change in the sum of inverse distances (closeness) when a node is removed. (iii) SwanConnectivity: calculates the loss of connectivity (i.e., the number of connected node pairs) after a node is removed. Robustness was assessed by progressively removing nodes under multiple removal strategies: (i) random removal to estimate baseline resilience and (ii) targeted removal based on centrality rankings (degree and betweenness) to model vulnerability to the loss of highly connected taxa and quantify cascading effects. Robustness metrics were compared between Bd-detectable and Bd-non-detectable networks, and where applicable, across species-specific networks). For robustness parameters, group comparisons were performed using one-way ANOVA, followed by Tukey’s post hoc tests when normality assumptions were met (Shapiro–Wilk test); otherwise, non-parametric alternatives were used.

3. Results

Across all samples, the skin bacteriome was dominated by Proteobacteria, with Actinobacteria consistently representing the second most abundant phylum. This pattern was observed in both Bd-detectable and Bd-non-detectable salamanders (Figure 1, top panel). Visual inspection of relative abundance profiles did not reveal marked shifts in the dominant phyla between conditions, suggesting a broadly comparable high-level community composition across Bd statuses.

At a finer taxonomic resolution, the relative abundances at the family level showed highly similar profiles between the Bd-detectable and Bd-non-detectable groups (Figure 1, middle panel). In both conditions, Vibrionaceae was the most abundant family, followed by Halomonadaceae, Pseudomonadaceae, and Alteromonadaceae. No conspicuous condition-specific enrichment or depletion was observed among these dominant families based on relative abundance patterns.

Consistent with the family level results, the genus-level profiles were largely shared across Bd statuses (Figure 1, bottom panel). Vibrio exhibited the highest relative abundance in both groups, followed by Halomonas, Pseudomonas, Alishewanella, and Stenotrophomonas. The dominant taxa remained stable across conditions at the phylum, family, and genus levels, and any differences, if present, were restricted to lower-abundance taxa rather than changes in the most prevalent lineages (Figure 1).

Diversity analyses indicated modest variation between Bd-detectable and Bd-non-detectable salamanders; however, these differences were not statistically significant in the pairwise comparisons (Figure 2). Ordination based on community dissimilarities (PCoA) revealed a partial overlap between conditions, with no consistent separation of Bd-detectable and Bd-non-detectable samples in the multivariate space (Figure 2). Collectively, these results provide no evidence that Bd detection is associated with broad shifts in the overall bacterial diversity or community composition in this dataset.

Community-level tests (PERMANOVA via adonis2) were used to evaluate whether Bd detection status was associated with global differences in community composition, complementing the taxon-level biomarker analyses. PERMANOVA revealed no significant differences in community composition among the groups (R² = 0.01, F = 0.32, p = 1.00, based on 9,999 permutations). These values suggest that there is no significant change in the overall bacterial composition between groups. However, specific differential taxa may be detectable by other methods, such as LEfSe or ANCOM-BC.

To identify taxa associated with Bd detection status, we performed biomarker discovery using LEfSe and validated candidate signals with ANCOM-BC. LEfSe detected a set of discriminative taxonomic features for each condition (Figure 3). In Bd-non-detectable (“healthy”) salamanders, LEfSe highlighted enrichment of Proteobacteria (phylum) and several lower-rank features including Solimonadaceae (family), Salinisphaerales (order), Enterobacteraceae (family), and Cryomorphaceae (family), as well as the genera Lelliottia and Elizabethkingia. In Bd-detectable salamanders, LEfSe identified a larger panel of discriminative features (n = 16), including higher-rank clades/taxonomic labels (e.g., CCD24, GOUTA6) and multiple bacterial families and genera. Among these, Burkholderiaceae (family) and the genera Paucibacter, Lautropia, Altererythrobacter, Methylorubrum, Aureimonas, Shinella, Chloronema, and Amnibacterium were detected as condition-associated features (Figure 3).

To corroborate these targeted compositional differences using a method that explicitly accounts for compositionality and biases in sequencing data, we applied the ANCOM-BC at multiple taxonomic ranks (Figure 4). Consistent with the LEfSe results, ANCOM-BC supported differential abundance of several taxa between Bd-non-detectable and Bd-detectable groups. All reported taxa showed p < 0.05 in LEfSe and exceeded the LDA threshold. This indicates both statistical significance and biological relevance.

In Bd-non-detectable salamanders, ANCOM-BC identified Enterobacteraceae, Solimonadaceae, and Cryomorphaceae (families), Salinisphaerales (order), and the genera Lelliottia and Elizabethkingia as differentially abundant features contributing to group separation. In Bd-detectable salamanders, ANCOM-BC highlighted Chloroflexales (order) and families including Microbacteraceae, Roseiflexaceae, CCD24, and Burkholderiaceae, together with multiple genera distinguishing this condition, including Methylobacterium, Shinella, Amnibacterium, Xylophilus, Aureimonas, Paucibacter, Altererythrobacter, Lautropia, and Chloronema (Figure 4).

ANCOM identified the genera Phreatobacter, Melittangium, Deinococcus, Candidatus omnitrophilus, and Enhydrobacter as having high W = 4 statistics. This indicates consistent differences in log-ratios across multiple taxon comparisons between groups. The bacteria enriched in Bd- are Alterococcus, Angustibacter, Bilophila, Candidatus rhabdochlamydia, Cardiobacterium, Fusobacterium, Verrucomicrobium y Gemella. The bacteria enriched in Bd+ are Candidatus anammoximicrobium, Hypnocyclicus, Lachnoclostridium y Sphingopyxis. In addition, bacteria depleted in Bd- are Stella, Pseudorhizobium, Geobacillus y Kouleothrix. In addition, bacteria depleted in Bd+ are Acidobacter, Aridibacter, Flavimaricola, Methylocella and Nitrincola.

ANCOM-BC identified a set of taxa showing statistically supported differential abundance between Bd-detectable (Bd+) and Bd-non-detectable (Bd−) salamanders (Figure 4 and Figures 4 and S3–S4). Across taxonomic ranks, significant features included multiple genera and higher-level clades such as Phreatobacter, Melittangium, Deinococcus, Candidatus omnitrophus, Enhydrobacter, Myroides, Sumerlaea, Solirubrobacter, Coxiella, Dyadobacter, Aquicella, Paucibacter, Lautropia, Roseomonas, Chloronema, Thermomonas, and additional taxa reported by the classifier (e.g., Candidatus ovatusbacter), as well as higher-rank groups including Rhodothermaceae, Morganellaceae, Fimbriimonadaceae, and others.

When stratifying the results by condition, ANCOM-BC indicated that Bd-detectable samples were characterized by the differential abundance of several higher-rank groups, including the order Azospirillales and the families Bacillaceae, Enterococcaceae, and Leptotrichiaceae. At the genus level, taxa differentiating the Bd+ condition included Asticcacaulis, Candidatus anammoximicrobium, Geobacillus, Lachnoclostridium, Pseudorhizobium, Sphingopyxis, and Stella (Figure 4 and Figures 4 and S3–S4).

In contrast, Bd-non-detectable salamanders exhibited differential abundance of several families, including Cardiobacteriaceae, Chromobacteriaceae, Coleofasciculaceae, Fusobacteriaceae, Gemellaceae, Labraceae, Oscillospiraceae, Parachlamydiaceae, Rhodocyclaceae, Thermomicrobiaceae, and Woeseiaceae. Genera associated with the Bd− condition included Alterococcus, Angustibacter, Arsenicitalea, Bilophila, Candidatus rhabdochlamydia, Cardiobacterium, Crenobacter, Fusobacterium, Gemella, Intestimonas, Rurimicrobium, Spirochaeta 2 (as reported), Suttonella, Tabrizicola, Verrucomicrobium, and Williamsia (Figure 4 and Figures 4 and S3–S4).Overall, ANCOM-BC supports the conclusion that Bd detection is associated with specific taxon-level shifts across multiple ranks, consistent with targeted compositional differences rather than broad restructuring of the dominant community membership.

3.1. Network-Level Differences Between Bd− and Bd+ Salamanders

To characterize community organization beyond taxonomic composition, we inferred genus- and ASV-level association networks separately for Bd-non-detectable (Bd−) and Bd-detectable (Bd+) salamanders (Figures S5 and S6). In the Bd− network, the nodes were predominantly assigned to Proteobacteria (42.25%), followed by Bacteroidota (15.86%), and Actinobacteriota (13.48%). Additional phyla contributed smaller fractions, including Chloroflexi (4.90%), Verrucomicrobiota (4.90%), Myxococcota (4.11%), and Acidobacteriota (3.89%). The Bd− network exhibited relatively strong pairwise associations (most retained edges had absolute association strengths in the ~0.8–1.0 range) and a comparatively even distribution of node connectivity, with most nodes showing approximately 20–25 connections (Figure S5). Visually, the Bd− network was organized into approximately five major modules/clusters, each containing a mixture of phyla, suggesting that dominant lineages contribute across multiple modules rather than being confined to a single compartment.

In contrast, the Bd+ network displayed a more segmented modular organization with more than five clusters/modules (Figure S5). Across modules, Proteobacteria remained prominent and appeared to be overrepresented within most clusters. The phylum-level node composition in the Bd+ samples differed moderately from that in the Bd− samples, with Proteobacteria representing 36.05% of the nodes, followed by Bacteroidota (19.52%) and Actinobacteriota (11.13%). Firmicutes increased to 6.21%, while Verrucomicrobiota (4.52%), Acidobacteriota (4.35%), and Myxococcota(4.35%) contributed smaller fractions. Despite these shifts in modular structure and node composition, the overall node connectivity remained in a range similar to that of Bd−, with typical node degrees of approximately 20–25 (Figure S5). Retained associations in Bd+ samples were also strong (edges largely around ~0.9 in absolute strength), indicating that the Bd detection status was associated more with the reorganization of associations and module structure than with a global increase in network density.

3.2. Positive vs. Negative Associations and Shared ASVs

Next, we compared the balance of positive and negative associations between the Bd− and Bd+ networks (Figure S6; Table S4). The two networks differed in the number of ASVs and taxonomic identities of the condition-specific nodes. Notably, only one ASV was shared across both networks (ASV1), taxonomically assigned to Proteobacteria within the family Vibrionaceae, whereas the remaining ASVs were condition-specific (Figure S6; Table S4).

Bd− salamanders exhibited a higher number of positive associations overall, consistent with a larger set of co-occurrence/co-abundance links within the inferred community structure than Bd + salamanders. In contrast, although Bd+ salamanders showed fewer positive edges, the strength of positive associations tended to be higher, suggesting tighter coupling among subsets of taxa within Bd-detectable communities. Negative associations were observed in comparable numbers across conditions, indicating that putative antagonistic/exclusionary relationships were preserved in terms of overall frequency, even as the positive association structure was reorganized (Figure S6; Table S4).

3.3. Negative Associations and Taxa Linked to Vibrionaceae

Across both inferred networks, taxa assigned to Vibrionaceae exhibited the highest number of negative associations with other network components (Figure 5 and Figure 6). This pattern suggests that Vibrionaceae-associated nodes act as major “antagonistic” or mutually exclusive partners in the inferred association structure, although the ecological interpretation of negative edges should be made cautiously (i.e., they may reflect competitive exclusion, niche partitioning, or compositional effects rather than direct inhibition).

In the Bd-non-detectable (Bd−) network, the strongest/most frequent negative associations involving Vibrionaceae co-occurred with a set of condition-specific ASVs, including ASV73, ASV455, ASV745, ASV809, ASV132, and ASV34, taxonomically assigned to the phyla Myxococcota, Fibrobacterota, Patescibacteria, Armatimonadota, Acidobacteriota, and Bacteroidota, respectively (Figure 5 and Figure 6).

In the Bd-detectable (Bd+) network, negative and positive associations involved a broader set of ASVs (Figure S6; Table S4), including ASV1, ASV9, ASV12, ASV36, ASV72, ASV82, ASV84, ASV87, ASV129, ASV175, ASV179, ASV196, ASV206, ASV266, ASV400, ASV705, ASV859, ASV875, ASV966, and ASV1195. These ASVs mapped to multiple families spanning diverse ecological strategies, including Vibrionaceae, Sanguibacteraceae, Streptococcaceae, Weeksellaceae, Methylacidiphilaceae, Bdellovibrionaceae, Solibacteraceae, Phaselicystidaceae, Spirochaetaceae, Vermiphilaceae, Gemmataceae, Desulfovibrionaceae, Nitrospiraceae, Abditibacteriaceae, and Trueperaceae (Figure S6; Table S4). Together, these results indicate that Bd detection is associated with a more taxonomically diverse set of nodes participating in the signed association structure, even though Vibrionaceae remains a central “negative-edge” contributor in both conditions.

3.4. Species-Level Composition and Shared Taxa Across Bd Status

We compared genus-level community composition across the three salamander species, stratified by Bd detection status, to assess shared (“core”) taxa within each condition. Bd-detectable (Bd+) animals shared 68 genera across species (Tables S1-3), including representative taxa such as Vibrio, Pseudomonas, Halomonas, Stenotrophomonas, Pseudoalteromonas, Aeromonas, Lelliottia, Elizabethkingia, and Mycobacterium (Figure 5). Bd-non-detectable (Bd−) animals shared 70 genera across species (Table S1-3), with recurrent members including Vibrio, Pseudomonas, Halomonas, Stenotrophomonas, Delftia, Microbacterium, Roseomonas, Methylobacterium/Methylorubrum, Elizabethkingia, and Martelella (Figure 5). Overall, both conditions displayed a substantial shared community membership across host species, consistent with a common set of taxa recurrently detected on salamander skin, while network-level analyses (above) suggest that Bd detection is associated with changes in how these taxa are connected within the inferred association structure.

To explore the condition-associated differences in the taxa represented within the modular (community) structure of the inferred genus-level networks, we extracted the list of genera present in each condition-specific network (Bd-detectable vs. Bd-non-detectable) and summarized them by salamander genus (Figure 5; Tables S1–S3). Rather than presenting exhaustive inventories, we report below a set of representative genera for each group; complete taxon lists are provided in Supplementary Material (S1 and S2).

3.5. Eurycea

In Eurycea, the Bd-detectable network contained a subset of genera that were not observed in the Bd-nondetectable network. Representative examples of Bd-associated genera included Treponema, Oceanotoga, Wautersiella, Crossiella, Agrococcus, Iodobacter, and Chitinimonas, along with additional low-frequency taxonomic labels reported by the classifier (Tables S1–S3). In contrast, the Bd-non-detectable network included condition-specific taxa from diverse lineages (Supplementary Tables). Notably, a substantial set of genera was shared across both conditions in Eurycea, including widely detected skin-associated members, such as Vibrio, Pseudomonas, Halomonas, Stenotrophomonas, Pseudoalteromonas, Aeromonas, Lelliottia, Elizabethkingia, Mycobacterium, Delftia, and Microbacterium (Tables S1–S3).

3.6. Desmognathus

For Desmognathus, the Bd-detectable network included several anaerobe-associated and environmentally derived genera as representative members, including Peptococcus, Romboutsia, Anaerotruncus, Terrisporobacter, Lachnoclostridium, Truepera, and Rhodococcus, in addition to additional genera spanning multiple phyla (Tables S1–S3). Bd-non-detectable animals showed a distinct set of taxa (Supplementary Tables), while retaining overlap with common skin-associated genera shared across conditions (e.g., Vibrio, Pseudomonas, Halomonas, and Stenotrophomonas), consistent with a stable core membership with condition-specific additions and losses.

3.7. Notophthalmus

In Notophthalmus, both Bd-detectable and Bd-non-detectable networks shared a broad range of recurrent genera. Representative shared taxa included Vibrio, Pseudomonas, Lelliottia, Elizabethkingia, Martelella, Mycobacterium, Enterococcus, Microbacterium, Halomonas, and Stenotrophomonas (Tables S1–S3). In addition, Bd-non-detectable networks contained taxa such as Roseomonas and methylotroph-associated genera (e.g., Methylobacterium/Methylorubrum) as representative members. In contrast, Bd-detectable networks included additional genera spanning diverse ecological strategies (including taxa related to nitrogen cycling and anaerobic metabolism), suggesting a condition-dependent rearrangement of network membership (Supplementary Tables S1–S3).

Across salamander genera, the inferred network membership displayed a consistent set of widely shared skin-associated taxa (e.g., Vibrio, Pseudomonas, Halomonas, Stenotrophomonas, and Pseudoalteromonas), accompanied by condition-specific genera that differed between Bd-detectable and Bd-non-detectable networks. These patterns support the interpretation that Bd detection is associated with selective gains/losses in network membership and potential reorganization of the network modules. Simultaneously, a substantial subset of taxa remained common across conditions.

Complete lists of each salamander genus and Bd status are provided in Tables S1–S3.

Subsequently, we quantified network modularity for the genus-level association networks to evaluate how strongly each community was partitioned into modules (i.e., groups of genera that were more densely connected than the rest of the network). Bd-non-detectable (Bd−) salamanders showed consistently high modularity across the three host species, indicating a pronounced modular organization. Specifically, modularity values for Eurycea bislineata, Desmognathus (genus), and Notophthalmus viridescens were 0.732, 0.742, and 0.814, respectively (Figure 6). The Bd− networks of Eurycea and Desmognathus were highly similar in modularity, whereas Notophthalmus exhibited a more strongly compartmentalized structure.

In Bd-detectable (Bd+) salamanders, the modularity values were more heterogeneous across host species. Desmognathus and Notophthalmus retained substantial modular structure, with modularity values of 0.695 and 0.851, respectively. In contrast, the Bd+ network for Eurycea bislineata yielded a modularity value of 0, consistent with the absence of detectable modular partitioning under the applied network construction and community-detection settings (Figure 6). This outcome can occur when the inferred network is extremely sparse, lacks sufficient connectivity for stable community detection, or collapses into a trivial structure under the thresholding/selection criteria. Therefore, it should be interpreted cautiously as a structural breakdown rather than a precise estimate of “true” modularity.

Comparing Bd+ vs. Bd− within species, modularity decreased in Bd-detectable animals for Eurycea (0.732 → 0) and Desmognathus (0.742 → 0.695), suggesting reduced modular compartmentalization associated with Bd detection in these hosts. In contrast, Notophthalmus viridescens showed a slight increase in modularity from 0.814 (Bd−) to 0.851 (Bd+), corresponding to a change of +0.037 (Figure 6). Taken together, these results indicate that Bd detection is associated with species-dependent effects on network modular organization, with evidence of modularity loss in some hosts and modest strengthening of the modular structure in Notophthalmus.

3.8. Module Membership Differences Across Conditions

To compare community organization at a finer scale, we analyzed the module (community) structure of each inferred genus-level network. The communities were detected using the Louvain algorithm on the unweighted graph. Modules were identified using a community detection algorithm and labeled by color for visualization. For each network, we extracted the list of genera within each module and compared the module membership between Bd-non-detectable (Bd−) and Bd-detectable (Bd+) conditions.

Across species, the module composition differed markedly between conditions. As an initial standardized comparison, we focused on the first module (colored pink) in each network and summarized the overlap using Venn diagrams across the three host species and two Bd conditions (Figure 6 and Figure 7). We repeated the same comparison, restricting the analysis to Bd− salamanders, to evaluate cross-species similarities in module membership under non-detectable conditions (Figure 6 and Figure 7). These analyses indicated a limited overlap in module composition across conditions, consistent with the condition-dependent reorganization of network communities.

To illustrate this pattern in detail, we examined three representative modules (pink, orange, and green) in Desmognathus monticola networks. Comparing the Bd− and Bd+ networks for this species showed that the genera comprising these modules were largely condition-specific, with minimal to no shared nodes between Bd− and Bd+ within the same module color label (Figures S7–S8; Tables S1–S3). Together, these results suggest that Bd detection is associated with substantial re-assortment of taxa among network communities, even when the overall taxonomic profiles appear broadly similar.

3.9. Core Network Analysis Under Stringent Association Thresholds

Because network topology can depend on the association threshold used to retain edges, we evaluated core networks constructed under more stringent criteria that retained only the strongest associations. The core network was determined by filtering the original network using a fixed correlation threshold greater than 0.9. Additionally, nodes with more than 20 edges were selected. These core networks aim to capture subsets of taxa linked by high-strength associations (i.e., a putative “backbone” of the community).

Core networks differed significantly between Bd− and Bd+ salamanders (Figure 8). In Bd− (non-detectable) salamanders, the core network contained stronger associations (approximately 0.9–1.0), a higher edge count (280 connections), and more nodes (65 genera) than Bd + salamanders. The node degrees were relatively homogeneous, indicating a more evenly connected core structure. In contrast, the Bd+ core network was smaller and sparser, with 96 connections and 35 nodes, and edge strengths were generally lower (approximately 0.7–0.8), as reflected by the reduced edge intensity in the visualization. Additionally, the node degree varied more widely in Bd+, indicating a less uniformly connected core and suggesting a greater concentration of connectivity among a smaller subset of taxa (Figure 8).

Module-level comparisons and core network analyses consistently support condition-dependent restructuring of the inferred association architecture, with Bd detection linked to changes in both community partitioning (modules) and the density/strength of the core association backbone.

Edges are displayed with color gradients reflecting association strength (|ρ| ≈ 0.9–1.0): Bd-detectable (Bd+) networks use a blue gradient, and Bd-non-detectable (Bd−) networks use a pink gradient, with darker/thicker edges indicating stronger associations. For visualization clarity, nodes with fewer than three retained connections (degree < 3) were excluded from the plotted graphs.

To compare the taxonomic membership of the core networks, we constructed a Venn diagram using a list of genera retained in each condition-specific core (Bd− vs. Bd+). The two circles represent the genera present in the Bd− and Bd + cores, and the intersection represents the genera shared by both conditions; the numbers indicate the counts of condition-specific and shared genera (Figure 8).

The core network composition showed minimal overlap between the conditions. Only five taxa were shared between the Bd− and Bd+ cores: SH3-11, Blastococcus, Ruminococcus, Candidatus paracaedibacter, and Deefgea (Figure 8). The Bd− core contained a distinct set of genera, including Oligoflexus, Phreatobacter, Crenobacter, Nocardia, Cardiobacterium, Ruminiclostridium, Fusobacterium, Haemophilus, Syntrophobacter, and Rheinheimera (complete list in Table S2). In contrast, the Bd+ core was characterized by a different suite of taxa, including Candidatus anammoximicrobium, Peptococcus, Gemmata, Kineococcus, Gordonibacter, Robinsoniella, Rubrivivax, Anaerotruncus, Achromobacter, and Macellibacteroides (complete list in Table S1). Together, these results indicate that stringent “core” filtering reveals a condition-specific association backbone, with Bd detection linked to a pronounced turnover in core network membership.

We examined whether the distributions of network centrality and topology metrics differed (i) between Bd detection conditions (Bd+ vs. Bd−) and (ii) among salamander hosts with contrasting reported susceptibility profiles to Bd and Bsal. Specifically, we compared centrality measures across species pairs representing tolerant vs. susceptible outcomes (e.g., Desmognathus vs. Notophthalmus) and tolerant vs. tolerant hosts with different Bsal outcomes (e.g., Eurycea vs. Desmognathus), stratifying all comparisons by Bd status.

For the tolerant–susceptible comparison, Bd-detectable (Bd+) networks of Desmognathus and Notophthalmus exhibited significant differences across multiple metrics, indicating that the Bd+ community organization and node roles differ between these hosts (Figure 9). Notably, several of these differences were reduced or absent when the same species were compared under Bd-non-detectable (Bd−) conditions, suggesting that Bd detection amplifies host-specific differences in the network topology (Figure 9). Among the metrics showing the clearest condition dependence were betweenness centrality and clustering coefficient. High-betweenness nodes often occupy “bridge” positions that connect modules and mediate connectivity across networks. In contrast, the clustering coefficient reflects the tendency of nodes to form tightly connected local neighborhoods (triangles), providing a proxy for community cohesion. Together, changes in these metrics are consistent with Bd-associated reorganization of (i) putative connector taxa and (ii) the local modular structure of the community (Figure 9 and Figure 10). We also analyzed and compared the topological properties and robustness of the networks using the same statistical method with a different pair of susceptible versus tolerant salamanders (see Supplementary Figures S9–S14).

When comparing the two Bd-tolerant hosts (Eurycea and Desmognathus), Bd+ networks also showed significant differences in several topological properties (Figure 10), supporting the idea that Bd detection is associated with species-dependent restructuring of network organization, even among tolerant hosts. Across both tolerant and susceptible comparisons, betweenness and clustering patterns were broadly conserved as informative descriptors of topology; however, their distributions depended strongly on Bd status (Bd+ vs. Bd−). For example, the degree (number of connections per node) was similar between species under Bd− conditions, whereas under Bd+ conditions, degree distributions diverged, with Eurycea showing a significant increase in node connectivity (Figure 11, Figure 12 and Figure 13). This shift suggests a Bd-associated change in the overall network connectedness in Eurycea, potentially reflecting tighter coupling among subsets of taxa or the emergence of more densely connected modules under infection-detectable conditions in the host.These analyses indicate that Bd status and host identity jointly shape the architecture of inferred microbial association networks. Bd detection is linked to changes in the roles of putative connector taxa (betweenness) and local community cohesiveness (clustering), and these effects differ across salamander species with distinct susceptibility profiles (Figure 9, Figure 10, Figure 11, Figure 12 and Figure 13).

Finally, we assessed network robustness to evaluate whether Bd detection was associated with increased vulnerability of the inferred skin microbiome association structure. Robustness was quantified in R using NetSwan by measuring the network connectivity loss as nodes were progressively removed. Networks were analyzed by condition (Bd-detectable, Bd+; Bd-non-detectable, Bd−) and host species using genus-level networks inferred separately for each condition and host species. For each network, we simulated node removal under four strategies: (i) random removal(baseline tolerance to non-targeted loss), (ii) degree-based removal (targeted attack on highly connected nodes), (iii) betweenness-based removal (targeted attack on bridge/connector nodes), and (iv) cascade removal, in which the removal of a focal node triggers additional secondary losses by eliminating nodes/edges that become disconnected or otherwise fail as a consequence of the initial perturbation (i.e., cascading failure). Robustness curves were plotted as the fraction of remaining connectivity (y-axis) against the fraction of removed nodes (x-axis), enabling a direct visual comparison across networks and removal strategies (Figure 11, Figure 12 and Figure 13).

Across conditions, robustness curves showed the expected pattern that targeted removals (degree- and betweenness-based) and cascade perturbations produced a substantially faster decline in connectivity than random removal (Figure 11). In these plots, a steeper early slope indicates more rapid fragmentation and, therefore, lower robustness; conversely, a slower decline under random removal indicates higher tolerance to stochastic losses.

For Bd-detectable networks, degree-based and cascade removal produced broadly similar trajectories, with connectivity decreasing sharply after the removal of a relatively small fraction of nodes. In particular, the removal of fewer than ~25% of nodes resulted in a significant reduction in connectivity (approximately three-quarters of connectivity lost), consistent with a strong dependence on a subset of highly connected or structurally critical taxa. As node removal progressed beyond ~50%, the curves approached a near-zero plateau, reflecting a near-complete collapse of network connectivity(Figure 11).

Bd-non-detectable networks displayed qualitatively similar sensitivity to targeted attacks: degree-based and cascade strategies led to rapid fragmentation beginning with the removal of approximately 25% of nodes. In contrast, under random removal, the Bd− network retained connectivity for a longer period, indicating a comparatively higher robustness to non-targeted node loss. Nevertheless, both Bd+ and Bd− networks were highly vulnerable to targeted removal strategies, implying that their inferred association structures depend disproportionately on a limited set of central/connector nodes (Figure 11).

Robustness analyses indicate that while networks tolerate random loss to some extent, they are structurally fragile under targeted attacks, particularly those aimed at high-degree or high-betweenness nodes, and under cascading failure scenarios. These results are consistent with a topology in which a subset of central taxa disproportionately supports global connectivity and provides a quantitative framework to compare vulnerability across Bd conditions and host species (Figure 11, Figure 12 and Figure 13).

3.10. Susceptible vs Tolerant Hosts (Notophthalmus vs Desmognathus)

In the degree-based removal scenario (targeting highly connected nodes), Bd− networks differed significantly between Notophthalmus viridescens (susceptible) and Desmognathus monticola (tolerant), indicating species-dependent sensitivity to the removal of high-degree taxa even in the absence of detected Bd (Figure 12). This between-species difference was also evident in Bd+ networks, suggesting that vulnerability to high-degree node loss is a consistent discriminator across Bd conditions (Figure 12).

In contrast, random removal did not yield significant differences between groups, consistent with the broadly similar tolerance to stochastic node loss across species and Bd conditions (Figure 12). This pattern supports the interpretation that the most pronounced contrasts emerge under targeted perturbations rather than under nontargeted losses.

For cascade removal, connectivity loss differed by species and Bd condition. Within Desmognathus, Bd+ networks exhibited greater connectivity loss than Bd− networks under cascading failures, whereas Notophthalmus showed the opposite pattern, with Bd+ networks displaying greater vulnerability than Bd− (Figure 12). In absolute terms, cascade removal produced substantial fragmentation in both species, with Desmognathus typically showing ~70–75% connectivity loss and Notophthalmus reaching ~75–80% loss under comparable removal fractions (Figure 12).

The most substantial divergence between tolerant and susceptible hosts emerged under betweenness-based removal(targeting bridge/connector nodes). Notophthalmus Bd+ networks collapsed most rapidly, approaching a complete loss of connectivity (~100%), whereas Notophthalmus Bd− networks retained more structure (~75% loss). Incontrast, Desmognathus networks showed a comparatively stable response across Bd conditions, typically losing ~70–75% connectivity in both Bd+ and Bd− cases (Figure 12). Together, these results indicate that Bd detection is associated with a markedly increased vulnerability of the connector-node architecture in susceptible hosts. In contrast, tolerant hosts maintain more consistent robustness profiles across conditions.

Overall, these patterns provide quantitative evidence that susceptible and tolerant hosts differ in network robustness, particularly under the targeted removal of structurally important taxa (degree and betweenness) and cascade failures (Figure 12).

3.11. Tolerant vs Tolerant Hosts (Eurycea vs Desmognathus)

In the comparison between two Bd-tolerant hosts (Eurycea bislineata and Desmognathus monticola), the robustness profiles were broadly similar across removal strategies and Bd conditions, with typical connectivity losses in the ~70–80% range (Figure 13). Most removal methods did not produce significant between-group differences, suggesting that robust responses are more conserved among tolerant hosts than in tolerant–susceptible comparisons. The primary exception was degree-based removal, in which a significant difference was detected, driven by greater connectivity loss in the Desmognathus Bd− networks (Figure 13).

Taken together, robustness analyses indicate that Bd-associated changes in network vulnerability are host-dependent: the tolerant–susceptible contrast reveals pronounced differences, especially for connector-node (betweenness) removal, whereas tolerant–tolerant comparisons show largely overlapping robustness profiles with only limited differences (Figure 12 and Figure 13).

4. Discussion

Chytridiomycosis caused by Batrachochytrium dendrobatidis (Bd), together with disease caused by B. salamandrivorans (Bsal), remains one of the most critical drivers of amphibian population declines and extinctions worldwide. Increasing evidence indicates that the cutaneous microbiome plays an active role during infection, influencing host resistance, tolerance, and immune function through mechanisms that include competition for space and nutrients, the production of antifungal metabolites, and the modulation of host immune responses [52]. Probiotic strategies aimed at augmenting antifungal bacteria have shown partially favorable outcomes in some systems [45,53]; however, the results are inconsistent across host taxa and environmental contexts. In salamanders, field surveys often report a low prevalence of “classic” anti-Bd taxa and, instead, a substantial representation of commensal and opportunistic microorganisms that may be neutral or even facilitate pathogen establishment under certain conditions [44]. This is consistent with the broader view that amphibian immune responses alone may be insufficient to clear chytrid infection and that durable acquired immunity is limited or context dependent [54].

4.1. Community Composition and Diversity: Stability at High Taxonomic Ranks

In this study, salamanders of the genera Eurycea, Notophthalmus, and Desmognathus exhibited broadly similar high-level community profiles, with Proteobacteria, Bacteroidota/Bacteroidetes, Actinobacteriota/Actinobacteria, and Firmicutes among the predominant phyla across conditions. This pattern resembles that reported in previous studies on salamanders (e.g., Salamandra salamandra). This supports the notion that cutaneous communities frequently converge on a set of dominant bacterial lineages shaped by the shared ecological constraints of amphibian skin (hydration, osmotic conditions, antimicrobial peptides, and environmental exposure) [44]. Importantly, we did not detect significant differences in alpha diversity indices or broad community ordination between the Bd-detectable and Bd-non-detectable groups, consistent with studies reporting that Bd infection in salamanders can occur without major shifts in overall diversity or global composition [44]. Together, these results suggest that Bd detection is not necessarily accompanied by community-wide restructuring at a coarse taxonomic resolution and that disease-relevant changes may instead manifest as targeted shifts in specific taxa and/or in the architecture of microbial associations.

4.2. Differential Taxa and the Ecological Meaning of “Biomarkers”

Although dominant taxa were broadly shared between conditions, both LEfSe and ANCOM-BC identified condition-associated taxa (at the family/genus levels and at some database-specific clades). Such taxa should be interpreted as candidates of ecological relevance rather than direct causal agents. Several differentially abundant genera detected here have been reported in amphibian systems in association with lesions, stress, or opportunistic infections (e.g., Elizabethkingia, Klebsiella, and related Enterobacterales), although species-level resolution and functional verification are required before attributing pathogenicity or protective roles [57]. Conversely, taxa such as Pseudomonas and other Proteobacteria can participate in biofilm formation, metabolite production, and wound-associated processes in other amphibian models, potentially contributing to tissue repair or microbial antagonism, depending on the strain context [58]. These observations reinforce a central challenge in microbiome interpretation: many genera contain both beneficial and opportunistic strains, and functional effects are often strain-specific, context-dependent, and mediated by interactions with the broader community and host physiology.

Environmental exposure is also likely to contribute to the abundance of detected taxa. Amphibian skin communities are strongly influenced by habitat microbial pools (water, soil, and leaf litter), and shifts in temperature, humidity, pollutants, and land use can remodel the available microbial reservoirs and thus the community assembly process [59,60,61]. This environmental dependence may explain why microbiome signatures of Bd detection can vary across host species, locations, and between experimental and natural infections [44,55]. The presence of taxa previously associated with skin damage or ulceration, such as Acinetobacter in some amphibian contexts [63], is compatible with the idea that Bd-induced disruption of the stratum corneum could alter colonization opportunities; however, the present dataset cannot determine whether these taxa are drivers, passengers, or markers of altered skin physiology.

4.3. Network Analyses: Reorganization of Associations Despite Similar Composition

A key contribution of this study is the use of microbial association network analysis to move beyond “who is there” toward “how taxa are organized.” Across conditions, the inferred networks displayed modular organization, and the distributions of several centrality measures were broadly similar at the pooled condition level. This suggests that the community maintains comparable global architectural constraints, even when some taxa shift in relative abundance. Nevertheless, multiple lines of evidence indicate condition-dependent reorganization.

• Module membership turnover: Even when focusing on the same nominal module (e.g., the first/pink module across networks), overlap in node membership across Bd− and Bd+ networks was minimal, implying that Bd detection is associated with a reshuffling of taxa among modules. This is consistent with a community in which different taxa can occupy similar ecological “roles” or niches under varying conditions.
• Core network differences: Under stringent thresholds defining “core” associations, Bd− communities exhibited a larger and denser core than Bd+ communities, with limited overlap in core membership (only a small set of shared taxa). This pattern is consistent with Bd detection being associated with the loss (or reduced detectability) of strongly coupled taxa and/or the weakening of high-strength associations among community members.
• Modularity changes were host-dependent: Bd− networks showed high modularity across host taxa, whereas Bd+ networks showed modularity responses that differed by species. The modularity value of zero observed for Bd+ Eurycea warrants cautious interpretation: it may reflect the collapse of the community structure, but it can also arise from methodological factors (e.g., sparse networks after filtering, small sample size, disconnected components, or algorithmic resolution limits). Nonetheless, considering the broader evidence of module membership turnover and core network shrinkage, our results are consistent with the hypothesis that Bd detection can disrupt modular organization in a host-specific manner.

A critical point is that network edges inferred from compositional sequencing data represent associations rather than direct interactions between nodes. Negative edges may reflect competition, niche differentiation, or compositional constraints, whereas positive edges may reflect co-occurrence driven by a shared environment or metabolic coupling. However, consistent differences in module structure and core density across conditions suggest that Bd detection is linked to changes in community coupling and potential functional redundancy.

4.4. Central Taxa, Clustering, and Implications for Tolerance vs Susceptibility

Comparisons of centrality distributions between tolerant and susceptible hosts suggest that Bd detection may amplify differences in the roles of “connector” taxa (high betweenness) and local cohesiveness (clustering coefficient). Nodes with high betweenness often connect modules and can disproportionately influence network connectivity during perturbations. Simultaneously, clustering captures the local neighborhood structure and is commonly interpreted as a proxy for the capacity to form cohesive sub-communities. The observation that specific centrality differences are stronger in Bd+ than in Bd− networks supports the view that infection status can reshape the taxa that act as bridges and how tightly local neighborhoods are organized. Importantly, these patterns appear more pronounced in tolerant–susceptible contrasts than in tolerant–tolerant comparisons, suggesting that tolerance may be associated with a more stable network organization across conditions. In contrast, susceptibility may coincide with larger infection-associated shifts in the roles of nodes.

4.5. Robustness: Vulnerability Under Targeted Attacks and Ecological Interpretation

Robustness analyses indicated that both Bd− and Bd+ networks were relatively tolerant to random node loss but were vulnerable to targeted removal strategies (degree-, betweenness-based, and cascade failures), consistent with dependence on a subset of highly connected or bridging taxa. In tolerant hosts, robustness profiles were generally more consistent across Bd conditions, whereas in susceptible hosts, they showed stronger Bd-associated changes, particularly under bridge-node removal. These findings align with the hypothesis that host tolerance may be related to maintaining a structurally resilient microbial association architecture, potentially buffering functional processes on the skin during pathogen exposure. Ecologically, multiple perturbations can mimic node-loss scenarios, including host antimicrobial peptides, antibiotics, toxic compounds, or competitive displacement among microbes. However, translating robustness curves into ecological causality requires functional validation and time-resolved sampling.

4.6. Limitations and Future Directions

This study has some limitations. First, this study relied on cross-sectional 16S rRNA gene amplicon data, which limited taxonomic resolution and did not directly capture function or strain-level variation. Second, network inference from compositional data is sensitive to preprocessing choices (filtering, normalization, and thresholding) and should be interpreted as hypothesis-generating, rather than definitive interaction maps. Third, Bd detection status (Bd+ vs. Bd−) does not necessarily capture infection intensity, lesion severity, or temporal stage, all of which could influence microbiome dynamics. Finally, host species differ in their ecology and physiology, and environmental microbial pools may confound the effects of these conditions.

Future work should prioritize (i) incorporating Bd load and disease severity metrics, (ii) longitudinal sampling to distinguish cause from consequence, (iii) shotgun metagenomics or targeted culturing to resolve strain-level functional traits, (iv) explicit testing of antifungal activity and community assembly mechanisms, and (v) integrating environmental microbiome baselines. Together, these approaches will clarify whether the observed condition-dependent network reorganizations reflect adaptive restructuring, pathogen-driven disruption, or environmental filtering, and whether particular taxa or interaction motifs can be leveraged for conservation efforts.

5. Conclusions

This study evaluated whether Bd detection is associated with changes in the salamander skin bacteriome using 16S rRNA amplicon data and a complementary complex network framework. Across Eurycea bislineata, Notophthalmus viridescens, and Desmognathus monticola, Bd-detectable (Bd+) and Bd-non-detectable (Bd−) samples showed no significant differences in broad community composition or standard diversity metrics, indicating that wholesale shifts in dominant taxa or overall richness/evenness do not necessarily accompany Bd detection. Simultaneously, differential abundance analyses identified condition-associated taxa consistent with Bd-linked changes that are specific and taxon-resolved rather than global.

Network-based analyses revealed that Bd detection is associated with the reorganization of the microbial association architecture. Module membership exhibited limited overlap between Bd+ and Bd− networks, and stringent “core” networks differed strongly in density and composition, with only minimal shared memberships. Modularity patterns further suggested host-dependent responses, including a marked loss of detectable modular structure in Bd+ Eurycea, whereas other hosts retained strong compartmentalization. Comparisons of centrality and robustness supported the idea that host susceptibility profiles are reflected in network topology, with tolerant–susceptible contrasts showing larger differences under targeted perturbation scenarios than tolerant–tolerant comparisons.

Bd-related microbiome signatures in salamanders may be subtle at the level of bulk composition but are detectable in the structure of microbial network associations and core connectivity. By integrating differential abundance with network modularity, centrality, and robustness, this study provides a reproducible analytical template for identifying structural biomarkers of dysbiosis and resilience in amphibian skin. This perspective can inform future conservation strategies by prioritizing not only “which taxa are present,” but also how microbial communities are organized, which may ultimately shape stability, tolerance, and responsiveness to microbiome-based interventions.