Experimentally Induced Metamorphosis in Axolotl (Ambystoma mexicanum) Under Constant Diet Restructures Microbiota Accompanied by Reduced Limb Regenerative Capacity

Axolotl (Ambystoma mexicanum) is a critically endangered salamander species and a model organism for regenerative and developmental biology. Despite life-long neoteny in nature and in captive-bred colonies, metamorphosis of these animals can be experimentally induced by administering Thyroid hormones (THs). However, biological consequences of this experimental procedure, such as host microbiota response and implications for regenerative capacity, remain largely unknown. Here, we systematically compared host bacterial microbiota associated with skin, stomach, gut tissues and fecal samples based on 16S rRNA gene sequences, along with limb regenerative capacity, between neotenic and metamorphic Axolotls. Our results show that distinct bacterial communities inhabit individual organs of Axolotl and undergo substantial restructuring through metamorphosis. Drastic restructuring was observed for skin microbiota, highlighted by a major transition from Firmicutes-enriched to Proteobacteria-enriched relative abundance and precipitously decreased diversity. Remarkably, shifts in microbiota was accompanied by a steep reduction in limb regenerative capacity. Fecal microbiota of neotenic and metamorphic Axolotl shared relatively higher similarity, suggesting that diet continues to shape microbiota despite fundamental transformations in the host digestive organs. The results provide novel insights into microbiological and regenerative aspects of Axolotl metamorphosis and will establish a baseline for future in-depth studies.


Introduction
limb regenerative capacity between two developmental stages. The results provide novel 112 insights into relationship between perturbed microbial communities due to host factors 113 and putative implication of microbiota in regeneration capacity. 114 115

