Compositional distinction of gut microbiota between Han Chinese and Tibetan populations with liver cirrhosis

Ethics statement

All individuals were recruited from the Department of Digestive Medicine, Hospital of Chengdu Office of People’s Government of the Tibetan Autonomous Region. The study was approved by the Medical Ethics Committee of the Hospital of Chengdu Office of People’s Government of Tibetan Autonomous Region in Chengdu, China (approval number: (2017) study No. 10). Written informed consent was obtained from each subject in this study. Inclusion criteria: patients with an identified risk factor for cirrhosis; portal hypertension and liver dysfunction; evidence of LC on abdominal CT or color doppler ultrasound imaging (Acharya & Bajaj, 2017; Acharya & Bajaj, 2019).

Human subjects

A total of 55 individuals were divided into four groups: Tibetan healthy control (THC) group (n = 10), Han Chinese healthy control (HHC) group (n = 9), Tibetan patients with LC group (n = 20), and Han Chinese patients with LC group (n = 16).

Fecal sample collection and genomic DNA extraction

Fecal samples of all the subjects in this study were collected and immediately frozen at −80 °C for further analyses. Total genomic DNA was extracted from the frozen fecal samples using the DNeasy PowerSoil Kit (Tiangen Biotechnology Company, Beijing, China) following the manufacturer’s instructions. Quality and quantity of DNA were verified with the NanoDrop spectrophotometer and agarose gel. Extracted DNA was diluted to a concentration of one ng/μl and stored at −20 °C until further processing.

16S rRNA gene analysis of fecal samples

The 16S rRNA gene analysis of the fecal samples was performed as previously described (Qiu et al., 2019). Briefly, the diluted DNA obtained above was used as a template for polymerase chain reaction (PCR) amplification of the bacterial 16S rRNA gene with V3-V4 region-specific primers: 343F (5′-TACGGRAGGCAGCAG-3′) and 798R (5′-AGGGTATCTAATCCT-3′). An Illumina sequencing adaptor was added to the primers and the reverse primer contained a sample barcode. The amplicon quality was visualized using gel electrophoresis before the PCR products were purified. Next, concentrations were adjusted for sequencing on an Illumina MiSeq PE300 system (OE Biotech Co., Ltd, Shanghai, China).

Data analyses

Bioinformatics analyses were performed as previously described (Derakhshani et al., 2016). Raw sequencing data was generated in the FASTQ format. Paired-end reads were preprocessed using Trimmomatic (v0.35) (Bolger, Lohse & Usadel, 2014) to detect and cut off ambiguous bases Low quality sequences were cut off using the sliding window trimming approach and the paired-end reads were assembled using FLASH (v1.2.11) (Reyon et al., 2012). Then using QIIME (v1.8.0) (Caporaso et al., 2010) and UCHIME (v4.2) (Edgar et al., 2011), the sequences were denoised as follows: reads with ambiguous sequences, homologous sequences, or below 200bp were abandoned and reads with 75% of bases above Q20 were retained. All sequences with mismatches or ambiguous calls in the overlapping region were discarded. Clean reads were subjected to primer sequences removal and clustering to generate operational taxonomic units (OTUs) using VSEARCH (v2.4.2) (Rognes et al., 2016) software with a 97% similarity cutoff (Edgar, 2013). All representative reads were annotated and blasted against the SILVA database (v123) (or Greengenes) (16srDNA) using the RDP classifier (v2.2) (confidence threshold was 70%) (Wang et al., 2007). The phylotypes were computed as percent proportions based on the total number of sequences in each sample. Alpha- and beta-diversity indexes were calculated using QIIME as previously described (Bindels et al., 2016). To visualize α diversity, we employed QIIME software to calculate α estimators for each sample including the Shannon index, Chao1 index, Simpson’s Diversity Index, and Observed Species. In addition, β diversity between the communities was assessed by weighted UniFrac for Principal Coordinate Analysis (PCoA).

