Supplementary Information
SupplementaryMaterials and Methods
Denaturing gradient gel electrophoresis (DGGE) analysis of gut microbiota
A. PCR amplification and DGGE
DNA isolated from each faecal sample was used as template in the amplification of the V3 region of the 16S rRNA gene using the universal bacterial primers P2 (5’-ATTACCGCGGCTGCTGG-3’) and P3 (5’-CGCCCGCCGCGCGCGGCGGGCGGGGCGGGGGCACGGGGGGCCTACGGGAGGCAGCAG-3’) and the hot-start touchdown protocol described by Muyzer et al(Muyzer et al., 1993)ina thermocycler PCR system (PCR Sprint, Thermo electron, Corp., UK). Each 25μl PCR reaction mixture contained 1.5 U of rTaq DNA polymerase (Takara, Dalian, China), 2.5μl of the corresponding 10×buffer (Takara, Dalian, China), 0.2 mmol/L each deoxynucleoside triphosphate (dNTP),25 pmol of each primer, and 20 ng of total faecal DNA. After the initial amplification, a reconditioning PCR method was performed to decrease heteroduplexes formation(Thompson et al., 2002). Parallel DGGE was performed using a Dcode Systemapparatus (Bio-Rad) in an 8% (w/v) acrylamide gel with a gradient from 27% to 52% and electrophoresis in 1×Tris-acetate-EDTA (TAE) buffer at the constant voltage of 200 V and a temperature of 60oC for 240 minutes. After electrophoresis, the gels were stained with SYBR GreenⅠ(Amresco, Solon, Ohio) and visualized on a UVI gel documentation system (UVltec, Cambridge, United Kingdom).
B. Sequence analysis of DGGE bands
Important DGGE bands were excised from the original gel and incubated in 100 μl sterile distilled water at 4°Covernight. A 1-μl aliquot of elution was used for PCR amplification of the DNA fragments from the excised gelwith corresponding primers, P2 and P3. PCR products was excised from a 1.0% agarose gel and purified with a DNA Gel Extraction Kit (V-gene, Hangzhou, China). The products were ligated into the pGEM-T easy vector (Promega, Madison, WI) and transformed into competent E.coli DH5α cells (TIANGEN, China). Inserted DNA was amplified using corresponding primers and resolved byDGGE to verify the position of the original band. Then three clones migrating to the same position of the original band were sequenced (Invitrogen, Shanghai, China).
The sequences of excised DGGE bands were submitted to the RDP database (version 9.33) to determine their closest relatives with length greater than1,200 nt. The sequences obtained are available in the GenBank database under accession numbers EU584214-EU584231.
Terminal restriction fragment length polymorphism (T-RFLP)analysis of faecal samples
A.16S rRNA gene amplification and purification
DNA isolated from each faecal sample was used as template in the amplification of the 16S rRNA gene using the universal bacterial primers,8F (5’-GAGAGTTTGATCCTGGCTCAG-3’), 5’ end labelled with D4, and 1492R(5’-GGC/TTACCTTGTTACGACTT-5’)(Hayashi et al., 2002). Each 25μl PCR reaction mixture contained 1.5 U of rTaq DNA polymerase (Takara, Dalian, China), 2.5μl of the corresponding 10×buffer (Takara, Dalian, China), 0.2 mmol/L each deoxynucleoside triphosphate (dNTP),25 pmol of each primer, and 10 ng of total faecal DNA. A 20 cycles PCR program (Eckburg et al., 2005) was performed with a thermocycler PCR system (PCR Sprint, Thermo electron, Corp., UK). Each 100μl PCR product was digested by 1 U Mung bean nuclease (Promega, Madison, WI) to remove the single stranded extensions. Digestion products were purified witha DNA purification Kit (V-gene, Hangzhou, China), according to the manufacturer’s instruction.
B. T-RFLP
A 50-ng sample of 16S rRNA gene amplification product from each mouse sample was digested at 37°C for 3 h with 5 U of the following restriction endonucleases: AluⅠ, HaeⅢ, HhaⅠ, MspⅠ, or Csp6Ⅰ(Promega, Madison, WI). The efficiency of restriction digestion was confirmed by agarose gel electrophoresis, and the digested fragments were separated on a CEQTM 8000 genetic analysis system (Beckman Coulter). The sizes of the fluorescently labelled fragments were determined by comparison with the internal size standard (GenomeLabTM DNA Size Standard 600 Kit, Beckman Coulter). Fluorescence intensity data were automatically collected and subsequently analyzed by the fragment analysis software provided with the CEQTM 8000 system. Relative peak areas of each TRF were determined by dividing the area of the peak of interest by the total area of peaks within the following threshold values: a lower threshold at 60 nt and an upper threshold at 640 nt. A threshold for relative peak area was applied at 3%, and only TRFs with higher relative abundances were included in the remaining analyses. The two peaks were identified as the same if the difference of the peak size was less than 1nt. T-RFLP was repeated three times from each 16S rRNA gene amplification.
C. Statistical analysis of T-RFLP data
In the matrix used for statistical analysis, the data of three replicates from each animal were as independent objects and the relative peak height of TRFs from every restriction endonuclease was as the variable. This matrix was used for PCA. Multivariate analysis of variance (MANOVA) is a generalized form of analysis of variance (ANOVA) methods to cover cases where there is more than one (correlated) dependent variable and where the dependent variables cannot simply be combined, and the number of variables should less than the number of samples. Because the TFRs far outnumber our animals and PC1 to PC9 from PCA,which accounts for 98% of total variations, PC1 to PC9 were used for Multivariate ANOVA test. The clusters are computed by applying the single linkage method to the matrix of Mahalanobis distances between group means. H returns a vector of handles to the lines in the figure. All the above methods were implemented in Matlab® (ver. 7.1, The MathWorks, Inc.).
Pyrosequencing of 16S rRNA gene V3 region
- PCR amplification of 16S rRNA gene V3 region
For faecal samples from each mouse, the extracted DNA was used as template in the amplification of the V3 region of 16S rRNA gene. The forward primer was 5’-NNNNNNNNCCTACGGGAGGCAGCAG-3’, and the reverse primer was5’-NNNNNNNNATTACCGCGGCTGCT-3’, where the underlined sequence is the universal bacterial primer P1 and P2. The NNNNNNNN is the unique eight base barcode used to distinguish PCR product from different samples. Reaction conditions were as follows: Each 25μl PCR reaction mixture contained 0.25 U of Platinum® Pfx DNA polymerase (Invitrogen, USA), 2.5μl of the corresponding 10×Pfx amplification buffer (Invitrogen, USA), 0.5 mMof MgSO4 (Invitrogen, USA), 0.3 mmol/L each deoxynucleoside triphosphate (dNTP),and6.25 pmol of each primer, and 20 ng of total faecal DNA. PCR reactions were run in a thermocycler PCR system (PCR Sprint, Thermo electron, Corp., UK)using the following program: 3 minutes denaturing at 94°C followed by 20 cyclesof 1 minute at 94°C (denaturing), 1 minute for annealing (1°C reduced for every 2 cycles from 65°C to 57°C followed by 1 cycle at 56°C and 1 cycle at 55°C), and 1 minute at 72°C (elongation), with a final extension at 72°C for 6 minutes. Three independent PCR reactions using different barcoded primers were performed for two animals from each group.
- Gel purification and pyrosequencing
Each PCR product was isolated with the Gel/PCR DNA Fragments Extraction Kit (Geneaid, UKAS). 30 ng of each purified PCR product was mixed, and the mixture was purified from 1.2% agarose gel with the Gel/PCR DNA Fragments Extraction Kit (Geneaid, UKAS) and used as the DNA library for pyrosequencing using a GS20 platform (Roche),as described previously(Margulies et al., 2005). Based on several previous reports describing sources of errors in 454 sequencing runs(Margulies et al., 2005; Sogin et al., 2006; McKenna et al., 2008), standards wereused for quality control, as follows: if a sequence (a) shows no mismatch to the barcode and 16S rRNA gene primer at sequencing end, (b) is more than 100nt in length, (c) has no more than two undermined bases in the sequence read, and (d) finds more than 75% mach to a previously determined 16S rRNA gene sequence, then it will be regarded as usable. The unique sequences obtained in this study are available in the GenBank database under accession numbersFJ032696- FJ036862.
C.Bioinformatic and Statistical analysis
The unique V3 sequences of 16S rRNA gene from pyrosequencing were aligned using NAST multi-aligner with a minimum template length of 100 bases and a minimum percent identity of 75% (DeSantis et al., 2006). The resulting alignments were imported into the ARB (Ludwig et al., 2004) for construction of a neighbour-joining tree. The phylogenetic tree was then used for online UniFrac ( with abundance weighting (incorporating abundance data).
A distance matrix of unique sequences from ARB was imported into DOTUR for phylotype binning (Schloss and Handelsman, 2005). The abundance information of each unique sequence was added to the OTU results from DOTUR, and measured for coverage (rarefaction analysis with software Past) and diversity (Shannon index with R software (
OTUs were defined using a threshold of 97% identity, which wasa criterionfor species level delineation in previous studies(Huse et al., 2007). One sequence randomly selected from each OTU was BLAST searched against the RDP database (version 9.33) to identify the taxonomic group and inserted into pre-established phylogenetic trees of full length 16S rRNA gene sequences database from GreenGenes using ARB (hypervariable regions masked with the lanemaskPH filter) (Ludwig et al., 2004).
Errors in pyrosequencing may occur at a rate of about 0.25%, suggesting that most of the 200nt sequences will contain either 0 or 1 error. Using the Poisson model, we would expect only 1.029e-6 of the reads to contain the six errors that would be required to form a new species-level at the 97% OTU identity. Thus, it is unlikely that a single OTU in analysis was generated through that mechanism.
PLS-DA was used to test if two groups can be separated based on OTUs data. Martens’ uncertainty testwas usedto select significant OTUs which can discriminatethe different treatment groups(different diets, host genotypes, or cages). One-wayANOVA was further performed to validate these differential variables. All statements of significance are atP< 0.05. The correct prediction rate of the PLS-DA model was performed withleave-one-out CV.
Quantitative Real-time PCR of Bifidobacterium spp.
The plasmid containingthe 16S rRNA gene of Bifidobacterium spp. was prepared with 3S Spin plasmid Miniprep Kit V3.1 (Shanghai Biocolour BioScience & Technology Company) and linearized with ScalⅠlinear. The plasmid was then diluted from 2.39E+9 to 2.39E+2 (copies/μl) to make a standard curve. Real-time PCR amplification and detection were performed in the DNA Engine Opticon 2 system (MJ Research). The primers were the Bifidobacterium specific primers: Bif164-f (5’-GGGTGGTAATGCCGGATG-3’) and Bif662-r (5’-CCACCGTTACACCGGGAA-3’) (Satokari et al., 2001). The reaction mixture (25 μl) was composed of 1.5 U rTaq DNA polymerase (Takara, Dalian, China), 12.5μl of 2×SYBR greenⅠmix (Shanghai Biocolour BioScience & Technology Company),25 pmol of each primer, and 20 ng of total faecal DNA from each mouse or 1μl stander plasmid. The amplification program was as follows: one cycle of 94°C for 5 min, then 40 cycles of 94°C for 30 s, 62°C for 20 s, and 72°C for 40 s, and finally one cycle of 88°C for 10 s. The fluorescent product was detected at the last step of each cycle. Following amplification, melting temperature analysis of PCR products was performed to determine the specificity of the PCR. The melting curves were obtained by slow heating at 0.5°C/s increments from 70 to 95°C, with continuous fluorescence collection. Each sample had three replicates of quantitative analysis. The data detected by the system was analyzed using Opticon Monitor (Version 1.1).
Supplementary Results
T-RFLP analysis of gut microbiota
In the PCA based on T-RFLP data, thevariables influenced PCs 1,2, and 3are shown in Figure S8.
Pyrosequencing of 16S rRNA gene V3 region
A total of 29,343 high quality sequences of pyrosequencing were first selected based on the following criteria: if a sequence (a) shows no mismatch to the barcode and 16S rRNA gene primer at sequencing end, (b) is more than 100nt in length, (c) has no more than two undermined bases in the sequence read. All the sequences were aligned using NAST multi-aligner with a minimum template length of 100 bases and a minimum percent identity of 75% (DeSantis et al., 2006). In all, 29 sequences had no neighbours with higher than 75% homology and were discarded. A total of 29,314 useable sequence reads were then obtained. There were 4,156 unique sequences in all the samples, 3,145 of which were detected only once. Becausetwo animals from each treatment group had three replicates, we used 36 different barcodes. All barcodes were well populated, with an average of 814 sequences per sample tested and only two samples had less than 500 reads (Table S 3).
UniFrac provides a suite of tools for the comparison of microbial communities using phylogenetic information. The UniFrac analysis based on the unique sequences of all the fourtreatment groups revealed the significant impact of diet on gut microbiota (Fig.S2a), but the difference between mice having different genotypes was not obvious. These results confirmed that diets have a more significant influence on gut microbiota than genotypes. Hierarchical Clustering Analysis based on UniFrac analysis also indicated that the three replicates of the same animal were more similar to each other than to another animal (Fig.S2b). In other words, this 454 barcoded technologyshows satisfactory reproducibility.
In order to do phylotype binning and diversity estimation of our pyrosequencing data, the ARB distance matrix of sequence reads were imported to DOTUR. When sequences were condensed under 99% identity, 1,360 different operational taxonomic units (OTUs) were obtained. But the number of OTUs was reduced to 516 under 97% identity (Fig.S1a). The rarefaction analysis and Shannon Diversity Index were calculated for each of the 20 mice (Fig.S1b and 1c). For all the samples, the rarefaction curves did not reach a stable value, indicating that the actual numbers of OTUs in the samples are larger, echoing some recent reports based on pyrosequencing that the diversity of community of macaque gut or deep sea has been underestimated (McKenna et al., 2008)(Sogin et al., 2006). The curves of Shannon Diversity Index of all the samples had reached stable values. There was no significant difference of Shannon Diversity Index among the fourgroups.
Analysis of the taxa present in the mice gut communitiesindicated that bacteria of Firmicutes and Bacteroidetes were dominant, which was followed by Actinobacteria and Proteobacteria. Phylum-wide changes associated with IGT/Obesity werenot observed. More than 20 different families were found in the 20 mousefaecal samples, showing significant difference with the family level composition of gut bacteria in human (Ley et al., 2005). The relative abundance distribution of families showed individual-specific difference among the 20 mice (Fig.S5). The proportion of Erysipelotrichaceaewas the highest in most mice.
PLS-DA was used to find patterns out of this complex dataset. Using 97% identity to define OTUs, a PLS-DA scores plot with the first two components showed that animals with different diets, genotypes, or health phenotypes were separated into different classes (Fig.S3). The PLS-DA modelof Apoa-Ⅰ-/- mice fed on different diets yielded a 99% correct classification rate in leave-one-out cross validation when 2 PLS components were used. The correct classification rate of wildtype mice fed different diets, animals with different genotypes fed HFD, and animals with different genotypes fed NC are as follows: 99% using 1 PLS component, 80% using 2 PLS components, 90% using 2 PLS components. When the mice were divided into two groups based on the diet, the correct classification ratewas90% with 1 PLS component. When the mice were divided into two groups based on genotype, the correct classification ratewas75% with 2 PLS components. When the mice were divided into two groups based on health status, the correct classification ratewas95% with 2 PLS components.
Sixty-five OTUs were selected using Martens’ uncertainty test key variables for the classifications. One random sequence from each of the key OTUs was inserted into the pre-established phylogenetic tree (Fig.S4). 21 OTUs (pink marked in the tree) increased in the mice group fed onHFD compared with their counterparts fed on NC,but 26 OTUs (blue) showed oppositebehaviour. We also found 4OTUs (yellow) withgenotype-dependent reactions. For example, one phylotype, OTU 64 in Class Clostridia,washigh in theApoa-Ⅰ-/-/NC group but low in Wt/NC animals, showing a response to genotype difference (Table S1). However, it disappeared completely in HFD groups regardless of genotype, confirming the results with DNA fingerprinting that high fat diet can diminish genotype-related differences. There was 1 OTU (purple) increased in the knockout mice, which may be only associated withgenotype. There were12 OTUs (green) reduced and 1 OTU (orange) abundant in mice with IGT. Four lineages were found in Class Erysipelotrichi, M1 (9 OTUs), M2 (6 OTUs), M3 (14 OTUs), M4 (1 OTU). M1 wereabundant in healthy Wt/NC animals, but significantly reduced in theother three groups with IGT. M2 and M4 werepredominant in HFD/obese groups and M3 hadmuch higher population levels in NC/lean animals.
When one sequence randomly selected from each OTU under 97% identity was BLAST searched against the RDP database (version 9.33) to identify the taxonomic group and inserted into pre-established phylogenetic trees of full length 16S rRNA gene sequences using ARB, more OTUs were found inthe fourlineages in Class Erysipelotrichi. The total OTUs in the four lineages M1-4 were36, 30, 44, and 4, respectively. Those OTUs which were not selected as key variables for separating classes weremostly rare phylotypes. However, most rare phylotypes showed similar behaviour to the identified ones in the same lineage (Fig.S6).