Detection and comparative analysis of cutaneous bacterial communities of farmed and wild Rana dybowskii (Amphibia: Anura)

Abstract Rana dybowskii Gunther, 1876 is a dominant amphibian species in northeast China. In order to understand the composition and structure of the cutaneous bacterial communities of farmed and wild R. dybowskii, two experimental groups (farmed and wild) were investigated in this study. Following DNA extraction, the V3-V4 hypervariable region of the 16S rRNA gene was targeted and analyzed by high-throughput sequencing. The cutaneous bacterial community diversity was investigated and compared through analysis of alpha and beta diversity. A total of 524,852 valid sequences and 1,023 operational taxonomic units (OTUs) were obtained from these two experimental groups. The number of shared OTUs was 333, while there were 603 unique OTUs in the farmed group and 87 unique OTUs in the wild group. The Chao, ACE and Shannon indices of the farmed group were significantly higher than those of the wild group (p < 0.05). The composition and abundance of the dominant bacteria at the phylum and genus levels were different. The dominant phyla in the farmed group were Firmicutes, Proteobacteria, Actinobacteria, Chloroflexi, Fusobacteria, Bacteroidetes, and Cyanobacteria. The dominant phyla of the wild group were Proteobacteria, Bacteroidetes, Firmicutes, and Actinobacteria. Furthermore, the alpha and beta diversities of the cutaneous bacterial communities of farmed and wild R. dybowskii exhibited significant differences. This study provides a theoretical basis for comprehensively understanding the composition and abundance of bacteria on the skin of farmed and wild R. dybowskii, which helps to develop its breeding industry and gain economic benefits.


