Changes of intestinal microbiota in the giant salamander (Andrias davidianus) during growth based on high-throughput sequencing

Despite an increasing appreciation of the importance of host–microbe interaction in healthy growth, information on gut microbiota changes of the Chinese giant salamander (Andrias davidianus) during growth is still lacking. Moreover, it is interesting to identify gut microbial structure for further monitoring A. davidianus health. This study explored the composition and functional characteristics of gut bacteria in different growth periods, including tadpole stage (ADT), gills internalization stage (ADG), 1 year age (ADY), 2 year age (ADE), and 3 year age (ADS), using high-throughput sequencing. The results showed that significant differences were observed in microbial community composition and abundance among different growth groups. The diversity and abundance of intestinal flora gradually reduced from larvae to adult stages. Overall, the gut microbial communities were mainly composed of Fusobacteriota, Firmicutes, Bacteroidota, and Proteobacteria. More specifically, the Cetobacterium genus was the most dominant, followed by Lactobacillus and Candidatus Amphibiichlamydia. Interestingly, Candidatus Amphibiichlamydia, a special species related to amphibian diseases, could be a promising indicator for healthy monitoring during A. davidianus growth. These results could be an important reference for future research on the relationship between the host and microbiota and also provide basic data for the artificial feeding of A. davidianus.


Introduction
The Chinese giant salamander (Andrias davidianus), widely distributed in China, is the largest and most primitive urodele amphibian alive worldwide. It has been listed as a national class II protected species in China and the Appendix I of the Convention on International Trade in Endangered Species of Wild Fauna and Flora (CITES), 2008 (Zhang et al., 2003;Jiang et al., 2015;Turvey et al., 2021). With the breakthrough of artificial breeding technology, the A. davidianus aquaculture industry is rapidly developing. Giant salamander aquaculture could not only meet the consumption needs of modern life but also promote economic revitalization in remote mountains, which is listed as a key project for agricultural industrialization. In China, there were ~2 million A. davidianus artificially reproduced each year (Meng et al., 2014). The breeding modes of this species mainly included entire captive, bionic captive, and primordial ecological breeding (Liang, 2007). Under artificial breeding conditions, A. davidianus generally could grow to the adult stage in 2-3 years. However, the morbidity of infectious diseases is increasing due to intensive culture and artificial changes in the living environment. Infectious diseases mainly included bacterial, viral, parasitic, and other diseases, which could cause heavy economic losses (Jiang et al., 2015;Gui et al., 2018). At present, it is difficult to evaluate whether A. davidianus grows healthily and to choose reasonable measures to prevent relative diseases at the early stage. Hence, it is important to understand the health status and select an accurate monitoring index to guide the breeding of A. davidianus.
The gut microbiota plays important roles in host metabolism and immunity, making it an effective indicator for monitoring the response of animal organisms to environmental changes (Tilg and Kaser, 2011). The gut microbes are diverse, mainly live in the second half of the digestive tract, and consist of nearly 200 common species and about 1,000 uncommon species (Ley et al., 2008). Factors such as host diets, genetic background, and immune status could influence microbiota composition (Turnbaugh et al., 2009;Benson et al., 2010). The population and abundance of gut microbes in different habitats vary greatly, and the status also is reflected in the giant salamander under different temperatures and ages Zhu et al., 2021). There is a close relationship between intestinal flora and host health , which mainly reflected that microbes are mainly colonized in the intestinal tract and play key roles in the digestive system (Gerritsen et al., 2011). The complex microbes constitute a microbial community, and the balance contributes to maintaining the host's gut function, including energy uptake, metabolite production, immune system development, and gastrointestinal diseases (He, 2012).
According to statistics, <1% of microorganisms in nature could be obtained by the traditional media method in vitro. Highthroughput sequencing technology is widely used to detect the gut microecosystem considering the complex and culture-independent feature (Turnbaugh et al., 2010;Siezen and Kleerebezem, 2011). However, there are insufficient profiles and changes of gut microbiota in the Chinese giant salamander, which results in a lack of comprehensive understanding of the relationship between the microbiome and host considering diverse species.
This study mainly focused on the monitoring indices of A. davidianus during the breeding process. Intestinal microbes of healthy A. davidianus in different growth periods were collected for 16S rRNA sequencing to understand the changes in intestinal microflora from larvae to adult stages. The results contributed to establishing a gut microbiota database for A. davidianus and provided much information for the construction of an accurate monitoring system and standardized artificial breeding of A. davidianus.

