17β-Estradiol supplementation changes gut microbiota diversity in intact and colorectal cancer-induced ICR male mice

The composition of the gut microbiota is influenced by sex hormones and colorectal cancer (CRC). Previously, we reported that 17β-estradiol (E2) inhibits azoxymethane/dextran sulfate sodium (AOM/DSS)-induced tumorigenesis in male mice. Here, we investigated whether the composition of the gut microbiota is different between male and female, and is regulated by estrogen as a secondary outcome of previous studies. We established four groups of mice based on the sex and estrogen status [ovariectomized (OVX) female and E2-treated male]. Additionally, three groups of males were established by treating them with AOM/DSS, and E2, after subjecting them to AOM/DSS treatment. The mice were sacrificed at 21 weeks old. The composition of the gut microbiota was analyzed using 16S rRNA metagenomics sequencing. We observed a significant increase in the microbial diversity (Chao1 index) in females, males supplemented with E2, and males treated with AOM/DSS/E2 compared with normal males. In normal physiological condition, sex difference and E2 treatment did not affect the ratio of Firmicutes/Bacteroidetes (F/B). However, in AOM/DSS-treated male mice, E2 supplementation showed significantly lower level of the F/B ratio. The ratio of commensal bacteria to opportunistic pathogens was higher in females and E2-treated males compared to normal males and females subjected to OVX. Unexpectedly, this ratio was higher in the AOM/DSS group than that determined in other males and the AOM/DSS/E2 group. Our findings suggest that estrogen alters the gut microbiota in ICR (CrljOri:CD1) mice, particularly AOM/DSS-treated males, by decreasing the F/B ratio and changing Shannon and Simpson index by supply of estrogen. This highlights another possibility that estrogen could cause changes in the gut microbiota, thereby reducing the risk of developing CRC.