Introduction
In recent decades, amphibians have become the most endangered group of vertebrates with their numbers sharply decreasing throughout the world. The major causes of the decline in amphibians are global warming, habitat loss, diseases and various human factors, including extermination and use of pesticides. Almost onethird of amphibian species are threatened with extinction worldwide, and the conservation of amphibians is necessary (Ochoa-Ochoa et al. 2009). Chytridiomycosis, which is caused by a fungus named Batrachochytrium dendrobatidis (Bd), has received increasing attention among the factors that are causing the decline of amphibians. Chytridiomycosis can affect the function of amphibian skin and cause a large number of amphibian deaths in a short time (Ryan et al. 2008;Van Rooij et al. 2015). Chemicals released by Bd can cause host pathology, even in the absence of infection (Mcmahon et al. 2013).
Rana dybowskii is taxonomically classified as Amphibia, Anura, Ranidae, Rana. Rana dybowskii is one of the dominant amphibian species in northeast China and is mainly found in the northeastern provinces of Heilongjiang, Jilin, Liaoning and Inner Mongolia. Rana dybowskii has been classified as a near-threatened species on the Red List of China's Vertebrates. In recent years, a large number of R. dybowskii have been killed, as their oviducts are used in expensive traditional Chinese medicine, thus, they have a very high economic value. In order to protect wild R. dybowskii while still gaining economic benefits, artificial breeding of R. dybowskii is gradually emerging as an industry. Artificial breeding is an important method for preventing the extinction of wild animals; however, there are many difficulties with this process, such as a reduced survival rate and susceptibility to disease. One possible cause of these difficulties is the change of resident microorganisms. It has been experimentally shown that the cutaneous bacterial communities of wild species change after artificial breeding ).
An animal's skin is the first level of protection against potential environmental pathogens. The function of the skin is more than a physical obstacle. Over the years, skin has been proven to have a critical immunological role, not only being a passive protective barrier but also a highly sophisticated network of effector cells and molecular mediators known as the "skin immune system" (SIS). Some microorganisms on the surface of an animal's skin can enhance host immunity, whereas some are opportunistic pathogens that are harmful to the body. Microorganisms that inhabit the skin of vertebrates dynamically change throughout the host's life cycle. Metabolites produce by these microorganisms combine with active compounds secreted by the host's glands to form a chemical defense (Sabino-Pinto et al. 2016). For amphibians, the skin provides a protective barrier against pathogens and plays an important role in the immune response of the host. In addition, the amphibian's skin has other functions. It is rich in glands (mainly mucous and granular glands) that secrete different bioactive compounds, as alkaloids and antimicrobial peptides, to enhance defense (Huang et al. 2016).
Amphibians have a large number of bacteria on their skin. In addition to participating in the physiological processes of the host, symbiotic bacteria also play an important role in disease resistance and can directly defend against pathogen invasion (Familiar López et al. 2017). Symbiotic bacteria can increase protection by competing with pathogens for resources, supporting the host immune system, or through secretions (Bates et al. 2018). Studies have identified Janthinobacterium lividum and Serratia marcescens on the skin of Atelopus zeteki, both of which secrete antibacterial active substances that increase the host's ability to resist pathogens (Woodhams et al. 2018). The study of Malagasy frogs found that skin-associated bacteria mainly belong to Actinobacteria, followed by Proteobacteria and Bacteroidetes, with many able to inhibit Bd (Bletz et al. 2017). The living environment will, to some extent, affect the bacterial community structure in amphibians (Harris et al. 2006Becker et al. 2009). Medina et al. (2017) demonstrated that the bacterial community on the skin of amphibians is altered in direct relation to the environment (Medina et al. 2017). Becker et al. (2014) found that the species diversity, phylogenetic diversity and the structure of cutaneous bacterial communities were significantly different between captive and wild Panamanian golden frogs . Furthermore, changes in diet have a great impact on the bacterial community of amphibians' skin (Antwis et al. 2014).
Parasitic bacteria on the skin are closely related to the health of R. dybowskii; however, some bacteria may play a role in resistance to external pathogens. Therefore, it is important to know the structure of the bacterial community on the skin of R. dybowskii. The community of parasitic bacteria on the skin of R. dybowskii is subject to many environmental factors such as humidity, temperature, and pH. As farmed and wild R. dybowskii live in different environments, the skin parasitic bacteria change accordingly. In addition, the lack of natural environmental bacterial reservoirs can also alter the composition of the bacterial community associated with the host. In order to get a better understanding of differences in bacterial communities on the skin of R. dybowskii from different environments, wild and farmed R. dybowskii in northeastern China were sampled. In this paper, the 16S rRNA gene of cutaneous bacteria was subjected to high-throughput sequencing to determine the composition and abundance of cutaneous bacterial communities on farmed and wild R. dybowskii. This study focused on alpha diversity and comparative analysis between the two groups, whereby the community composition and relative abundances were measured to understand the differences in these cutaneous bacteria communities.

Sample collection
Two groups of R. dybowskii were collected from northeastern China. The wild individuals (n = 7; weight, 4.89 ± 0.24 g) were collected from a natural R. dybowskii habitat in Luobei County (47°65′24″N, 130°46′24″E; 98 m alt) in June 2017. The farmed individuals (n = 7; weight, 5.14 ± 0.35 g) were collected from a farm in Huanan County (46°44′54″N, 130°6 9ʹ32″E; 80 m alt) in July 2017. The individuals in both groups were sampled on the same day and at the same location. They were sacrificed immediately upon arriving in the laboratory to prevent changes to the cutaneous bacterial community due to the laboratory environment.
Cultivating conditions: The young frogs used for this experiment were robust second instar R. dybowskii, which end hibernation in early May. The culture density was forty individuals per m 2 and the R. dybowskii was fed mealworms. The frogs were housed in a greenhouse (40 m length × 8 m width) with sparse and low vegetation planted on the ground. Water spray equipment and shading nets were installed in the greenhouse. The ground humidity was 25%-35% and there was usually no accumulation of water. The soil temperature was basically the same as the temperature in the greenhouse, but there was a certain degree of hysteresis. Iodophors were sprayed every 3-5 days in the greenhouse for disinfection.
For this study, seven individuals per group were randomly collected from the same area to reduce pseudo-replicates. Exploratory analyses of our dataset did not yield differences in bacterial community composition or diversity between sexes Kueneman et al. 2014;Rebollar et al. 2016); therefore, data from male and female individuals were combined prior to analysis. All individuals were sacrificed and the skin of the back, venter, and limbs were peeled off and collected. The skin samples were loaded into 5 mL sterile tubes and cryopreserved at −80°C.
Ethical approval: Before sample collection, all animal protocols were approved by the Institutional Animal Care and Use Committee of the Northeast Agricultural University, China. All experiments were performed in accordance with the approved guidelines and regulations. IACUC#2015-035.