Animals breeding and sample collection
A total of 80 healthy A. davidianus were cultured in the Experimental Animal Center at Chongqing University of Arts and Sciences. The growth data of A. davidianus, such as body length and weight, were recorded during whole growth periods. According to the growth characteristics, the guts of A. davidianus were collected from five growth stages, which include tadpole stage (ADT), gills internalization stage (ADG), 1 year age (ADY), 2 year age (ADE), and 3 year age (ADS). Three A. davidianus were randomly selected for euthanasia in each group, and the gut was excised from the abdominal cavity. The separated guts were transferred to a sterilized kraft paper and knotted with cotton rope to decrease the crosspollution in the different intestinal segments. The tissue samples were then immersed in a 4% paraformaldehyde solution for fixation. The contents at the end of the rectum were immediately collected and stored at −80°C until further high-throughput sequencing analysis.

Histomorphological observation of guts
The tissue samples of the small intestine, large intestine, and rectum collected at different growth stages were sectionized using the conventional paraffin ultra-thin sectioning method (Schmitt et al., 2019). This process mainly includes embedding, sectioning, and HE staining steps. After completion, the morphological changes of intestinal tissue at each stage were observed with a microscope and photographed using a microscopic imaging system.

DNA extraction
To analyze the composition of the bacterial community, genomic DNA was extracted from the aforementioned gut content samples using the CTAB method (Wilson, 2001). DNA extraction operation was quickly performed after these gut samples were fully mixed. Agarose gel (0.8% w/v) electrophoresis was then performed to evaluate the purity and concentration of the extracted DNA. According to the concentration detected, an appropriate amount of DNA was taken and diluted to 1 ng/μL with sterile water.

16S rRNA amplification and sequencing
Using diluted genomic DNA as a template, specific primers with barcodes were designed according to the selection of sequencing regions. The 16S rDNA target region (V3/V4) was amplified by PCR with primers 338F: 5′-ACTCCTACGGGAGGCAGCA and 806R: GGACTACHV GGGTWTCTAAT-3′. More detailed parameters of the PCR reaction were performed as described previously (Song et al., 2019). PCR products were detected by electrophoresis with agarose gel (2% w/v). The PCR products were purified by magnetic beads, quantified by a microplate reader, and then mixed in equal amounts according to concentration. After full mixing, PCR products were detected by electrophoresis with agarose gel (2% w/v). DNA fragments from the agarose gels were recovered using a QIAquick Gel Extraction Kit. The purified PCR products were used for constructing the sequencing library using TruSeq ® DNA PCR-Free Sample Preparation Kit. Prior to sequencing, the sequencing libraries were subjected to quantification using Qubit and Q-PCR. The qualified libraries were subjected to high-throughput sequencing using an Illumina NovaSeq 6000.

Bioinformatics and statistical analysis
The paired-end sequences from high-throughput sequencing were assigned to the corresponding samples according to the primer and barcode information. After amputation of barcode and primer sequences, the reads of each sample were spliced to obtain raw tags using FLASH (version 1.2.7; Magoč and Salzberg, 2011) and were then strictly filtrated to obtain clean tags (Bokulich et al., 2013). The quality of raw tags was evaluated using Quantitative Insights into Microbial Ecology (QIIME) software (version 1.9.1; Caporaso et al., 2010). For raw tag interception, the threshold value of continuous low-quality bases is 19, and the base length is 3. These tags would be filtered when the base length is <75% of whole tags with continuous high quality. Finally, the effective tags were obtained when raw tags, such as ambiguous bases, chimeras, and mismatched primers in the reads, were filtered by initial quality screening. The obtained high-quality effective tags for all samples were clustered in operational taxonomic units (OTUs) based on 97% identity using the UPARSE algorithm (version 7.0.1001.). Prior to homogenization, the representative sequences for OTUs were annotated and classified based on the Mothur method and SSU rRNA database and were phylogenetically analyzed by MUSCLE (version 3.8.31; Quast et al., 2013). Alpha diversity analysis was performed using QIIME (version 1.9.1) to calculate the values of observed-OTUs, Chao1, Shannon, Simpson, and so on. Beta diversity was used to identify similarities and differences between different samples using the same QIIME software. In addition, rarefaction curves were used to evaluate the rationality of sequencing depth. Linear discriminant analysis effect size (LEfSe) was generated to identify significantly differential biomarkers among groups. Based on the microbiota composition, the functional enrichment of the KEGG pathway was further predicted by PICRUSt2 (version 2.3.0-b; Langille et al., 2013). R software (version 2.15.3) was applied to statistical analysis and plotting. The criterion of significance was conducted at p-values of <0.05, and the data were expressed as means ± SD.
The histological changes of the A. davidianus gut at different growth periods were observed using a microscope after staining the tissue section with hematoxylin-eosin (HE). In the small intestine, both muscular thickness and villi length gradually increased with the increase in age, as shown in ADY, ADE, and ADS stages ( Figures 1A-C). Lymphoid follicular accumulation began to appear in the submucosa in the ADS stage ( Figure 1D).
Similar changes were observed in the large intestine in the three stages. The number of goblet cells gradually increased in the villous epithelium with the growth of the giant salamander (Figures 2A-C). Gland-like structures were seen in the lamina propria of villi in the ADS sample ( Figure 2D). In addition, there were a large number of lymphocytes aggregated in this stage.
For rectum samples, lymph nodes were observed in the submucosa during the ADG period ( Figure 3A). Partial epithelial cell necrosis occurred in the villous epithelium in the ADY stage ( Figure 3B). The intestinal adenoid structure was observed in the villi during the ADE stage ( Figure 3C). The entire intestinal wall thickens during ADS, and the villi become shorter compared to those in other periods ( Figure 3D).