116
Details of the experimental design were described in the methods section (Fig. 1). Within 117 2-3 weeks of hormonal treatment of the animals, we observed weight loss, progressive 118 disappearance of the fin and decrease in the gills size; and in approximately two months 119 all animals showed characteristics of accomplished metamorphosis (Fig. 2a). We first 120 performed comparative analysis of microbiota between neotenic and metamorphic 121 Axolotl organs. We then attempted to answer the question whether shifts in microbiota 122 might be accompanied by a change in host regenerative capacity since previous studies 123 provided evidence that limb regenerative capacity is reduced in metamorphic Axolotl 33 . 124 We reproduced these observations in this study that regeneration is indeed impeded in the 125 limbs of metamorphic Axolotl (Fig. 2b). We observed that all neotenic animals (n=15) 126 regenerated (day 64) the limb in a miniaturized form with four digits. In contrast, slower 127 blastema formation and regeneration process were apparent in metamorphic animals. We Sequencing of the V3-V4 region of the 16S rRNA gene produced approximately 3.7M 139 reads generated from 27 samples (24 samples from Axolotls and 3 aquarium water 140 samples, hereinafter "Aqua"). The sequences were clustered into 14451 high quality, 141 singleton, chloroplast, and mitochondria removed and chimera-checked Operational 142 Taxonomic Units (OTUs). Average amplicon sequences per sample was 139059 ± 49159 143 sequences). Our data included 12224 (85% of total) de novo OTUs (OTU IDs that begin 144 with "New.Reference" or "New.Cleanup.Reference"). We then classified a representative 145 sequence of these OTUs using RDP classifier (v. 2.2) at 70% bootstrap cutoff. We 146 identified 621 OTUs that did not find hits in the RDP database even at the phylum level 147 ("Unclassified Bacteria"). We thus used MOLE-BLAST to determine their identity. Species richness and diversity were analyzed using a variety of alpha-diversity metrics 155 across neotenic and metamorphic samples (Fig. 3a- Methylobacterium, Elizabethkingia, Vogesella, Chryseobacterium, and Zoogloea were the 235 top scoring genera in the metamorphic skin while Limnohabitans and several taxa that 236 were classified at higher taxonomic ranking were differentially abundant in the neotenic 237 skin samples (Fig. 6a). These genera were also detected in the indicator species analysis 238 in several OTUs. For example, OTU89, OTU141, OTU1052, OTU1301, OTU1450, 239 OTU2081, were classified as Pseudonocardia and OTU89 had the highest abundance of metamorphic samples, respectively. We also collected water samples ("Aqua") from the 257 Aquarium of Axolotl to identify OTUs present in the water that may have colonized the 258 Axolotl skin (Fig. 7a). We found 89 OTUs were shared by the neotenic and metamorphic 259 skin and Aqua samples; 38 OTUs were present only in the neotenic and Aqua samples 260 while 105 OTUs were shared by metamorphic and Aqua samples. However, relative 261 abundances of all these OTUs were mostly less than 1%. We also compared unique and 262 shared OTUs between gut tissue and fecal samples (Fig. 7b). Exceptionally, the number 263 of OTUs in the metamorphic gut was 4999, representing substantial increase from 2961 264 OTUs found in the neotenic gut. And the two gut tissue samples shared only 227 OTUs. 265 In stark contrast, the number of unique OTUs decreased to 1365 OTUs whereas neotenic 266 fecal sample had 3408 unique OTUs. Notably, the shared number of fecal and gut OTUs, The main purpose of this study was to comparatively characterize bacterial microbiota of 295 Axolotl in neotenic and metamorphic life stages since microbiota of this important 296 biological model has not been reported before. We also asked if differential profiles of 297 Axolotl microbiota in the examined life stages might be accompanied by a change in host 298 regenerative capacity. Our results show that (a)-substantial shifts occurs in the structure 299 and composition of microbiota, particularly in the skin but also in digestive organs; and 300 that (b)-regenerative capacity in the limbs of metamorphic Axolotl steeply diminishes 301 while neotenic Axolotl maintains normal regenerative capacity under identical conditions 302 and that the shifts in the skin microbiota appears to be related to the reduced limb 303 regenerative capacity. Although the scope of this report does not include pinpointing 304 molecular mechanisms accounting for the reduced limb regeneration and its putative 305 association with microbiota, our observations align well with recently published 306 experimental evidence supporting the conclusion that limb regenerative capacity of 307 Axolotl diminishes upon metamorphosis 33,39 . Crucially, our ongoing experiments suggest 308 that reduced regenerative capacity is not observed for some other types of tissues in 309 metamorphic Axolotl (Demircan et al., unpublished data). 310

311
We detected drastic expansion of Proteobacteria in metamorphic Axolotl skin relative to 312 the neotenic skin and noted steep decrease of bacterial diversity in the fecal and skin 313 samples (see Fig. 3). Evidence of skin dysbiosis and the linkage of members of 314 Proteobacteria to reduced regenerative capacity comes from the study by 40 ; in this study 315 pathogenic shifts in microbiota and infections were implicated in reduced regeneration in 316 planaria, a model for tissue regeneration; Although the abundance of Vogesella, 317

Pseudomonas, and
Sphingomonas, all a member of Proteobacteria, and 318 Chryseobacterium within Bacteroidetes phylum were virtually not detected in neotenic 319 planarian skin samples, the abundance of these genera considerably increased in the skin 320 of the metamorphic planarian and these genera demonstrably impeded TAK1/MKK/p38 321 signaling pathway of the innate immunity of the planarian host. Strikingly, our analysis 322 revealed that abundance of these genera and some other members of Proteobacteria 323 known to cause nosocomial infections were, too, highly significantly increased in the 324 metamorphic skin samples from Axolotl (Fig6a and 6b); Vogesella (q=7.4E-11); 325 Chryseobacterium (q=5.4E-12); Undibacterium (q=7.01E-11), and Sphingomonas 326 (q=4.03E-09). The average abundance of Pseudomonas, too, increased (from 0.1% in 327 neotenic skin to 2.2% in metamorphic skin) but the increase was barely significant after 328 FDR correction for multiple testing (q=0.056). Furthermore, infiltration of macrophages, 329 as key players of innate immune system, to the site of amputation and cellular signaling We sequenced pools of water samples from the Axolotl aquarium ("Aqua") to identify 389 water-borne bacterial taxa acquired by Axolotls. Surprisingly, most abundant genera of 390 Aqua samples (e.g. Acidovorax, Armatimonas, Flectobacillus) were not associated with 391 Axolotl organs but only a subset of low-abundance bacteria was detected. For example, 392 Aquabacterium abundance in metamorphic skin and neotenic stomach were 3.9% and 393 7.3%, respectively while its abundance in Aqua sample was 0.1%. Our results are 394 consistent with previous work reporting amphibian skin may select rare taxa from the 395 environment 51-53,57,58 . In this study, we found that 89 out of 509 OTUs in Aqua samples 396 were present in neotenic skin samples while 105 OTUs shared between Aqua and 397 metamorphic skin samples (see Fig. 6), and 150 shared OTUs with Aqua samples, the 398 highest among others, with metamorphic stomach samples. Consequently, our results 399 support the notion that both host and external factors shape the host microbiota but host 400 genetics applies selective filter. 401 402 Conversely, diet is another crucial factor strongly influencing structure of gut microbiota 403 of animal host and even dominate host genotype 59 . Although host genotype and diet are 404 constant in this study, metamorphosis is likely to cause remodeling of epigenetic 405 landscape in the host genome, which in turn is expected to reshape microbiota. Notably, 406 fecal samples from the neotenic and metamorphic Axolotls clustered together, albeit 407 richness in metamorphic fecal samples significantly decreased. We observed that feeding 408 behavior of the metamorphic Axolotl change during metamorphosis, the animals tend to 409 eat less often (low appetite), which might account for the decreased fecal diversity. 410 Although no major restructuring of intestine via metamorphosis was apparent as 411 described before 16 , we observed a higher number of goblet cells and thicker mucus layer 412 in the metamorphic gut tissue compared to neotenic gut ( Supplementary Fig. S8). Taken 413 together, relative influence of diet and host epigenetics seems to be compartmentalized; 414 diet appears to influence strongly the bacterial diversity in the fecal microbiota in the gut 415 lumen whereas the host epigenetics (and the resulting changes in transcriptome due to 416 metamorphosis) seems to play a greater and selective role in gut tissues (crypts). Some 417 genera, often associated with symbiosis such as Alistipes and Elusimicrobium, were 418 differentially enriched in the metamorphic gut tissues whereas neotenic gut tissues were 419 represented by two genera with unclassified Clostridiaceae and unclassified 420 Enterobacteriaceae. Suprisingly, the abundance of Akkermansia considerably increased 421 in the metamorphic gut tissues (16%) relative to neotenic gut (9%). A. muciniphila within 422 this genus is known to be a mucin degrading and symbiotic bacterium and the increase in 423 abundance of this genus is meaningful with increasingly thicker mucus layer we observed 424 in the metamorphic gut tissue staining. 425

426
Our analysis revealed that Axolotl microbiota to certain extent was "humanized" in 427 captivity as manifested from similarity of gut and skin microbiota with the human 428 microbiota, which prompted us to compare predicted functions using PICRUSt 60 , which 429 is optimized for human microbiota. We found that all samples included in microbiota 430 analysis were highly similar based on the abundance of predicted microbial genes (40%). 431 Though further study using shotgun metagenomics technology is warranted, our findings 432 raise the possibility that deficiency in the microbial community function in a given organ regulations. The experimental design followed in this study is depicted in Fig. 1 Nephele pipeline, primers were removed using cutadapt program 63 and the pipeline was 556 stringently configured to perform the following steps; reads below average quality scores 557 (q<30) and read-length >450 bp were eliminated. After joining the pair-end reads, the 558 reads were clustered into operational taxonomic units (OTUs) using open reference OTU-559 picking strategy 64 . The open-reference approach initially runs a closed-reference step. 560 Sequences that fail closed-reference assignment are then clustered as de novo OTU based 561 on pairwise similarity among all sequences in the data set. Open reference clustering was 562 performed based on the 97% clustered SILVA reference (SILVA 123 release; ) database 563 65 and SortMeRNA combined with SUMACLUST algorithms 66 . Non-matching reads to 564 closed reference were subsequently clustered de novo. After obtaining the OTU table 565 from the pipeline, UCHIME (v.4.2) program (http://drive5.com/uchime) integrated into 566 the Mothur (v.1.39.5) tool 67 was separately run to remove chimeric reads. Additionally, 567 The RDP classifier 68 (v. 2.2) was locally run to assign taxonomy for each OTU at a 568 confidence greater than 70% cutoff. Reads that could not be classified at the genus level 569 were sequentially assigned to higher taxonomic hierarchy up to the kingdom level. 570 Unclassified reads at the kingdom level ("Unclassified Bacteria") were extracted from the 571 OTU-representative sequences and searched for nearest neighbor method using MOLE-572 BLAST 69 . This tool computes a multiple sequence alignment (MSA) between the query 573 sequences along with their top BLAST database hits, and generates a phylogenetic tree. 574 Species richness and diversity were estimated by QIIME with the following alpha 575 diversity metrics: OTU richness, Chao1, Shannon, Simpson E, Inverse Simpson, and 576 Faith's phylogenetic diversity (PD). We assessed normality of alpha diversity data using 577 Shapiro-Wilk tests and compared the metrics between neotenic and metamorphic 578 microbiota using unpaired two-tailed t-test. Venn diagrams were constructed using jvenn, 579 a web based tool (http://jvenn.toulouse.inra.fr) 70 . 580 581

Multivariate analysis of community structures and diversity 582 583
Bray-Curtis similarity index 71 and Jaccard index of similarity 72 were used to obtain 584 distance matrix after standardizing by the column sums and transforming (square-root) 585 the read abundance data. Similarities in microbial community structures among samples 586 were first displayed using principal coordinate analysis (PCO) (unconstrained). 587 Differences in community structure related to metamorphosis were displayed using a 588 constrained ordination technique, Canonical Analysis of Principal coordinates (CAP). 589 Tests of the multivariate null hypotheses of no differences among a priori defined groups 590 were examined using PERMANOVA and the CAP classification success rate. CAP uses 591 PCO followed by canonical discriminant analysis to provide a constrained ordination that 592 maximizes the differences among a priori groups and reveals patterns that cannot be 593 unraveled using unconstrained ordinations 73