Illumina MiSeq sequencing
Purified amplicons were pooled in equimolar concentrations and paired-end sequenced (2 × 300) on the Illumina MiSeq platform (Illumina, San Diego, USA) according to the standard protocols.
Raw fastq files were demultiplexed, qualityfiltered using Trimmomatic and merged with FLASH as follows: (i) reads were truncated at any site that received an average quality score < 20 over a 50 base pairs (bp) sliding window; (ii) primers were exactly matched allowing 2 nucleotide mismatching and reads containing ambiguous bases were removed; and (iii) sequences overlapping by < 10 bp were merged accordingly.

Statistical analysis
MOTHUR (version v.1.30.1 http://www.mothur.org/ wiki/Schloss_SOP#Alpha_diversity) was used to create rarefaction data and alpha diversity indices. Wilcoxon rank-sum tests were used to analyze the difference between alpha diversity indices and the relative abundance of phyla and genera in these two groups. Bray Curtis was used for the analysis of similarities (ANOSIM). Weighted and unweighted Unifrac distance matrices were used to calculate beta diversity and were visualized by principal coordinates analysis (PCoA) (Rebollar et al. 2016). The linear discriminant analysis (LDA) effect size (LEfSe) method was used to determine the bacterial taxa that most likely explained differences between the two groups. Following the precedent of previous studies, phyla and genera with LDA scores > 3.5 were considered informative (Lokesh & Kiron 2016).

Results
In this study, the V3-V4 hypervariable region of the bacterial 16S rRNA gene was sequenced to determine the composition and abundance of cutaneous R. dybowskii bacterial communities using highthroughput sequencing.

Sample sequencing data
A total of 14 samples were sequenced for the two study groups, resulting in 524,852 valid sequences (average sequence length = 446) and 233,963,272 bases after sequence filtration and double-stitched splicing. After clustering to 97% similarity, 1,023 total OTUs were obtained, which belonged to 700 species, 471 genera, 218 families, 113 orders, 52 classes, and 26 phyla. The rarefaction curve tended to be flat (Supplementary material 1), indicating that the samples were sequenced to a reasonable amount and reflected the vast majority of the species in the samples.

Alpha diversity
With respect to microbial community ecology, alpha diversity can reflect the abundance and diversity of microbial communities, including a series of statistical indices to estimate species abundance and diversity within the community. Sobs, which represents the actual number of observed OTUs in each sample, confirms that the species richness was higher in the farmed group (Table I). The coverage index estimates how well the sequenced community represents the actual microbial community of the sample. The experimental results show that the coverage index of the samples was greater than 99% (Table I), indicating that the sequencing results were credible. The Chao and ACE indices, which estimate the total number of species and reflect bacterial abundance, were significantly higher (p < 0.01; Chao: p = 0.0022; ACE: p = 0.0022) for the farmed group than for the wild group (farmed group: Chao = 672.33 ± 52.32, ACE = 661.65 ± 50.57; wild group: Chao = 253.40 ± 31.14, ACE = 247.66 ± 25.74). The Shannon (H) and Simpson (D) indices reflect the diversity of the bacterial community. The higher the Shannon index, the higher the community diversity; whereas the higher the Simpson index, the lower the community diversity. The diversity of the farmed group was significantly higher (p < 0.01; H: p = 0.0022; D: p = 0.0049) than that of the wild group (farmed group: H = 3.5540 ± 0.8557, D = 0.1434 ± 0.1168; wild group: H = 1.5773 ± 0.3887, D = 0.3872 ± 0.0961).
Of the 15 most dominant bacteria, 14 phyla exhibited significant differences in the relative abundances between the wild and farmed R. dybowskii (Wilcoxon rank-sum test, two-tailed test, FDR correction, p < 0.05; Figure 1(a)), 13 that were extremely significant (p < 0.01). At the genus level, there were significant differences in the relative abundances between the two groups for 12 of the genera  Figure 1(b)), 11 of which were extremely significant (p < 0.01).

