Dorsal Skin Microbial Assemblages in Free-ranging Hellbenders, Cryptobranchus Alleganiensis, Associated With Subspecies and Chronic Toe Lesions, but Not With Chytrid Fungal Infection.

Skin microbiomes are important components of skin health and have been shown to contribute to immunity in amphibians, especially against the chytrid fungus, Batrachochytrium dendrobatidis (Bd). Hellbenders (Cryptobranchus alleganiensis) are large aquatic amphibians (Order Caudata) native to the eastern United States that have experienced population declines of both the Ozark and eastern subspecies, C. a. bishopi and C. a. alleganiensis, respectively. In addition, ulcerative non-healing toe lesions have become increasingly prevalent in C. a. bishopi, in Arkansas (AR) where populations are now reduced to a single river. To evaluate the potential impacts of both chronic toe lesions and Bd on hellbender health, we compared dorsal skin microbial assemblages based on 16S rRNA gene amplicons between a declining Ozark hellbender population in Arkansas (AR) presenting lesions and a reference, recruiting, lesion-free, population of eastern hellbenders in eastern Tennessee (ETN). We further evaluated effects of mass and life stage across both subspecies, as well as toe lesion severity and Bd infection status within AR, to better understand the associations between microbiomes and disease in a wild salamander. Results demonstrate that eastern hellbenders of ETN have richer and less dispersed dorsal skin bacterial assemblages compared to Ozark hellbenders of AR. Furthermore, we suggest that increased severity of toe lesions may be linked to systemic changes resulting in skin microbial dysbiosis, independent of Bd infection. Although lesions remain to have an unknown etiology, this study is another step towards understanding skin bacterial microbiomes in hellbenders, and their potential associations with chronic disease.


Background
The skin is the rst line of defense against invading pathogens and despite a highly developed immune system, thousands of bacterial species are resident on all vertebrate skin surfaces ( Evaluation of altered skin microbiomes may be especially important in amphibians. They rely on their skin for osmotic regulation, which is subsequently more permeable than other vertebrate hosts (Shoemaker and Nagy 1977). In addition, microbial community shifts may be major factors in amphibian host susceptibility to pathogens such as the amphibian chytrid fungus Batrachochytrium dendrobatidis (Bd), which is a fungal dermal pathogen responsible for the extinction of approximately 90 amphibian species worldwide (Scheele et al. 2019).
Intrinsic factors related to the amphibian host are important for colonization and maintenance of these skin bacterial communities (Jiménez and Sommer 2017). In apparently healthy amphibians, bacterial skin communities associate to species across different isolated ponds (McKenzie et al. 2011). Amphibian skin communities can also resist colonization from repeatedly applied bacterial skin washes (Küng et al. 2014). However, other studies reveal amphibian skin bacteria can in fact be shaped by the surrounding environment. Amphibian skin microbial community metrics have been linked to climate and biosphere (Kueneman et al. 2019; Jani 2019), environmental substrate (Fitzpatrick and Allison 2014), and microhabitat type (Bletz et al. 2017a). This dynamic relationship between amphibian host and environment on skin microbiome composition is an important consideration during an assessment of wild amphibian health (Jani and Briggs 2018;Bernardo-Cravo et al. 2020). Associations between environmental degradation and altered host-associated microbial communities are documented in other wildlife species such as black howler monkeys (Amato et al. 2013) and corals (Krediet et al. 2013). In corals particularly, these host-associated microbial community changes have been further linked to decreased coral immunity (Glasl et al. 2016).  Nickerson et al. 2011). In Ozark hellbenders of Arkansas (AR), over 93% of individuals had progressive distal limb lesions that ranged from toe swelling to ulceration and digital necrosis (Hardman et al. 2020a). Toe lesions are documented to occur and worsen over several years. There is no identi ed etiology but disease is likely manifesting from a multifactorial process (Hardman et al. 2020a). These lesions may be linked to environmental factors, although, host and infectious factors are hypothesized to contribute there is a paucity of research on how these surface communities may be related to wound healing. Only one study has evaluated foot lesions in a small sample of Ozark hellbenders and found bacterial assemblages on toe lesions were richer and had differential abundance of core taxa when compared to the dorsal skin (Hernández-Gómez et al. 2017b). It will be important to also understand how lesion presence and severity may be associated with overall changes in the skin microbiome to better mitigate this disease. Furthermore, hellbenders are long-lived animals (> 20 years) and present an opportunity for us to understand the long-term effect of chronic disease in a wild amphibian species, particularly related to the association between chronic non-healing wounds and skin microbiomes. Therefore, we evaluated dorsal skin bacterial assemblages of hellbender populations from two rivers of different conservation status: a stable population of eastern hellbenders in eastern Tennessee (ETN) and a declining and diseased population of Ozark hellbenders in Arkansas (AR). We further evaluated whether variation within the unhealthy AR population could be explained by toe lesion severity or Bd infection status. Lastly we evaluated how life stage and mass could help explain any observed variation.