Quality assessment and OTU classification of intestinal microbiota
We initially performed a quality screening for high-throughput sequencing data of intestinal microbiota to eliminate erroneous and questionable sequences, which contributed to verifying sequence reliability. A total of 2,824,829 high-quality reads were produced from the data with an average of 62,774 reads per sample (40,489) and an average length of 252-256 bp ( Table 1).
The qualified reads were composed of 5,772, 1,033, 1,066, 479, and 393 OTUs in ADT, ADG, ADY, ADE, and ADS based on 97% nucleotide-sequence identity, respectively ( Figure 4A). The curves of rarefaction and rank abundance per sample were relatively flat and displayed saturate tendency, which suggested that the depth and evenness of sequences meet the requirements for sequencing and further analysis ( Figures 4B-D).

Analysis of microbial community diversity
The alpha diversity of gut microbiota samples showed the goods coverage estimates varied from 99.3% to 99.9% for all samples, which indicated excellent coverage ( Table 2). The average Chao1 indices for experimental groups ADT, ADG, ADY, ADE, and ADS were 1594. 56, 488.91, 512.81, 354.01, and 262.45, and the corresponding ACE indices were 1589. 90, 492.12, 517.26, 350.47, and 264.35, respectively. Moreover, the averages of Shannon indices for these five groups were 7.43, 5.60, 5.55, 4.58, and 3.03, respectively. The Chao1, ACE, Shannon, and Simpson indices for these five groups displayed gradually downward trends, which indicated that the abundance and diversity of the intestinal microbial community reduced as growth. Remarkably, the three diversity indices (ACE, Chao1, and Shannon) of the initial ADT group were much higher than those of other groups. In contrast, significant differences in gut microbiota abundance and diversity were observed between the ADT and other groups. The α-diversity indices revealed a significant difference in the diversity and richness of gut microbiota among different growth groups. Both the weighted and the unweighted principal coordinate analysis (PCoA) plots revealed that the samples in most groups were clustered separately except for some similarity between ADG and ADY, which Frontiers in Microbiology 04 frontiersin.org indicated that the differences existed in gut microbiota for most of the comparative samples ( Figure 5).
The top 10 genera were Cetobacterium, Lactobacillus, Candidatus Amphibiichlamydia, Hydrogenoanaerobacterium, Akkermansia, Methylobacterium-Methylorubrum, Bacteroides, Desulfovibrio, Bilophila, and Parabacteroides, which accounted for 17.53%-57.96% of the taxonomic groups identified ( Figure 6B). The genus Cetobacterium had a significantly higher abundance than other genera, which comprised 3.93%, 19.49%, 24.13%, 13.72%, and 57.52% of the overall bacterial composition in ADT, ADG, ADY, ADE, and ADS groups, respectively. The levels of this genus in the ADS group were significantly higher than in the other groups. Bacteroides was the second most abundant at 2.90%, 3.59%, 6.13%, 2.02%, and 0.20% for ADT, ADG, ADY, ADE, and ADS groups, respectively. The other dominant genera were Hydrogenoanaerobacterium and Candidatus Amphibiichlamydia, which represented 5.57% (ADG) and 5.41% (ADY) of the overall bacterial composition, respectively. In the ADS group, the most numerous genus was Cetobacterium at 57.52%, whereas other genera except these top 10 were observed to be predominant for other groups at 51.77%-82.47% of the overall bacterial composition, respectively. The relative abundance of these bacteria was also displayed through a heatmap produced by clustering analysis. The distribution of bacterial genera in each sample could also be observed in the heatmap (Supplementary Figure S1).