Results
Effect of sex and CRC on clustering and phylum compositions following supplementation with estrogen. The experimental groups were distinguished into two categories, which were as follows: sex and OVX as Group 1, and colon cancer in males as Group 2. Group 1 was further divided into 4 subgroups, which were as follows: normal male control mice (n = 5), E2 supplemented male mice (n = 12), normal female control mice (n = 11), and ovariectomized (OVX) female mice (n = 7) (Fig. 1a). Similarly, group 2 was divided into 3 subgroups, which were as follows: normal male control mice (n = 5), AOM/DSS-treated male mice (n = 12), and AOM/DSS-treated male mice supplemented with E2 (n = 12) (Fig. 1b). Female mice were subjected to a sham surgery or ovariectomy at 4-week-old. To induce colitis-associated CRC, 5-week-old male mice were injected intraperitoneally with 10 mg/kg AOM and were given 2.5% DSS in the drinking water for 7 days starting 1 week following the injection of AOM. Estrogen supplemented groups were received a daily intraperitoneal injection of E2 for 7 days, and other groups were received olive oil as a vehicle (Fig. 1). Their feces were freshly collected and immediately frozen in liquid nitrogen.
After genomic DNA preparation from stools, high-throughput 16S rRNA metagenome sequencing was carried out to determine whether differences in the gut microbiome are influenced by sex of the individual and treatment with E2 with or without AOM/DSS treatment in male ICR mice. The rarefaction curve of all 59 samples formed a plateau (Supplementary Fig. 1). Variations occurring in the gut microbial communities among different groups were assessed based on their generalized UniFrac distances. Results of the principal coordinates analysis (PCoA) showed that the bacterial communities could be distinguished based on supplementation with estrogen ( Fig. 2a,  Supplementary Fig. 2a). Except the Male_E2 samples belonging to Group 1, which were representative of the sex and OVX groups, no clear divergence based on the differences in sex or OVX and AOM/DSS or estrogen treatments was observed in samples obtained from members of Group 1 and Group 2, which is indicative of the Scientific RepoRtS | (2020) 10:12283 | https://doi.org/10.1038/s41598-020-69112-w www.nature.com/scientificreports/ development of colon cancer in males, respectively ( Fig. 2a,b, Supplementary Fig. 2a, b). PERMANOVA results demonstrated the beta set-significance between Group 1 (p = 0.001) and between Group 2 (p = 0.017) (Fig. 2a,b). In Group 1, different compositions of the gut microbiota at the levels of both phylum and family taxa were observed (Fig. 2c,e, Supplementary Table 1, 2). At the phylum level, the abundance ratio of Cyanobacteria decreased significantly by estrogen supplementation in the groups involving male members (p < 0.001 based on the Kruskal-Wallis test; p < 0.001 for Male_Con vs. Male_E2) (Fig. 2c, Supplementary Table 1). Furthermore, the abundance ratio of Verrucomicrobia increased significantly in the absence of estrogen (OVX) in the groups with females (p = 0.054 based on the Kruskal-Wallis test; p = 0.007 for Female_con vs. Female_OVX) (Fig. 2c, Supplementary Table 1). Differences in the gut microbiota of members belonging to Group 2 were observed only at the family taxa level (Fig. 2d,f, Supplementary Table 1, 2). Furthermore, the ratio of Firmicutes to Bacteroidetes (F/B ratio) did not differ significantly among the members of Group 1 (Fig. 2g). Interestingly, the F/B ratio was significantly lower in the AOM/DSS group supplemented with E2 than that of their AOM/DSS counterparts in Group 2 (p = 0.061 for the Kruskal-Wallis test, p = 0.050 for Male_AOM/DSS vs Male_AOM/DSS/E2) (Fig. 2h).
The counts and abundance ratios of all operational taxonomic units (OTU) are presented in Supplementary Dataset 1 and 2. All p-values and q-values of the false discovery rate (FDR) obtained by comparing the results of Group 1 and Group 2 are shown in Supplementary Dataset 3. The result of LDA effect size obtained by comparing between Group 1 and between Group 2 are shown in Supplementary Dataset 4. The microbial compositions at the genus levels are presented in Supplementary Fig. 3a-b. Members of the genus Alistipes were the most dominant organisms detected in the male and female control groups, accounting, for 26% and 24% of the total, respectively ( Supplementary Fig. 3a). In the estrogen supplemented male (Male_E2) and estrogen eliminated female (Female_OVX) groups, Bacteroides was the most prominent genus identified, accounting for nearly 21% and 14% of the total, respectively ( Supplementary Fig. 3a). In Group 2, Alistipes decreased in the Male_AOM/DSS group (17%) compared to Male_Con group (26%). Interestingly, the decrement of Alistipes was slightly recovered by estrogen supplementation (Male_AOM/DSS/E2) (20%) in Male_AOM/DSS group ( Supplementary Fig. 3b).
Influence of sex and CRC following estrogen supplementation on the diversity of the gut microbiota. It was observed that, in Group 1, the OTU count increased significantly in both, estrogen supplemented male and female control groups, compared to that observed in the male control group (p = 0.002 for Male_Con vs. Male_E2, p = 0.008 for Male_Con vs. Female_Con) (Fig. 3a, Table 1). Similar to the OTU count, species richness of the gut microbiota, as indicated by the Chao 1 index, was also found to be significantly increased in estrogen supplemented male and female control groups compared to that of the male control (p = 0.002 for Male_Con vs. Male_E2, p = 0.008 for Male_Con vs. Female_Con) (Fig. 3c, Table 1). However, the observed OTU count and Chao 1 index did not vary post ovariectomy compared to those of the female control group (Fig. 3a,c, Table 1). Interestingly, significant alterations in the alpha diversity indices were observed in the OVX group compared to that of the female control (Shannon index, p = 0.013; Simpson index, p = 0.006 for Female_Con vs. Female_OVX) (Fig. 3e,g, Table 1). We observed that, in Group 2, treatment of male mice with AOM/DSS resulted in significant alteration in the OTU count (p = 0.006 for Male_Con vs. Male_AOM/DSS) and values of the Chao 1 index compared to those determined in the male control group (Fig. 3b,d, Table 1). On the other hand, the alpha diversity indices (Shannon and Simpson) was significantly changed upon supplementation with E2 in the AOM/DSS male group (Shannon index, p = 0.033 and Simpson index, p = 0.028 for Male_AOM/ DSS vs Male_AOM/DSS/E2) (Fig. 3f,h, Table 1).  www.nature.com/scientificreports/ Influence of sex and CRC due to estrogen supplementation on the composition of the gut microbiota. We performed the linear discriminant analysis (LDA) effect size (LEfSe) analysis to identify the taxonomic biomarker in the gut microbiota. According to the results of the LEfSe analysis, the proportions of few genera constituting the gut microbiota were higher or lower significantly in both Group 1 and Group 2 ( Fig. 4). In the male control mice, the relative abundances of opportunistic pathogens including Neglecta, Eubac-terium_g23, Ruminococcus, and Eubacterium_g17 (p = 0.022, 0.039, 0.034, and 0.026, respectively) decreased in response to estrogen supplementation (Fig. 4a). In the female control mice, the relative abundances of opportunistic pathogens, including Eubacterium_g8 and Akkermansia (p = 0.022 and 0.007, respectively), increased while that of the commensal bacteria Alistipes (p = 0.013) decreased in response to the estrogen elimination (OVX) (Fig. 4b). In the male control mice, the relative abundance of the opportunistic pathogen Neglecta (p = 0.019) was higher and that of commensal bacteria including KE159571_g and KE159538_g (p = 0.041 and 0.006, respectively) was lower compared to that of the female control mice (Fig. 4c). Interestingly, abundance of the commensal bacteria Butyricimonas (p = 0.047) was lower in the female control mice compared to that observed in the males (Fig. 4c). In addition to these results, the ratio of commensal bacteria to opportunistic pathogens was found to be higher in Male_E2 and Female_Con groups than the Male_Con or Female_OVX groups (Fig. 4f). The relative abundance of the opportunistic pathogen Roseburia (p = 0.009) decreased in the Male_AOM/ DSS group compared to that in the Male_Con group (Fig. 4d). Surprisingly, the abundances of commensal bacteria including Butyrivibrio_g1, Muribaculaceae_uc, and KE159538_g (p = 0.022, 0.045 and 0.020, respectively) were higher in the Male_AOM/DSS group compared to that in the Male_Con group (Fig. 4d). Additionally, the abundances of all of commensal bacteria including PAC000664_g and Phocea (p = 0.007 and 0.001, respectively) and opportunistic pathogens including Pseudoflavonifractor and Neglecta (p = 0.012 and 0.003, respectively) diminished in the Male_AOM/DSS/E2 group compared to that in the Male_AOM/DSS group (Fig. 4e). Unexpectedly, the proportion of commensal bacteria to opportunistic pathogens was higher in the Male_AOM/DSS group than the Male_Con or Male_AOM/DSS/E2 groups (Fig. 4f).
Functional profiles of microbial communities in mice of different sex and with CRC following estrogen supplementation. PICRUSt analysis was conducted to investigate the potential function of the gut microbiota and predict their classifications based on the Kyoto Encyclopedia of Genes and Genomes (KEGG) Ortholog (KO) (www.kegg.jp/kegg/kegg1 .html). As a result, we identified differentially enriched functional biomarkers in Group 1 ( Table 2) and Group 2 ( Table 3). In Group 1, two orthologs were found to be higher in both Male_E2 and Female_Con. groups compared to those of the Male_Con. or Female_OVX groups, which were the ABC-2 type transport system ATP-binding protein (K01990) and TonB-dependent starch-binding outer membrane protein SusC (K21573) ( Table 2). In Group 2, 17 orthologs increased and 8 orthologs decreased in Male_AOM/DSS group compared with those in the Male_Con. or Male_AOM/DSS/E2 groups, respectively (Table 3).