Shared OTUs and beta diversity analysis
To determine the differences in bacterial communities between the farmed and wild groups, a heatmap was drawn using the relative abundances of the top 30 species (Figure 2). The heatmap shows that the relative abundance of most bacteria was different, although there are similarities between some of the samples in the two groups. There were 333 shared OTUs within the farmed and the wild groups. Although these OTUs were present in Figure 1. Wilcoxon rank-sum test bar plot at the phylum (a) and genus (b) levels. The ordinate indicates the species name at different classification levels, the abscissa indicates the percentage value of a species abundance of the sample, and the different colors indicate different groupings. *0.01 < p ≤ 0.05, **0.001 < p ≤ 0.01, *** p ≤ 0.001. most samples, their abundance was different between the groups. In contrast, there were 603 unique OTUs in the farmed group and 87 unique OTUs in the wild group. The shared OTUs accounted for 35.58% of the OTUs in the farmed group and 79.29% in the wild group. There were 13 shared OTUs with a high abundance (> 1% of the total OTUs), belonging to Firmicutes (n = 4 OTUs), Bacteroidetes (n = 2), Proteobacteria (n = 3), Actinobacteria (n = 3) and Fusobacteria (n = 1).
The beta diversity using both unweighted (Figure 3(a)) and weighted (Figure 3(b)) UniFrac analyses showed that the communities of the farmed and wild groups were phylogenetically distant from each other. The weighted method takes into account not only the presence/absence of each OTU but also the abundance. Clustering was performed by principal coordinates analysis (PCoA) of the UniFrac distance matrix, with the farmed and the wild groups separating along with the principal coordinates. Furthermore, ANOSIM showed that the difference between the groups was greater than the difference within each group and that the difference was significant (R = 0.771, p < 0.05).

Linear discriminant analysis (LDA) effect size
The LEfSe method (Rebollar et al. 2016) was used to identify bacteria that best explained the differences between the two groups ( Figure 4). LEfSe identified eight distinct phyla (Firmicutes, Actinobacteria, Chloroflexi, Fusobacteria, Cyanobacteria, Deinococcus__Thermus, Proteobacteria, and Bacteroidetes) with significant differences in abundance. Among them, Bacteroidetes and Proteobacteria were associated with the wild group and the rest were associated with the farmed group. The largest number of species with significant differences was Actinobacteria. Alphaproteobacteria, which belongs to Proteobacteria, was a significantly different class of the farmed group. Simultaneously, Proteobacteria was the most abundant phyla in the wild group whose dominance was driven by Betaproteobacteria and Gammaproteobacteria, as demonstrated with higher LDA scores. All significant species identified in Actinobacteria and Firmicutes belonged to the farmed group. The number of significantly different bacteria at the genus to the phylum level in the farmed group was far more than that in the wild group.