Microbial profile and core microbiota of intestinal microflora
Analysis of similarities (ANOSIM) was used to evaluate whether differences between groups (two or more groups) were significantly greater than those within groups (Table 3). The results showed Frontiers in Microbiology 05 frontiersin.org significant and remarkable differences among different groups (R > 0, p = 0.001), except the ADG-ADY (p = 0.127). Although there was no significant difference, the R-values for the ADG-ADY comparison were greater than zero (0.106), indicating potential differences between the ADG and ADY groups.
To understand the ecosystem of five samples, linear discriminant analysis effect size (LEfSe) analysis was used to uncover the complex system of microbial communities. Meanwhile, linear discriminant analysis was used to identify the differences in all samples (LDA threshold is 4; Figure 7). The results illustrated that there were 43 bacterial clades, consisting of three classes, 15 orders, and 18 families, which were crucial bacterial branches distinguishing giant salamander samples. According to Figure 7A, Lactobacillus (4.4) and Sediminibacterium (4.1) were significantly enriched in ADT. Hydrogenoanaerobacterium (4.5), Akkermansia (4.4), Parabacteroides (4.2), and Bilophila (4.2) were significantly enriched bacterium in ADG. In ADY, three realms showed significant enrichment, which were Candidatus Amphibiichlamydia (4.5), Bacteroides (4.5), and Desulfovibrio (4.4), respectively. Hafnia Obesumbacterium (4.5) and Alistipes (4.0) were significant in ADE. More bacteria had significant abundance in ADS, such as Cetobacterium (5.4) and Clostridium-sensu-stricto-1 (4.3). In the cladogram, these circles represented different taxonomic levels from phylum to genus, and each small circle represented a classification at that level ( Figure 7B). The relative abundance of microbes is proportional to the diameter size of the circle. The significant differences were marked with different colors consistent with corresponding levels except for the yellow which represented no significant difference. These differential microbial groups that play important roles were visually displayed at different levels in the cladogram.

Function prediction
The Phylogenetic Investigation of Communities by Reconstruction of Unobserved States (PICRUSt) bioinformatics software package was used to predict the metagenomic function of marker genes based on the KEGG database. As shown in Figure 8, the genes obtained from 16S rRNA sequencing were mainly enriched into 41 KEGG pathways and were classified into seven categories. There were 4, 3, 4, 6, 12, 8, and 4 pathways involved in cellular processes, environmental information processing, genetic information processing, human diseases, metabolism, organismal systems, and unclassified categories, respectively. The membrane transport pathway had the most annotated genes including 8,313,360 genes, which belong to environmental information processing. Moreover, the metabolism category had the most pathways and relatively abundant genes, such as carbohydrate metabolism, amino acid metabolism, and energy metabolism had 73,944,375, 66,615,423, and 40,419,909 genes, respectively. Other enriched pathways, such as replication and repair (55,545,160) and poorly characterized (35,422,177), belonged to genetic information processing and unclassified categories, respectively.