Discussion
A direct interaction between gut microbiota, sex hormones, and the incidence of diseases has been determined based on the data obtained from many studies using animal models, thus suggesting that the potential differences in the composition of the gut microbiota between different sex 34 . The aim of previous studies was to investigate whether tumorigenesis is reduced by estrogen supplementation in AOM/DSS-induced CRC male or OVX female mice model 31,35 . We reported that E2 inhibited the initiation of CRC development by upregulating Nrf2 in AOM/ DSS-treated male ICR mice 31 . Furthermore, AOM/DSS-induced proximal colon cancer after ovariectomy was reduced by E2 supplementation 35 . The stool used in this study was taken from the previous experimental sets 31,35 and the current study was conducted as a secondary outcome. In this study, we investigated whether the composition of the gut microbiota is altered differently between sex and its correlation with the presence of estrogen  Members of Firmicutes and Bacteroidetes were the most common bacterial phyla constituting the human microbiome. Therefore, the F/B ratio is used as a representative index to compare between different microbial communities 36 . Many studies have supported the idea that a low F/B ratio signifies a healthy condition 36,37 . Additionally, a corresponding decrease in the F/B ratio was observed in patients with CRC, which could be an important marker for intestinal dysbiosis 38,39 . A higher F/B ratio at a BMI of 33 was observed in men compared to that in women 40 . Contrarily, a significantly lower F/B ratio was observed in men than that in women in the group with BMI > 33 40 . Furthermore, evidence from many studies has suggested that the F/B ratio was higher in normal range women than in men. Varied meanings of the F/B ratio between different investigations have been reported. We have previously reported that colitis associated colon carcinogenesis in AOM/DSS-treated mice model was induced more severely in male mice than female by way of the inflammatory mediators such as IL-1β and MPO 30 . Furthermore, we found that E2 inhibited the initiation of CRC by up-regulating nuclear factor erythroid 2-related factor 2 (Nrf2), a transcriptional factor-related pathways in the AOM/DSS-treated male ICR mice 31 . The stool from the previous experimental sets 31,35 was taken and the current study was conducted. The target protein expression such as IL-1β and MPO was strongly affected by estrogen in acute treatment condition. However, the modest changes on gut microorganisms observed in this study may be originated from too short of intervention, which is a limitation of this study. In the current study, no significant difference in the F/B ratio