Statistical analysis

Statistical analyses were performed using SPSS software version 20.0 (SPSS Inc., Armonk, NY, USA) and GraphPad Prism (V.6.0). Principal coordinate analysis (PCoA) was applied on the resulting distance matrices to generate plots using the default settings of the PRIMER-6 software (Derakhshani et al., 2016). Comparison of OTUs and taxonomy abundances was calculated using the Kruskall–Wallis test or Mann Whitney analysis. Following statistical analyses with multiple comparisons, P values were corrected using the Benjamini–Hochberg method to control the false discovery rate (FDR) (Freedberg et al., 2015). The resultant P-values were FDR corrected with a significance threshold of 5%. Furthermore, the linear discriminant analysis (LDA) and the linear discriminant analysis effect size (LEfSe) measurements were used to find unique bacterial taxa among different groups (Segata et al., 2011) (LDA > 2, FDR-p < 0.05) and P-values < 0.05 were considered statistically significant.

Community richness and biodiversity

In this study, the related characteristics of patients and healthy controls (Table 1) indicated that there was no difference in age, sex, and BMI (body mass index). A total of 55 fecal samples were collected from the four groups and underwent PCR amplification of the V3-V4 hypervariable regions of the bacterial 16S rRNA gene. Raw data for the 16S rRNA gene has been deposited in the NCBI SRA database with the accession number PRJNA603295. After filtration by QIIME quality filters an average of 50,770, 19,914, 48,338 and 19,963 clean tags were obtained from the THC, HHC, TLC and HLC groups, respectively and among these clean tags, an average of 39,053 (76.92%), 15,143 (76.04%), 36,872 (76.28%) and 15,006 (75.17%) tags were valid in their respective groups (Table S1). The lengths of all the valid tags were 428 ~ 431 bp (Table S1). All valid tags were clustered to operational taxonomic units (OTUs) using VSEARCH software at a 97% similarity threshold and a total of 1,579 unique OTUs were generated. Rarefaction curves of the healthy and liver cirrhosis humans displaying the richness and evenness of the gut microbiota in the samples can be found in Fig. 1. As the sequences increased, the amounts of Observed Species gradually stabilized. Additionally, the Good’s coverage indexes were over 99% in all samples (Table S2). These results collectively indicated that sequencing depth in our study was sufficient, and the sequencing results were reasonable and reliable enough for further analyses.

Compositional distinction of gut microbiota of all the groups between Han Chinese and Tibetan populations

After obtaining valid tags, the α diversity and β diversity metrics showed differences among the four groups. In general, the richness of α diversity, as determined by the Chao1 index or Observed Species, was significantly higher in Han Chinese populations compared to Tibetan populations (P < 0.01), while no significance was found in either the Shannon index or Simpson’s Diversity Index representing diversity and community evenness (P > 0.05). In the HHC group, the Chao1 index and Observed Species were higher compared to the THC group (P = 0.0015 and P = 0.0019 respectively); this trend was similar in patients with LC. In the HLC group, the Chao1 index or Observed Species were significantly higher than that of the TLC group (P = 0.0006 and P = 0.0001 respectively). As for β diversity analysis, in weighted UniFrac PCoA the first principal coordinate (PC1) occupied 41.92% of inter sample variance and revealed that inter-individual variation is greater than intra-individual variation, indicating a clear distinction between the gut microbiota of the THC and HHC groups (P = 0.011) (Fig. 2B). These findings suggested that ethnic factors have an important influence on intestinal flora.

The V3–V4 sequences from the THC and the HHC groups were predominately assigned to eleven phyla: Bacteroidetes, Firmicutes, Proteobacteria, Actinobacteria, Fusobacteria, Tenericutes, Verrucomicrobia, Lentisphaerae, Spirochaetae, Cyanobacteria and Acidobacteria; the relative abundance bar plot of the gut microbiota community at the phylum level was displayed in Fig. 2D. The results showed that the dominant phyla were Firmicutes and Bacteroidetes in both the THC and the HHC groups. Collectively, these results suggested that there was an ethnicity-based diversity of gut microbiota between Tibetan and Han Chinese populations as previously reported (Lai-Yu et al., 2014).