Discussion
In recent years, a great deal of attention has been given to cutaneous bacterial communities in animals because it is thought to be a key factor affecting the animal status, including health, growth, and disease (Mckenzie et al. 2012;Walke et al. 2014;Avena et al. 2016). The study of microorganisms has become relatively facilitated due to the latest advances in high-throughput sequencing, enabling culture-independent analysis. Thus, highthroughput sequencing was used in this study to explore the composition and structure of cutaneous bacterial communities of farmed and wild R. dybowskii. The results show that the bacterial community structures between the two groups living in different environments were significantly different.
The dominant phyla of the farmed group were Firmicutes, Proteobacteria, Actinobacteria, Chloroflexi, Fusobacteria, Bacteroidetes, and Cyanobacteria. The dominant phyla of the wild group were Proteobacteria, Bacteroidetes, Firmicutes, and Actinobacteria. When comparing the two groups, there was a similarity in bacterial composition at the phylum level, but the abundances were quite different. The abundance of Firmicutes and Actinobacteria were significantly higher (p < 0.05) in the farmed group compared with the wild group. Actinobacteria have been detected in soil samples from many places (Elbendary et al. 2018) and they are also common bacteria on R. dybowskii skin (Kueneman et al. 2014). They can produce antibiotics (Barka et al. 2015) and some of them can be used to produce enzyme preparations (Le Roes-Hill et al. 2011). The abundance of Proteobacteria in the wild group was Figure 4. Cladogram with differently abundant taxa color-coded by classification (different groups). Yellow nodes represent species with no significant difference; bacteria associated with the different groups are represented using letters. Significant differences in abundance were identified by the non-parametric factorial Kruskal-Wallis (KW) sum-rank test, with an LDA score cut-off of 3.5. higher than in the farmed group. Proteobacteria are more competitive than other microorganisms because of their morphological and physiological diversity (Shin et al. 2015). Additionally, they have the capacity to adapt to harsh conditions of high temperature, and many of them are opportunistic pathogens (Jardine et al. 2017). Bacteroidetes was highly dominant in the wild group, second only to Proteobacteria. Compared to the wild group, the abundance of Bacteroidetes in the farmed group was very low (average relative abundance < 1%). Furthermore, the abundance of Firmicutes in the wild group was lower than that in the farmed group. Studies have shown that Firmicutes can survive under harsh environmental conditions, which may be related to their ability to produce endospores that are resistant to desiccation (Gomez-Montano et al. 2013). Bacillus is often used as a probiotic on farms to prevent disease; therefore, the high abundance of Firmicutes in the farmed group may depend on human factors. Bacteroidetes and Firmicutes are predominant bacteria on the skin of many amphibians (Abarca et al. 2017;Bates et al. 2018). The relative abundances of the phyla Proteobacteria, Bacteroidetes, Firmicutes and Actinobacteria in the two groups were significantly different, both of which belonged to the phylumlevel in LDA.
There were differences at the genus level between the farmed and the wild groups. The relative abundances of Pseudomonas, unclassified_f__Flavobacteriaceae, and Sphaerotilus in the wild group were significantly higher than those in the farmed group (p < 0.05). In contrast, the relative abundances of Staphylococcus, Salinicoccus, Nesterenkonia Citrobacter, Lactococcus, Rhodococcus, Glutamicibacter, norank_o__JG30-KF-CM45 and Enterococcus were significantly higher in the farmed group (p < 0.05). Pseudomonas was abundant in both the wild and farmed groups. Pseudomonas is a common opportunistic pathogen (Meynet et al. 2018) widely found in air, skin, water and soil ecosystems. The genus Pseudomonas is ubiquitous in moist environments and some are pathogenic to animals and humans (Moradali et al. 2017). It is also a resident bacterium on the skin of amphibians and studies have shown that Pseudomonas isolated from Hemidactylium scutatum has an inhibitory effect on Bd in vitro (Harris et al. 2006). The abundance of Ralstonia and Bacillus was similarly low in both groups. The relative abundances of twelve dominant bacteria in the farmed group were less than 1% in the wild group, while four of the dominant genera in the wild group had an average relative abundance of less than 1% in the farmed group. The experimental results showed that the average relative abundance of genera belonging to R. dybowskii grown in two contrasting environments was quite different.
Most OTUs were shared between the two groups. Nevertheless, the PCoA plot showed little overlap between the samples. Furthermore, the beta diversity between the bacterial communities was significantly different. The 603 unique OTUs observed in the farmed group, which was much higher than that of the wild group, may explain why the Shannon diversity index of the farmed group was higher than the wild group. The ANOSIM showed that the cutaneous bacterial communities were significantly different between the farmed and wild groups.
There are many factors that can affect the cutaneous bacterial community structure of amphibians, such as the living environment, different body parts, host species, temperature, the effects of pathogens, the developmental stage of the host, seasons and Bd infections (Longo et al. 2015;Van Rooij et al. 2015). Many other factors that affect the cutaneous bacterial communities remain to be explored. In this study, the living environment may primarily explain the changes in the data sets. Thus, among the many influencing factors, the cutaneous bacterial community structures may be affected most by the living environment.
Wild R. dybowskii prefers to inhabit vegetated environments with tall trees, shrubs and humid air, such as broad-leaved or conifer-broadleaf forests. Although the R. dybowskii farm is located in the natural distribution area of R. dybowskii, the farm environment is affected by human factors and many other differences compared with the natural environment. Studies have shown that after approximately about eight years of living in captivity, the offspring of captive golden frogs still shared 70% of their microbial communities with wild frogs. Hostrelated microbial communities changed significantly during captive management, but most community members can be retained ). However, differences in cutaneous bacterial community structure, species richness, and diversity between the wild and farmed R, dybowskii were observed in this research. Similar results have been reported in previous studies on other amphibians. For example, the cutaneous bacterial communities were significantly different between wild and captive Cynops pyrrhogaster (newts), with the wild individuals exhibiting a higher alpha diversity (Sabino-Pinto et al. 2016). Conversely, this study found that the alpha diversity of the farmed bacterial community significantly increased (p < 0.05), which is inconsistent with the previous studies. For most animals, wild individuals live in a complex environment, while cultured individuals live in a single environment (Sabino-Pinto et al. 2016). However, the ability of R. dybowskii to adapt to the living environment is not strong. This conclusion can be drawn from the limited distribution range of R. dybowskii (mainly distributed in Northeast China). The R. dybowskii farm needs to be built within the natural habitat of R. dybowskii, so that the natural survival environment of wild and farmed R. dybowskii is similar. At the same time, the farms are disturbed by human factors, which are more complex than the natural environment. This may be the reason for the increased bacterial community diversity of the farmed R. dybowskii. A possible reason for this is the complex farm environment where pollution is more serious. Studies have shown that mice in a farmhouse harbored a significantly more diverse and richer gut microbiota, due to their exposure to unhygienic conditions, compared with mice reared in a specific pathogen-free animal room (Zhou et al. 2016). The cultured R. dybowskii is susceptible to disease (Weng et al. 2016), and the treatment options are limited. If they are infected, their bacterial community will be affected. Furthermore, there is an interaction between surface and pathogenic bacteria, which is also a way for amphibians to resist pathogenic bacteria (Lam et al. 2010). Higher microbial diversity may be caused by infection, because the host may actively recruit beneficial microorganisms as a mechanism to resist infection (Longo et al. 2015). The differences in bacterial communities can also be caused by different diets. Wild R. dybowskii are free to prey upon many species while the food for farmed R. dybowskii is limited mainly to yellow mealworm (Tenebrio molitor). However, the food of wild R. dybowskii comes from their own environment, while the food of farmed R. dybowskii has been exposed to further complexities, such as transportation equipment and farm staff. Subsequently, the bacteria obtained from this environment may be the reason for the higher diversity of the bacterial community of farmed R. dybowskii .

Conclusion
This work can provide basic data for the study of R. dybowskii. The dominant phyla in the farmed group were Firmicutes, Proteobacteria, Actinobacteria, Chloroflexi, Fusobacteria, Bacteroidetes, and Cyanobacteria. The dominant phyla of the wild group were Proteobacteria, Bacteroidetes, Firmicutes, and Actinobacteria. There were significant differences in the alpha and beta diversities of the cutaneous bacterial communities of the farmed and wild R. dybowskii. During farming, R. dybowskii is susceptible to various diseases, including red-legged syndrome, rotten skin disease, fester disease and curved worm disease. Further studies are needed on the changes in cutaneous bacterial community structures of R. dybowskii after disease and drug application. The bacteria that inhabit the skin of the host have a certain influence on the status of the disease. Thus, an understanding of the cutaneous bacterial community can provide a theoretical basis for the prevention and treatment of R. dybowskii disease.