In AOM/DSS group Ortholog Definition LDA effect size p-value p-value (FDR)
Abundance ratio (mean, %)  www.nature.com/scientificreports/ was observed in the Group 1 samples based on the differences in sex and the presence or absence of estrogen. Interestingly, the F/B ratio was significantly decreased upon estrogen supplementation in the male AOM/DSS group relative to Male_AOM/DSS group. These results have supported that estrogen treatment in CRC models could decrease intestinal dysbiosis by decreasing the F/B ratio. Different from gender which is related to the social and cultural role of human, sex is biologically determined based on the physiology of the reproductive systems and functions derived from the types of chromosomes or hormones. The terms of 'male' and 'female' are used to describe the sex of human participants or other sex-related factors 41 . Despite differences in sex, most preclinical studies have adopted the use of male animals, and a few have not mentioned the sex of the animal used 42 . In the present study, we used both male and female mice and investigated the impact of sex hormones, particularly, the female sex hormone 'estrogen' , on changes in the composition of the gut microbiota. Furthermore, the female mice were subjected to an OVX operation to determine the effect of endogenous estrogen. The results derived from these experiments signify the importance of our research. According to an animal study reported by Org et al., independent analysis of the gut microbiota using 89 different inbred strains of mice was conducted and clear differences in the composition of the gut microbiota and diversity between the sex within each strain was observed 25 . Particularly, obvious differences in the results between sex were observed in the C57BL/6J and C3H/HeJ strains 25 . In the total animal experimental cohort, the members of the phyla Actinobacteria and Tenericutes were more abundant in the male mice than in female mice 25 . In our present study, abundance of the phyla Tenericutes was not correlated with sex. However, its proportion was higher in male mice than in male mice supplemented with E2. Furthermore, at the species level, the abundance of Bacteroides thetaiotaomicron was found to be increased in male mice relative to those supplemented with E2, which was similar to the results observed in a previous study conducted in 2008 involving Chinese family members 43 . In another large cohort study, it was determined that females in the Netherlands had a higher proportion of Akkermansia muciniphila even after correcting all potential risk factors such as diet, lifestyle, and medication for the normal health 27 . In a Japanese study, significantly higher levels of members belonging to the genus Akkermansia were observed in the females than in males 44 and higher levels of Akkermansia muciniphila were identified in C57BL/6 background female relative to male mice 45 . Additionally, an increased abundance of Clostridium bolteae was observed in humans after bilateral ovariectomy 27 . Unlike the compositional change of the gut microbiota in humans, which is dependent on sex and the presence of estrogen, in the present study, the abundance of Akkermansia muciniphila did not increase in female mice compared to that of the male mice. Contrarily, the abundance of Akkermansia muciniphila was found to be increased in OVX female mice compared with that of the intact female mice, which was consistent with the results presented in a previous report by Choi et al. 46 Akkermansia muciniphila is a mucin-degrading bacteria that can convert mucin to short-chain fatty acids 47 . In addition to its role as a commensal bacteria, it has been reported to be involved in pro-inflammatory pathways and the activation of chemotaxis and the complement cascade 48 . Furthermore, it has been reported that the organism had exacerbated an intestinal inflammation in mice infected with Salmonella typhimurium by disturbing host mucus homeostasis 49 . Although no change in abundance of Akkermansia muciniphila based on sex was observed, its increase in the OVX female group suggested that mucin homeostasis was disrupted due to estrogen removal.