Discussion
The intestinal tract is the most important digestive and absorption organ of animals. A large number of bacteria were colonized in the intestines and played indispensable roles in maintaining the overall health of the host (Singh et al., 2017). In our histomorphological results, numerous wrinkles in the lining and gradual increases in villi length were observed in A. davidianus intestine during growth periods, which contributed to the colonization of microorganisms. The intestinal floras were different among different individuals or the same one in different conditions, which were related to ages, diets, health status, and ranging profile (Korpela et al., 2021;Liu et al., 2021;Spencer et al., 2021;Sztandarski et al., 2022). High-throughput sequencing was used to investigate the gut and lung prokaryotic community profiles of adult Chinese giant salamanders at age 3 . In addition, the previous gut microbiota report of A. davidianus provided much microbial information related to age changes from age 1 to 4 . However, studies on earlier ages, histomorphological observation, developmental effects, and disease-related pathogen identification are still lacking. These intestinal microbes are comprised of both beneficial and harmful members. The balance maintenance of intestinal flora and an increase in probiotics proportion could effectively keep the host healthy (Aziz et al., 2013). Therefore, the acquisition of composition and abundance of intestinal flora is well helpful to disease prevention and treatment (Feng et al., 2018).
Gut microbiota plays an important role in the growth and development of host animals with huge abundance and complex structures. We found that the alpha diversity indices, such as OTU number, Shannon, Simpson, Chao1, and ACE, were decreased as A. davidianus grew, which indicated the species diversity and abundance of intestinal flora reduced. These values in the ADT stage were significantly higher than in other periods, which may be deduced that the number of beneficial bacteria increased while conditional pathogens were reduced. The phenomena of lymphoid follicles accumulation and adenoid structure that existed in submucosa first observed in ADG also contributed to the elimination of harmful microbes in guts. A previous study also revealed that phylogeny plays the most important role in the formation of microbial communities, rather than food and environment (Bai et al., 2021). Previous studies verified that aging is a multifactorial process and would influence many principal physiological systems, including the gastrointestinal system (Mabbott et al., 2015). In addition, histomorphological changes in the gut at different ages or stages were found in mice (Yang et al., 2019). In humans, both the composition and stability of gut microbiota were reported to change with age (Li et al., 2016). However, the relationship between gut microbiota and gut histomorphological needs further studies.
The proportion of Fusobacteriota in the ADS group was significantly higher than in other groups, and there was a gradually increased tendency during growth periods except for ADE. More detailed analysis showed that Cetobacterium was the dominant genus in this phylum. On the contrary, the levels of Bacteroidota, unidentified_Bacteria, and Kapabacteria displayed a gradual decline. Fusobacteriota, Firmicutes, and Actinobacteriota are the dominant intestinal phyla for all animals. These similarities suggested that these  microorganisms are important participants of host functions, such as normal digestion, absorption, and immune responses. Cetobacterium, an anaerobic bacteria belonging to the core microbiota of fish gut, was also the most dominant genus in A. davidianus. Lactobacillus, known as beneficial bacteria, was common in the gastrointestinal tract of most aquatic animals (Liu et al., 2021). They could convert large amounts of hexose substrates into pyruvate and then generate the final lactate via the glycolytic pathway (Dempsey and Corr, 2022). More importantly, the abundance of Candidatus Amphibiichlamydia genus in early growth periods (ADT, ADG, and ADY) was much higher than   in later phases (ADE and ADS). We speculated that it may be the result of A. davidianus development, as these salamanders develop, their small intestine increases in complexity, and the wrinkles and partmentalization may play key roles in the abundance changes of Candidatus Amphibiichlamydia. This species has attracted much attention as a special pathogen for amphibian diseases (Martel et al., 2013), which may be a promising potential marker for the prediction of A. davidianus diseases. Although none of the A. davidianus tadpoles showed signs of clinical disease, more experiments related to pathogenic conditions need to be performed considering a high prevalence of 71% in bullfrogs (Lithobates catesbeianus).
In beta diversity, the result of PCoA showed that there was a very similar species composition between ADG and ADY, which indicated a closer relationship in the two samples compared to others. This result was also verified by statistical analysis of ANOSIM and LEfSe, further indicating the differences of intergroup were greater than intragroup. These results were consistent with previous highthroughput sequencing results in zebrafish gut microbiota (Xiao et al., 2021). Host development overwhelmed environmental effects in governing fish gut microbial community succession from larvae to adult fish stages due to host genetics, immunology, and gut nutrient niches. This study is another example of reduced abundance and diversity of the intestinal microbial community as growth.
Maintaining a healthy gut is key to disease prevention in animals. Metabolic activities of microorganisms would generate many important nutrients, such as short-chain fatty acids, vitamins, and amino acids, which would affect host health. The succinate and secondary bile acids produced by Parabacteroides distasonis played key roles in the modulation of host metabolism, and disturbance of intestinal flora was closely related to the occurrence of obesity, diabetes, and hyperlipidemia . The enriched KEGG pathways of annotated genes were also oriented to multiple functions of intestinal microorganisms in the giant salamander. Different growth and developed stages could have differential intestinal structures, which may influence intestine flora. The intestinal flora, in turn, also could affect the intestinal development and immune function of A. davidianus. Keeping the balance of intestinal flora would effectively help the host maintain health status.

Conclusion
This study investigated the changes in the intestinal microbial community in A. davidianus from larvae to adult stages. The results revealed that the diversity and abundance of gut microbiota had significant alterations that were characterized by declining levels of growth. Most of the top 10 phyla and genera had significant differences among different groups, except for a similar microbial community in ADG-ADY. These genes annotated in intestinal microbes were mainly enriched into KEGG pathways including cellular processes and environmental information processing, which played important roles in growth metabolism, nutrient absorption, and immune regulation. In addition, Candidatus Amphibiichlamydia, a special species for amphibians, was a promising potential indicator of gut microbiota stability. These findings expanded our current understanding of the succession of gut microbiota across A. davidianus growth and also provided new insights into the breeding monitoring of other aquatic animals.

Data availability statement
The data presented in the study are deposited in the NCBI repository, accession number PRJNA903076. The data is publicly available with accession number of PRJNA903076 and website: https://www.ncbi.nlm.nih.gov/bioproject/PRJNA903076. Heatmap of the KEGG pathways annotated by PICRUSt.