Identification of gut microbiota significantly varied in Tibetan populations with liver cirrhosis patients

LEfSe analysis of taxon abundance was conducted to detect the differences in fecal bacteria between Tibetan patients with LC and Han Chinese patients with LC. Thus, the microbial signature in both two groups was identified. The relative abundances of Clostridiaceae_1, Clostridium_sensu_stricto_1, Barnesiella, Peptococcus, Ambiguous_taxa, Gastranaerophilales, Peptococcaceae, Melainabacteria, Succinivibrio, Succinivibrionaceae, and Aeromonadales as well as some uncultured bacterium were more prevalent in Tibetan healthy populations than t in Tibetan patients with liver cirrhosis (Fig. 3). However, more Prevotella_2, Streptococcus, Streptococcaceae, Lactobacillales, and Bacilli were found in Tibetan patients with liver cirrhosis (Fig. 3). Among these bacteria, Streptococcaceae has been reported to be prevalent in patients with cirrhosis (Yanfei et al., 2011). In addition, the LEfSe analysis between cirrhosis patients from both ethnic groups is in Fig. S1.

Comparison of gut microbiota between Han Chinese and Tibetan populations with liver cirrhosis

We investigated the microbial communities of Han Chinese and Tibetan patients with liver cirrhosis to explore whether different microbial composition exists between these two groups. The information of all the diversity indexes for Tibetan and Han Chinese patients with liver cirrhosis is in Table S2. As a result, the PCoA plots showed a clear separation between the HLC and TLC groups on the graph (Fig. 4A). The percentages of variation explained by PC1 and PC2 were 56.45% and 10.67%. Two distinct areas from samples of two groups on the graph indicated a significant difference between the gut bacterial compositions. All taxonomic assignments and their abundances are in Table S3.

Then, we analyzed gut microbiota composition at the phylum level in both groups and the relative abundance chart is in Fig. 4B. The top 15 most abundant bacteria in phylum level were Bacteroidetes, Firmicutes, Proteobacteria, Actinobacteria, Fusobacteria, Tenericutes, Acidobacteria, Gemmatimonadetes, Nitrospirae, Lentisphaerae, Verrucomicrobia, Spirochaetae, Chlorobi, Cyanobacteria and TA06, seven of which were in accordance with previous studies (Yanfei et al., 2011). Additionally, the relative abundance of Bacteroidetes and Firmicutes at the phylum level were significantly different between the TLC and HLC groups respectively (P < 0.001 and P = 0.008, Fig. 4C).

Then, a percentage stacked histogram was used to display the relative abundance of gut microbiota at the genus level; the top 30 most abundant bacteria at the genus level are shown in Fig. 5A. In both the TLC and HLC groups, Prevotella_9 was the predominant genus. In total, eight significantly different pairs were found in Alloprevotella (P < 0.001), Anaerostipes (P < 0.001), Bifidobacterium (P = 0.03), Blautia (P = 0.0035), Dorea (P =<0.001), Prevotella_2 (P = 0.0015), Prevotella_7 (P < 0.001) and Prevotella_9 (P = 0.0028). In the TLC group, the relative abundance of Anaerostipes, Bifidobacterium, and Blautia was higher than in the HLC group (Fig. 5B). However, the relative abundance of Alloprevotella, Dorea, Prevotella_2, Prevotella_7, and Prevotella_9 were decreased in the TLC group. Taken together, these results indicated a remarkable difference between Tibetan and the Han Chinese patients with LC in terms of constituent gut microbiota, suggesting the potential influence of ethnic and environmental factors in the development of LC.