Male_Con Male_AOM/DSS Male_AOM/DSS/E2
The etiology of CRC development has been investigated in human patients and experimental animal models. It was concluded that the genetic (e.g. diabetes 50 and obesity 51 ) and dietary risk factors (e.g. high fat diet 52 and processed and red meat 53 ) are involved in its development. However, according to a recent review article, it has been suggested that the intestinal microbial community may also be an important contributing factor in the initiation and development of CRC 54 . The intestinal microflora of DSS-treated C57BL/6 background mice has been characterized by a loss in microbial diversity and changes taking place in the bacterial composition, which included a decrease in the number of bacteria from the Lactobacillus group and an increase Enterobacteriaceae, Akkermansia and Desulfovibrio 16 . Furthermore, it is well known that Fusobacterium nucleatum is highly associated with CRC development 55 . In the present study, Fusobacterium nucleatum was not identified in the experimental CRC group. Contrary to the results generated in the study involving C57BL/6 mice 16 , at the species level, a higher abundance of Lactobacillus reuteri was observed in AOM/DSS-induced CRC mice compared to that present in intact male mice. Contrarily, a higher abundance of the species Lactobacillus murinus was observed in AOM/DSS-treated male mice supplemented with E2 compared with that of the AOM/DSS-treated male mice. In the present study, Lactobacillus reuteri is a well-studied probiotic bacterium capable of colonizing in the bodies of many mammals 56 . In humans, the species has been detected from different sites of the body, including the gastrointestinal tract, urinary tract, skin, and breast milk 56 . Lactobacillus reuteri can produce antimicrobial compounds such as organic acids, ethanol, and reuterin. Due to its antimicrobial activity, Lactobacillus reuteri is able to inhibit the colonization of pathogenic microbes and remodel the composition of commensal microbiota in the host 56 . Our data suggest that Lactobacillus reuteri and Lactobacillus murinus may differently act during CRC progression.
It would have been insightful if we could simultaneously verify the effects of testosterone and estrogen, which is a limitation of this study. Testosterone, which is the primary male sex hormone, plays a pivotal role in the development of male reproductive organs such as the testes and prostate and promotes the development of secondary sexual characteristics such as muscle strength and bone mass. Additionally, testosterone is involved in the health and wellbeing of the individual 57 . Therefore, to support the limitation of this study, we are conducting experiments in which animals are subjected to orchiectomy. To confirm the impact of testosterone on the gut microbiota, which consequently is impacted by differences in sex and CRC development, a study involving orchiectomy study is needed. Unfortunately, it was difficult to find models with consistent gut microbiota among sex and affected by sex hormones. However, few were available, which could be identified based on the functional profiles of microbial communities in mice of different sex and in which, development of CRC was impacted due to the presence or absence of estrogen. Currently, we are conducting experiments to analyze gut Scientific RepoRtS | (2020) 10:12283 | https://doi.org/10.1038/s41598-020-69112-w www.nature.com/scientificreports/ microbiota according to sex differences in another strain, C57BL/6 mice. We have not verified sex differences in gut microorganisms in the AOM/DSS group, which is a limitation of this study. This is also currently under study in the C57BL/6 strain. In addition, it was regrettable that there was no consideration for progesterone therapy and that the estrogen supplementation group was not added after ovariectomy.
In conclusion, we can state that estrogen altered the gut microbiota in ICR mice, particularly in males treated with AOM/DSS, by decreasing the F/B ratio. This suggests another possibility that estrogen could cause changes in the gut microbiota, thus reducing the risk of developing CRC.
After 1 week of acclimatization, the 4-week-old female mice were subjected to ovariectomy in order to create an animal model that could mimic the physiology during menopause through the elimination of endogenous female sex hormones. The normal female control mice were subjected to a sham operation. E2 (10 mg/kg; Sigma-Aldrich, E8876) was dissolved in olive oil, and Male_E2 mice group was received a daily intraperitoneal injection of E2 for 7 days. Except Male_E2 group, others were injected olive oil as a vehicle.
To induce colitis-associated CRC, 5-week-old male mice were injected intraperitoneally with AOM (10 mg/ kg; Sigma-Aldrich, St. Louis, MO, USA) and 2.5% (w/v) DSS (MP Biomedicals, Aurora, OH, USA) was administered through their drinking water for 7 days starting 1 week following the injection of AOM 58 . E2-treated mice received an intraperitoneal injection of E2 dissolved in olive oil daily for 7 days. The injections were administered at the time of DSS consumption. Their feces were freshly collected from 59 individual male and female mice in fifteen separate cages. All fecal samples were immediately frozen in liquid nitrogen and stored at − 80 °C. The animals were euthanized by CO 2 asphyxiation (Fig. 1). The stool from the previous experimental sets 31,35 was taken and the current study was conducted. Mice were randomly divided from shipment group. All mice were housed in the same room in filter top cages, with three to five mice per cage. Mice were marked so that individual mice could be followed for the duration of the experiments. All animal experimental procedures were approved by the Institutional Animal Care and Use Committee (IACUC) of the Seoul National University Bundang Hospital (BA1310-139/091-01). The procedures were carried out in accordance with the ARRIVE (Animals in Research: Reporting of in Vivo Experiments) guideline.
DNA extraction and metagenome sequencing of the 16S rRNA gene. Genomic DNA was extracted from the frozen fecal samples using a QIAamp DNA Stool Mini Kit (Qiagen, Valencia, CA, United States) according to the manufacturer's recommendations. PCR amplification was carried out using the primers targeting the V3 to V4 regions of the 16S rRNA gene of the extracted DNA to conduct metagenome sequencing, and were as follows 59 : 341F (5′-TCG TCG GCA GCG TCA GAT GTG TAT AAG AGA CAG CCT ACGGGNGGC WGC AG-3′) and 805R (5′-GTC TCG TGG GCT CGG AGA TGT GTA TAA GAG ACA GGA CTACHVGGG TAT CTA ATC C-3′). Next procedures were performed as described previously 59 . Briefly, the identity of the PCR products was confirmed using electrophoresis. They were purified using a QIAquick PCR Purification Kit (Qiagen, Valencia, CA, USA). The purified PCR products were tagged with Illumina indices and adapters from a Nextera XT Index Kit (Illumina, San Diego, CA, USA). Short DNA fragments were eliminated using a FavorPrep Gel/PCR Purification Kit (Favorgen, Taiwan). The PCR amplicons were quantified using a Quant-iT PicoGreen dsDNA Assay Kit (Thermo Fisher Scientific, Wilmington, DE, USA). After pooling the DNA (300 ng per sample), the PCR products were purified using a FavorPrep Gel/PCR Purification Kit (Favorgen, Taiwan). The quality assessment for confirming the integrity and product size of the DNA was conducted using a Bioanalyzer 2,100 (Agilent, Palo Alto, CA, USA) using a DNA 7,500 chip at ChunLab, Inc. (Seoul, South Korea). Metagenome sequencing was carried out using the Illumina MiSeq platform at ChunLab, Inc. (Seoul, South Korea).
The raw reads were processed by first checking the quality wherein low quality (< Q25) reads were filtered out using a Trimmomatic 0.32 60 . The paired-end sequence data was merged using PANDAseq 61 . The primer sequences were then trimmed using an inhouse program of the ChunLab, Inc. (Seoul, South Korea) with a similarity cut off of 0.8. Non-specific amplicons, which did not encode 16S rRNA, were identified by using the HMMER program hmmsearch based on the 16S rRNA profiles 62 . The sequences were denoised using a DUDE-Seq 63 , and non-redundant reads were extracted through UCLUST-clustering 64 . The EzBioCloud database was utilized for performing the taxonomic assignment using USEARCH (8.1.1861_i86linux32) 64 followed by more precise pairwise alignment 65 . UCHIME 66 and the non-chimeric 16S rRNA database from EzBioCloud were used to detect the chimaeras for reads containing a best hit similarity rate of less than 97%. The sequence data was then clustered using CD-HIT 67 and UCLUST 64 .
Rarefaction for the obtained OTUs was calculated using the BIOiPLUG (ChunLab, Inc.). The reads were normalized to 7,525 to perform the analyses. The microbial diversity (Observed OTU count, Chao1, Shannon, and Simpson indices) was examined using BIOiPLUG. To visualize differences in the samples, PCoA was assessed using generalized UniFrac. The clustering of samples was explained based on the values of the principal coordinate (PC). Additionally, an Unweighted Pair Group Method with Arithmetic mean (UPGMA) tree was created Scientific RepoRtS | (2020) 10:12283 | https://doi.org/10.1038/s41598-020-69112-w www.nature.com/scientificreports/ using BIOiPLUG (ChunLab, Inc.). The beta diversity distances were calculated using generalized UniFrac. The significance of the separation between the groups was calculated by permutational multivariate analysis of variance (PERMANOVA) test. Taxonomic summary bar charts were created to determine the OTU abundance ratio (%) at the phylum and family levels using GraphPad Prism (version 5.01). The enterotype classification of the gut microbiota at the genus level was calculated using the R package "clusterSim". The optimal cluster number was determined by maximizing the value of the Calinski-Harabasz (CH) index 68 .
Biomarker discovery through the linear discriminant analysis (LDA) effect size (LEfSe). Based on the relative taxonomic abundances, the taxonomic biomarker discovery and related statistical significance were assessed by the LEfSe method 69 . The criteria for conducting the LEfSe analysis were as follows: (1) an alpha value for the factorial Kruskal-Wallis H test between assigned taxa compared to that of the groups with less than 0.05, (2) an alpha value for the pairwise Wilcoxon test among the taxonomic compositions of less than 0.05, (3) a threshold of the logarithmic LDA score for discriminative features less than 2.0, and (4) a multi-class analysis set as all-against-all. After this process, we went further to simplify the LEfSe plot. (5) overlapped bacteria selection, (6) search and classify the bacterial characteristics with reference to reported papers: commensal bacteria, opportunistic pathogens, not characterized, and (7) elimination of non-overlapped bacteria only at 'not characterized' bacteria. In LEfSe plot, all of commensal bacteria and opportunistic pathogens and only overlapped 'not characterized' bacteria are included.
To predict the functional biomarker profiles of bacterial communities, we conducted their phylogenetic investigations through reconstructing of unobserved states (PICRUSt) 70 based on the Kyoto Encyclopedia of Genes and Genomes (KEGG) database (www.kegg.jp/kegg/kegg1 .html). Discovery of the functional biomarker and associated statistical significance were assessed based on the LEfSe.

Statistical analysis.
Statistical calculations except those of the pyrosequencing data were performed using PASW Statistics version 18.0.0 (SPSS Inc., 2009, Chicago, IL, USA). The groups were compared based on the results of the Kruskal-Wallis H test. Then, the two groups were compared using the Mann-Whitney U test (also known as the Wilcoxon rank sum test), and the p-value < 0.05 was considered to be significant. For the adjustment of multiple comparisons, corrected q-values of the false discovery rate (FDR) were calculated for significance less than 5%.

Data availability
The raw unprocessed gene datasets of 16S rRNA, which were generated during the current study, are available with the NCBI Sequence Read Archive (SRA Accession Number PRJNA590320), https ://www.ncbi.nlm.nih. gov/sra/PRJNA 59032 0. The processed data generated and analyzed during this study has been included in this published article under Supplementary Dataset files.