Results
We sampled 63 individuals, 37 from AR and 26 from ETN (Fig. 1). We recovered 2,804,796 sequences assigned to 911 phylotypes after sequence processing in MOTHUR. Sequence counts per sample ranged from 3 to 379,352.
Subsampling and quality control resulted in subsequent removal of 25 samples with less than 4000 sequences, leaving 40 samples remaining for analysis (15 AR and 25 ETN). Within the 15 AR samples, two were from individuals with mild toe lesions, seven with moderate lesions, and six with severe lesions. The 15 AR samples were also screened for the presence of Bd, using quantitative PCR (qPCR). This screening process identi ed six individuals that were Bd positive within the AR samples. We encountered and sampled juveniles and adults in ETN (nine and 16, respectively), but only adults were sampled from AR. Mass and subsequent size class data were available for all individuals sampled for both subspecies with eight small and eight medium adults in ETN and six medium and nine large adults in AR. Mass ranged from 3g -807g with a median of 320g.
Beta dispersion increased with increasing lesion severity for both Morisita-Horn and Sorenson Indices but was not signi cant (p = 0.09 and 0.12, respectively; Fig. 5).
There was no effect of Bd infection status on measures of alpha diversity (  (Table 3). There were three top indicators (IndVal Index stat > 0.9) for ETN which belonged to genus Ferruginibacter, an unclassi ed order of Cyanobacteria (incertae sedis), and an unclassi ed group of Verrucomicrobia (subdivision 3). Only one top phylotype (genus Limnohabitans) was identi ed for AR (Table 3).  Table 3 List of indicator taxa based on indicator species analysis signi cantly associated with dorsal skin of hellbenders within each comparison group based on subspecies (eastern hellbenders (Cryptobranchus alleganiensis alleganiensis) of eastern Tennessee (ETN) and Ozark hellbenders (Cryptobranchus alleganiensis bishopi) of Arkansas (AR)), life stage within ETN (juvenile vs adult), lesion severity within AR ( mild vs moderate vs severe) and Batrachochytrium dendrobatidis (Bd) infection status (Bd + vs Bd -). No taxa were identi ed as signi cantly associated with animals with severe lesions. Taxa with a stat > 0.9 were considered top indicator taxa and are noted in bold. A total of 136 phylotypes were shared between lesion severity groups within AR, with 17 unique to animals with mild lesions, 50 to animals with moderate lesions, and 48 to animals with severe lesions. Increased abundance in Firmicutes and Synergistetes was correlated with moderate and severe lesion groups for AR hellbenders (Fig. 7, Table 2). Indicator species analysis identi ed 27 phylotypes representative of individuals with mild toe lesions, only one for those with moderate lesions, and none for those with severe lesions (Table 3). Top indicator phylotypes (IndVal Index stat > 0.9) for mild lesions were species in the genera Opitutus, Emticicia, Ferruginibacter, Zavarzinella, and an unclassi ed of order Myxococcales. Five of the 27 indicator phylotypes for mild lesions within AR were also indicators for ETN, including three out of the ve top mild lesion indicator taxa Opitutus, Ferruginibacter, and an unclassi ed order Myxococcales. The only indicator species for animals with moderate lesions was an unclassi ed phylotype in the order Burkholderiales.
Although bacterial alpha diversity and assemblage structure was not signi cant between Bd positive and negative individuals, indicator species analysis identi ed one phylotype (Aquabacterium) indicative of Bd negative skin assemblages and ve for Bd positive skin (Table 3). For life stage analysis within ETN, we identi ed a total of 10 phylotypes as indicative of adult skin, and 10 of juvenile skin, none with a IndVal Index stat > 0.9 (Table 3)

Discussion
Our results showed differences in dorsal skin bacterial assemblage structure and diversity between two subspecies of hellbenders from populations not previously sampled. Furthermore, we showed that variation in these communities could be explained by subspecies, severity of chronic lesions in the feet, and partially by life stage. Bd infection status and adult size class (mass) were not correlated with dorsal skin bacterial assemblages. Several factors including both environment and host health may be contributing to these differences. Interestingly, we found a trend of decreased alpha diversity with increased lesion severity with the apparently unhealthy Ozark hellbender population of AR. This, along with subspecies differences, may suggest that skin bacterial diversity is negatively associated with hellbender skin health. This is in contrast to Hernández-Gómez et al.
(2017b) who found toe lesions to have richer assemblages than apparently healthy dorsal skin within a given individual. However, our study differed in that we compared dorsum between individuals of varying toe lesion severity, as opposed to direct comparison between lesioned and lesion-free skin within an individual. Furthermore, along with decreased diversity, we showed skin microbiomes of AR hellbenders, particularly those with severe lesions, had increased beta dispersion. This suggests severe lesions may be associated with uneven, random, and potentially dysbiotic, skin bacterial communities. In contrast to subspecies comparisons, dorsal skin bacterial assemblage differences between AR individuals of varying toe lesion severity were driven by relative abundance of more common phylotypes. A similar relationship was observed between dorsum and toe lesions in Ozark hellbenders of MO in that differences were driven by changes in relative abundance of shared OTUs (Hernández-Gómez et al. 2017b).
These distal limb lesions are progressive and can take several years to increase in severity (Hardman et al. 2020a), and the same factors affecting toe lesions may also be affecting dorsal skin communities. Changes we observed between AR hellbenders of varying toe lesion severity could represent a decreased ability to maintain an optimal microbial community as disease progresses in an individual. Another explanation could be that local microbial changes due to chronic lesions can eventually generate changes over the entire skin surface as has been seen around chronic ulceration due to leishmaniasis in humans (Gimblet et al. 2017). Alternatively, shifts in microbial communities could be the main mechanism by which lesions manifest, and lesion severity represent a sequela of this dysbiosis that manifests speci cally in an area of increased use (e.g. toes). The mechanism by which these patterns may arise is still speculative and it is possible these trends are not causal or are potentially spurious from limited sample size within AR (n = 15). More intensive within-river surveys alongside more detailed measurements of health (i.e., blood and immunological parameters, body condition score, and detailed lesion scoring) will be needed to better understand within watershed variation and drivers from host health and environment on changes in skin bacterial community parameters.
Unexpectedly, infection status of Bd had no effect on beta or alpha diversity. This was an unexpected result based on other studies of Bd infection being associated with skin microbial community shifts in amphibians ( Life stage also had minimal effect on the bacterial skin community, at least within the ETN population for which we were able to sample juveniles. Similar to subspecies differences, community membership but not relative abundance, was predictive of hellbender life stage. Our study was limited in the fact that only two juveniles were still gilled larvae (pre-metamorphosis). Evaluation of larval and post-metamorphic juvenile substages within a single site may reveal signi cant microbiome shifts after metamorphosis we were unable to evaluate.
We found commonalities in indicator taxa that overlapped with those found in previous hellbender studies. The top indicator phylotypes for AR (genera Limnohabitans and Acinetobacter) belonged to families Comamonadaceae and Moralleaceae, respectively. These taxa were previously identi ed with increased abundance in toe lesions compared to dorsal skin in Some phylotypes may associate with healthy hellbender skin. Several of the indicator taxa that distinguished dorsal skin assemblages of AR animals with mild versus moderate or severe toe lesions were the same that distinguished ETN from AR. One particular (genus Ferruginibacter) has also been identi ed as a core taxon on alpine newt (Ichthyosaura alpestris) skin (Bletz et al. 2017b).
No indicator phylotype for ETN eastern hellbenders from this study overlapped with taxa previously identi ed as signi cantly associated with eastern hellbender skin. As mentioned above, our study evaluated eastern hellbenders from ETN whereas the previous studies comparing hellbender subspecies evaluated those from MO. Populations from ETN and MO vary greatly in phylogeny, ecoregion, and population health. Conversely, AR Ozark hellbenders are phylogenetically very similar to MO Ozark hellbenders and individuals sampled from our study were from one of the same watersheds sampled in these previous studies. It is no surprise then that our indicator species analysis revealed similar results for Ozark hellbenders yet disparate results for eastern hellbenders. The importance of speci c taxa of bacteria and their function within a community context remains an understudied area in amphibian microbiome research. This study contributes to our understanding of the microbial assemblages of hellbenders and their correlation with skin health.

Conclusions
Our study demonstrated that bacterial assemblages on hellbender dorsal skin are not only associated with subspecies but could be shaped by skin health. Regardless of the immediate mechanism driving hellbender skin microbial communities, continued evaluation of their structure and diversity, when added to traditional amphibian disease surveillance, may provide a more comprehensive assessment of wild population health. We found a correlation between decreased dorsal skin bacterial assemblage similarity and richness within Ozark hellbenders compared to eastern hellbenders. This trend was further explained by severity of toe lesions in Ozark hellbenders. This pattern was not affected by presence of Bd on the skin and was likely related to other factors linked with chronic changes in the host and/or aquatic environment. These salamanders are long lived, river-dependent organisms with documented population declines associated with habitat changes, especially sedimentation. The skin environment is constantly exposed overtime to both this altered lotic environment as well as its secondary effects on host skin maintenance. Intensive sampling will be required in both healthy and unhealthy populations of both subspecies to disentangle factors affecting these important skin microbial communities, and to ultimately inform management to maintain and promote future healthy wild hellbenders.

Methods
We sampled hellbenders from 2014-2015 from mid-July to mid-August in two rivers, one from the Blue Ridge ecoregion of eastern Tennessee (ETN), the other in the Ozark Highlands of Arkansas (AR) (Fig. 1). We sampled in ETN for the eastern hellbender (C. a. alleganiensis) and in AR for the Ozark hellbender (C. a. bishopi). All sampling in Tennessee with approval from Tennessee Wildlife Resources Agency (TWRA) (Permit # 1877 (RHH)) and in AR under observation of the USFWS permit TE66039A-0 issued to KJI and Arkansas Game and Fish Commission. In TN, and shallower water in AR, we used standard snorkeling techniques to locate individuals. In AR, we also sampled hellbenders via arti cial nest boxes (Briggler and Ackerson 2012) and sampled deep water habitats (up to 4 meters) using a hookah dive system (i.e. gasoline powered air compressor with tethered air supply lines to dive regulators).
We captured any animals encountered under cover objects and nest boxes and placed them in a clean, soft cotton or mesh bag. We kept bags submerged in the river before and after animals were processed. (University of Tennessee IACUC protocol # 2481 − 0916). We changed dive gloves between animal captures to reduce pathogen transfer and contamination of samples. Skin microbiome samples were collected from a subset of hellbenders surveyed for disease and amphibian pathogen prevalence (see Hardman et al. 2020 a,b). We placed each individual in a new, cleaned plastic tub and rinsed the dorsum with approximately 200 mL sterile distilled water. We immediately rubbed a sterile cotton swab (Fisherbrand product# 23400111) over the dorsal skin surface for 30 seconds and placed it in a sterile 1.5 mL tube on dry ice for transport to a -80°C freezer. For individuals in AR only, we also obtained a swab for Bd load analysis for a concurrent pathogen prevalence study (see Hardman et al. 2020b), as well as, lesion severity score (scored from 0-7) for toe lesions of unknown etiology (Hardman et al. 2020a). All individuals sampled from AR for the microbiome analysis had toe lesions, but because of low sample size within AR, and subsequent inability to compare across all 8 lesion scores, we binned scores of 0-7 into three biologically relevant severity levels of mild (score 0-2; toe swelling only), moderate (3-4; toe swelling and some shortened) and severe (5-7; most toes shortened and/or missing).

Statistical Analysis
We performed all subsequent cleanup and analyses using RStudio v1.2.5019 (R Studio Team 2016). We removed rare phylotypes with less than ve total sequence reads or occurring in less than three samples. We also removed samples with less than 4000 sequence reads based on rarefaction curves performed in package Vegan (Oksanen et al. 2015).
To compare alpha diversity among samples we used packages Vegetarian (Charney and Record 2012) and Vegan to calculate an averaged 1000 iteration bootstrapped value generated from 1000 sequence subsamples for the following Hill numbers: relativized richness (q = 0), exponential of Shannon's index (q = 1) and the inverse Simpson's index (q = 2). We used Hill numbers as they provide better evaluation of changes in both richness and evenness (Chao et al. 2014). We performed the same bootstrap technique to calculate beta diversity distance matrices of the Sorenson Index and Morisita-Horn to evaluate structural changes based on phylotype presence/absence and relative abundance, respectively. We performed all statistical comparisons of both beta and alpha diversity in package Vegan.
For alpha diversity analysis we normalized each response variable (either natural log or cubed transformation) and performed an analysis of variance (ANOVA) across all samples to test for variation in alpha diversity by subspecies and mass. We performed a second ANOVA from AR samples only to test for variation in alpha diversity by Bd infection status and lesion severity within this population. A nal ANOVA was performed to evaluate variation in alpha diversity that correlated with life stage effects (juvenile vs adult) within ETN only where both adults and juveniles were encountered. We used PERMAVOVA (adonis − 999 permutations) to test for differences in average assemblage structure using Sorenson Index and Morisita-Horn distances and previously mentioned categorical predictors mentioned for the alpha diversity analysis. The only exception was that mass was binned into natural breaks of adult size class of small (150-299g), medium (300-449), and large (450+) adults. We also assessed beta diversity as a multivariate dispersion (community consistency among samples) by calculating the multidimensional area of minimum convex polygons t to clusters using the 'betadisp' function within vegan based on distances created from principal coordinates analysis (PCoA) from previously calculated Sorenson and Morisita-Horn indices.
For visualization of alpha diversity differences, we created boxplots for each Hill number for each parameter using package ggplot2 (Wickham 2016   Alpha diversity comparison of bacterial 16S assemblages from dorsal skin swabs of hellbenders (Cryptobranchus alleganiensis) using Hill numbers from order q=0 (relativized richness), q=1 (exponential of Shannon Index) and q=2 (inverse Simpson's Index). Comparison between eastern hellbenders (C. a. alleganiensis) of eastern Tennessee (ETN) and Ozark hellbenders (C. a. bishopi) of Arkansas (AR). ETN hellbenders are further grouped into juveniles and adults.

Figure 3
Alpha diversity comparison of bacterial 16S assemblages from dorsal skin swabs of Ozark hellbenders (Cryptobranchus alleganiensis bishopi) using Hill numbers from order q=0 (relativized richness), q=1 (exponential of Shannon Index) and q=2 (inverse Simpson's Index). Comparisons between individuals with varying toe lesion severity (mild, moderate, severe) in top row, and individuals with differing Batrachochytrium dendrobatidis (Bd) infection status (Bd + vs Bd -) in bottom row.  NMDS ordination plots of hellbender skin bacterial assemblages based on Sorenson and Morisita-Horn Indices. Each point represents a single skin bacterial assemblage from dorsal skin swabs of Ozark hellbenders (Cryptobranchus alleganiensis bishopi) of Arkansas (AR). Top plots grouped by toe lesion severity of mild (yellow), moderate (orange), or severe (red). Bottom plots groups by Batrachochytrium dendrobatidis (Bd) infection status of Bd + (green) or Bd -(gray).

Figure 6
Phylum abundance of the 28 bacterial phlya identi ed from 1000 sequence subsamples of 16S sequences obtained from dorsal skin swabs from eastern hellbenders (Cryptobranchus alleganiensis alleganiensis) captured in eastern Tennessee (ETN) and Ozark hellbenders (Cryptobranchus alleganiensis bishopi) in Arkansas (AR). ETN hellbenders are further grouped by life stage of juvenile (ETN-Juv) and adult (ETN-AR). Only adults were encountered for AR (AR-Ad). Figure 7