Pyrosequencing analysis of the microbial diversity of airag, khoormog and tarag, traditional fermented dairy products of mongolia.

Here, we used pyrosequencing to obtain a detailed analysis of the microbial diversities of traditional fermented dairy products of Mongolia. From 22 Airag (fermented mare's milk), 5 Khoormog (fermented camel's milk) and 26 Tarag (fermented milk of cows, goats and yaks) samples collected in the Mongolian provinces of Arhangai, Bulgan, Dundgobi, Tov, Uburhangai and Umnugobi, we obtained a total of 81 operational taxonomic units, which were assigned to 15 families, 21 genera and 41 species in 3 phyla. The genus Lactobacillus is a core bacterial component of Mongolian fermented milks, and Lactobacillus helveticus, Lactobacillus kefiranofaciens and Lactobacillus delbrueckii were the predominant species of lactic acid bacteria (LAB) in the Airag, Khoormog and Tarag samples, respectively. By using this pyrosequencing approach, we successfully detected most LAB species that have been isolated as well as seven LAB species that have not been found in our previous culture-based study. A subsequent analysis of the principal components of the samples revealed that L. delbrueckii, L. helveticus, L. kefiranofaciens and Streptococcus thermophilus were the main factors influencing the microbial diversity of these Mongolian traditional fermented dairy products and that this diversity correlated with the animal species from which the milk was sourced.

INTRODUCTION

Traditional fermented dairy products play an important role in the Mongolian diet because of their nutrient richness and medicinal potential [1, 2]. The nomads of Mongolia produce various kinds of traditional fermented dairy products. Airag is a mildly alcoholic, sour-tasting fermented drink that is usually made from the raw milk of mares; it is called Koumiss in Kazakhstan, Kyrgyzstan and Russia and Chigee in Inner Mongolia, China [3, 4]. Khoormog is a traditional fermented mild alcoholic beverage made from raw camel milk [5]. Tarag is a yogurt-type traditional fermented milk made from heat-treated cow, yak or goat milk [6]. These products are prepared by time-honored methods that were developed by the nomadic people in their gers (portable houses) and are naturally fermented by adding the milk to a traditional container without the use of commercial starters. Therefore, these products have unique microbial compositions depending on the individual houses in which they are prepared.

There have been numerous analyses of the lactic acid bacteria (LAB) and yeasts in the traditional fermented dairy products of Mongolia that have used culture-based methods [7,8,9]. However, in these studies, the species were identified on the basis of their phenotypic features, and the information obtained has not always contributed to an accurate and detailed picture of the microbial diversity of these products. Recently, some studies reported the detailed microbial composition of Mongolian dairy products by using a combination of culture-based method and molecular biological identification [3, 10]. In our previous study, we used culture-based isolation and molecular-based identification—random amplified polymorphic DNA (RAPD) PCR for strain typing and 16S rRNA gene sequencing for species identification—and reported on the diversity of LAB and yeasts in 22 Airag and 31 Tarag samples collected from various regions of Mongolia. A total of 367 LAB strains isolated from these samples revealed significant differences in LAB diversity. Specifically, we found that the predominant LAB species of Airag were Lactobacillus helveticus and Lactobacillus kefiranofaciens, and those of Tarag were Lactobacillus delbrueckii subsp. bulgaricus, Lactobacillus fermentum, L. helveticus, L. kefiranofaciens and Streptococcus thermophilus. However, in that study, differences in the LAB composition of the Tarag samples in relation to the animal species from which the milk was sourced were not observed [3].

Takeda et al. [10] reported a result similar to our previous result but with the slight difference that L. delbrueckii subsp. lactis and L. fermentum were the predominant species in Airag rather than L. kefiranofaciens.

To evaluate the features of Mongolian dairy products and their utility as probiotics, it is very important to conduct detailed analyses of their microbial diversity. To avoid the inherent disadvantages of culture methods, such as underestimation of the viable cell counts, culture-independent methods, such as clone libraries or denaturing gradient gel electrophoresis (DGGE), are often used. Although there have been a few studies [11, 12] on the bacterial diversity of Airag and Tarag using the DGGE approach, too few culture-independent approaches have been performed to clarify the microbial diversity of these targets. The 16S rRNA gene sequence-based pyrosequencing method enables a detailed, comprehensive and high-throughput analysis of microbial ecology, and this method has been applied to various traditional fermented food research studies [13,14,15]. In particular, the V1−V2 hypervariable region of the 16S rRNA gene has high frequencies of sequence variability and an outstanding ability to reproduce the full-length 16S rRNA gene-based taxonomic classification [16, 17].

In this study, we used the pyrosequencing method based on the sequence of the V1−V2 hypervariable region of the 16S rRNA gene to obtain a detailed analysis of the bacterial diversity of the 22 Airag, 5 Khoormog and 26 Tarag samples that we used in our previous study [3].

MATERIALS AND METHODS

Sample collection

All samples were collected from the Mongolian provinces of Arhangai, Bulgan, Dundgobi, Tov, Uburhangai and Umnugobi in July 2004 [3]. About 3 mL of samples were collected and stored at 4°C in a vehicle-mounted refrigerator. Subsequently, all of the samples were transported by air to the Yakult Central Institute for Microbiological Research, Tokyo, Japan, at below freezing and stored at –20°C until they were used for DNA extraction. We used 22 alcoholic beverages (Airag), 5 Khoormog samples and 26 Tarag samples made from milks of domestic animals. Although camel fermented milk samples were categorized as Tarag in our previous study, we recategorized them as Khoormog in this study in accordance with the proper Mongolian name [5], because raw camel milk was used for these samples (Table 1).

Table 1.

Properties, pyrosequencing data and alpha diversity scores of the samples used

No.	Fermented milk		Sampling location			Unfiltered sequence	Filter passed sequence	Chimeric sequence No. (Percentage)	Observed OTUs	Chao 1	Good’s coverage	PD whole tree	H’
	Type	Species origin	Province	City	Region
A01	Airag	Mare	Arhangai	Khashaat	Forest Steppe	3346	2635	1 (0.04)	15	20	99.8	0.85	0.56
A02	Airag	Mare	Arhangai	Tsagaannuur	Forest Steppe	3224	2460		14	34	99.7	0.77	0.56
A03	Airag	Mare	Arhangai	Tsagaannuur	Forest Steppe	3499	2756		11	12	99.7	0.66	0.28
A04	Airag	Mare	Arhangai	Tsagaannuur	Forest Steppe	3382	2718		10	10	100	0.61	0.39
A05	Airag	Mare	Arhangai	Tsagaannuur	Forest Steppe	3639	2909	2 (0.07)	17	18	99.6	0.90	0.41
A06	Airag	Mare	Arhangai	Tsagaannuur	Forest Steppe	3381	2678		18	19	99.9	0.89	0.98
A07	Airag	Mare	Arhangai	Tsenkher	Forest Steppe	3415	2729		11	16	99.9	0.68	0.20
A08	Airag	Mare	Arhangai	Tsenkher	Forest Steppe	3405	2680		11	16	99.9	0.62	0.40
A09	Airag	Mare	Arhangai	Tsenkher	Forest Steppe	3206	2575		21	23	99.9	1.14	1.32
A10	Airag	Mare	Arhangai	Tsenkher	Forest Steppe	3570	2828	1 (0.04)	18	18	99.8	1.01	1.14
A11	Airag	Mare	Arhangai	Tsenkher	Forest Steppe	3578	2883	1 (0.03)	16	17	99.9	0.78	0.71
A12	Airag	Mare	Arhangai	Tsenkher	Forest Steppe	3731	2905		19	20	99.9	1.00	1.60
A13	Airag	Mare	Bulgan	Khishig-Ondor	Steppe	3630	2952	1 (0.03)	19	30	99.8	1.04	0.45
A14	Airag	Mare	Tuv	Bayan-Onjuul	Steppe	3109	2471	1 (0.04)	16	17	99.8	0.86	0.68
A15	Airag	Mare	Tuv	Erdenesant	Steppe	3703	2992		16	19	99.9	0.86	0.27
A16	Airag	Mare	Uburhangai	Bat-Olzii	Forest Steppe	3584	2883		10	10	99.9	0.46	0.54
A17	Airag	Mare	Uburhangai	Bat-Olzii	Forest Steppe	3040	2347	1 (0.04)	26	34	99.7	1.35	2.25
A18	Airag	Mare	Umnugobi	Dalanzadgad	Gobi Desert	3490	2612	4 (0.15)	30	40	99.9	1.19	1.79
A19	Airag	Mare	Umnugobi	Dalanzadgad	Gobi Desert	3695	2796	2 (0.07)	29	31	99.6	1.33	1.76
A20	Airag	Mare	Umnugobi	Dalanzadgad	Gobi Desert	3535	2840		17	18	99.9	0.81	0.74
A21	Airag	Mare	Umnugobi	Dalanzadgad	Gobi Desert	3337	2613		24	45	99.9	0.91	1.11
A22	Airag	Mare	Umnugobi	Hanhongor	Gobi Desert	3406	2713		14	16	99.9	0.72	0.62
Subtotal for Airag						75905	59975	14 (0.02)	60 (17.4 ± 5.8)
K01	Khoormog	Camel	Dundgobi	Mandalgovi	Gobi Desert	3250	2573	2 (0.08)	16	18	100	0.95	0.81
K02	Khoormog	Camel	Dundgobi	Mandalgovi	Gobi Desert	3115	2459		19	20	99.7	0.78	1.18
K03	Khoormog	Camel	Umnugobi	Hanhongor	Gobi Desert	3294	2545		22	24	99.8	1.09	2.04
K04	Khoormog	Camel	Umnugobi	Hanhongor	Gobi Desert	3190	2418	15 (0.62)	20	26	99.6	0.91	1.49
K05	Khoormog	Camel	Umnugobi	Hanhongor	Gobi Desert	3176	2522		17	23	99.7	0.76	1.43
Subtotal for Khoormog						16025	12517	17 (0.13)	39 (18.8 ± 2.4)
T01	Tarag	Cow	Arhangai	Khashaat	Forest Steppe	4004	3087		7	7	99.9	0.68	0.50
T02	Tarag	Cow	Bulgan	Dasinchilen	Steppe	3830	3050		9	9	99.9	0.63	0.82
T03	Tarag	Cow	Tuv	Erdenesant	Steppe	4086	3267		9	9	99.9	0.64	0.78
T04	Tarag	Cow	Tuv	Erdenesant	Steppe	3447	2709		10	10	99.9	0.62	1.80
T05	Tarag	Cow	Tuv	Erdenesant	Steppe	4096	3203	1 (0.03)	14	14	99.9	0.78	1.07
T06	Tarag	Cow	Tuv	Erdenesant	Steppe	3562	2833	1 (0.04)	10	11	99.9	0.69	1.10
T07	Tarag	Cow	Tuv	Erdenesant	Steppe	3584	2861		10	11	99.8	0.70	1.35
T08	Tarag	Goat	Umnugobi	Dalanzadgad	Gobi Desert	3582	2865		14	21	99.8	0.76	0.56
T09	Tarag	Goat	Umnugobi	Dalanzadgad	Gobi Desert	4316	3420		9	10	100	0.56	0.48
T10	Tarag	Goat	Umnugobi	Dalanzadgad	Gobi Desert	4256	3386		6	6	100	0.40	0.68
T11	Tarag	Goat	Umnugobi	Dalanzadgad	Gobi Desert	4204	3305	1 (0.03)	9	9	99.9	0.59	0.83
T12	Tarag	Goat	Umnugobi	Dalanzadgad	Gobi Desert	3494	2754	2 (0.07)	25	27	100	1.14	1.35
T13	Tarag	Goat	Umnugobi	Dalanzadgad	Gobi Desert	3450	2716		13	15	99.7	0.77	0.77
T14	Tarag	Yak	Arhangai	Khashaat	Forest Steppe	4054	3214	2 (0.06)	11	13	99.9	0.68	1.17
T15	Tarag	Yak	Arhangai	Tsagaannuur	Forest Steppe	4007	3247	1 (0.03)	6	6	100	0.37	0.75
T16	Tarag	Yak	Arhangai	Tsagaannuur	Forest Steppe	3820	3021		8	9	99.9	0.49	0.89
T17	Tarag	Yak	Arhangai	Tsagaannuur	Forest Steppe	4043	3230	2 (0.06)	7	7	99.9	0.44	0.61
T18	Tarag	Yak	Arhangai	Tsagaannuur	Forest Steppe	3770	3021		7	7	100	0.49	1.05
T19	Tarag	Yak	Arhangai	Tsagaannuur	Forest Steppe	4120	3288	1 (0.03)	11	12	99.9	0.57	0.77
T20	Tarag	Yak	Arhangai	Tsenkher	Forest Steppe	4326	3493	2 (0.06)	7	7	100	0.48	0.91
T21	Tarag	Yak	Arhangai	Tsenkher	Forest Steppe	4326	3471	1 (0.03)	8	9	100	0.51	0.90
T22	Tarag	Yak	Arhangai	Tsenkher	Forest Steppe	4251	3403	1 (0.03)	7	7	100	0.48	0.60
T23	Tarag	Yak	Arhangai	Tsenkher	Forest Steppe	4164	3257	1 (0.03)	10	14	100	0.58	0.90
T24	Tarag	Yak	Arhangai	Tsenkher	Forest Steppe	3105	2439		10	10	100	0.76	1.00
T25	Tarag	Yak	Uburhangai	Bat-Olzii	Forest Steppe	3772	2981		10	10	100	0.61	0.79
T26	Tarag	Yak	Uburhangai	Bat-Olzii	Forest Steppe	3289	2572		12	15	99.9	0.67	2.00
Subtotal for Tarag						100958	80093	16 (0.01)	45 (10.0 ± 3.8)
Total						192888	152585	47 (0.03)	81 (13.9 ± 5.6)

DNA extraction

DNA was extracted as described previously [3]. Collected samples were centrifuged at 20,000 × g for 3 min. The pellet was suspended in 250 µL of extraction buffer (100 mM Tris–HCl, 40 mM EDTA, pH 9.0) and 500 µL of benzyl chloride; 0.7 g of glass beads (0.1 mm in diameter) were added to the suspension, and the mixture was shaken vigorously for 30 sec with a FastPrep FP120 (Qbiogene, Carlsbad, CA, USA) at a speed of 6.5 m/sec. Subsequently, 50 µL of 10% SDS was added to the suspension, which was then vortexed vigorously at 50°C for 20 min in a MicroIncubator M-36 (Taitec, Tokyo, Japan). The mixture was cooled on ice for 15 min after the addition of 150 µL of 3 M sodium acetate. After centrifugation of the mixture at 20,000 × g for 15 min, the supernatant was collected, and DNA was obtained by isopropanol precipitation. Finally, the DNA was diluted to 10 µg/ml with TE buffer (10 mM Tris–HCl, 1 mM EDTA, pH 8.0) and stored at –80°C until use.

Pyrosequencing based on the V1–V2 hypervariable region of the 16S rRNA gene

The PCR conditions were designed as previously described [14] with slight modifications. To amplify the V1–V2 hypervariable region, we used the universal primers 27F-mod (5′-AGRGTTTGATYMTGGCTCAG-3′) and 338R (5′-TGCTGCCTCCCGTAGGAGT-3′) [18]. For the first PCR step, 25 µL of the reaction reagent contained 10 ng of extracted DNA, 10 mM Tris-HCl (pH 8.3), 50 mM KCl, 2 mM MgCl2, 200 μM of each dNTP, 5 pmol of each primer, and 0.625 U Ex Taq® HS (Takara Bio, Shiga, Japan). The PCR conditions were as follows: 98°C for 2.5 min; 15 cycles at 98°C for 15 sec, 50°C for 30 sec and 72°C for 20 sec; and finally 72°C for 5 min. Subsequently, 27F-mod with 53 kinds of barcode-sequence tag, which comprised 10 bp nucleotides, provided by Roche Diagnostics (Basel, Switzerland), and 338R were used for the second PCR. Each primer also had an additional adapter sequence on its 5ʹ end, which was required for the subsequent pyrosequencing reactions. For the second PCR step, 50 µL of the reaction reagent contained 1 µL of the first step product, 10 mM Tris-HCl (pH 8.3), 50 mM KCl, 2 mM MgCl2, 200 μM of each dNTP, 5 pmol of each primer and 1.25 U Ex Taq® HS. The PCR conditions were as follows: 98°C for 2.5 min; 12 cycles at 98°C for 15 sec, 54°C for 30 sec and 72°C for 20 sec; and finally 72°C for 5 min. The PCR products were purified by using an AMPure® XP Kit (Beckman Coulter, Brea, CA, USA) according to the manufacturer’s protocol. The purified products were quantified with a Quant-iTTM PicoGreen® dsDNA Assay Kit (Invitrogen, Carlsbad, CA, USA). All samples were adjusted to 109 molecules dsDNA/μL with TE buffer, and equal volumes were mixed. The mixed samples were then purified and quantified again with a MinElute PCR Purification Kit (Qiagen, Hilden, Germany) and Quant-iTTM PicoGreen® dsDNA Assay Kit. The amplicon mixture was applied to a Genome Sequencer GS Junior Bench Top System (Roche Diagnostics).

Pyrosequencing data processing

Sequence data processing was performed by using the QIIME software package v1.6.0 [19]. Briefly, sequences were filtered according to the following limiting quality-check parameters: a minimum quality score of 25, a maximum forward primer mismatch of 0, a maximum reverse primer mismatch of 2, a minimum read length of 300 bp, a maximum of 6 homopolymers and a maximum number of ambiguous bases of 0. A batch of sequence reads was then sorted to each sample according to the barcode-sequence tag. The USEARCH algorithm [20] was used for sequence clustering to operational taxonomic units (OTUs) at the 97% similarity level and for chimera checks. Chimera analyses were performed by using de novo and reference-based chimera checks against the pre-built 16S rRNA reference database provided in the Microbiome Utilities (http://microbiomeutil.sourceforge.net). Then sequences that were determined to be chimeric were removed. The filtered and chimera-removed sequence reads were then used to calculate alpha diversity scores, Chao 1 indices [21], the phylogenetic diversity (PD) whole tree [22], Shannon indices (H′) [23] and Good’s coverage [24], which were in turn used for the subsequent analyses of bacterial diversities.

Analyses of bacterial diversities

Reference sequences for each OTU were classified at the phylum, family, genus and species level by using DNASIS® Taxon (Hitachi Solutions, Tokyo, Japan) against the bacterial 16S rRNA gene sequences of the Ribosomal Database Project database (Release 10, Update 31). At the species level, OTUs were assigned to the species names that showed the highest similarity score with a threshold of 97%, and OTUs that did not show over 97% similarity against any known species were described as unknown species (e.g., Lactobacillus sp. A, Lactobacillus sp. B, Lactococcus sp.). Any OTUs whose relative abundances were below 0.1% in each sample were placed in the “low abundance” cluster. For OTUs that were not assigned to a known species, phylogenetic trees were drawn with their closest related known species. Multiple alignment and construction of phylogenetic trees were performed with ClustalX ver. 2.0.12. [25]. Approximately 300 bp of the 16S rRNA gene was used to construct the phylogenetic trees by using the neighbor-joining method [26]. The statistical reliability of trees was evaluated by bootstrap analysis of 1000 replicates [27]. Bacterial diversities were analyzed based on the relative abundance of OTUs at the species level by using principal component analysis (PCA) in the BioNumerics software (Applied Maths BVBA, Sint-Martens-Latem, Belgium).

Statistical analyses

Tukey’s test and the χ2 test were used to compare average scores and detection rates, respectively.

RESULTS

Pyrosequencing data

The V1−V2 hypervariable region of the 16S rRNA gene was amplified from 22 Airag, 5 Khoormog and 26 Tarag samples, and a total of 192,888 sequence reads were obtained. The read number for each sample varied from 3040 to 4326, with the average being 3639. A total of 152,585 sequence reads passed the quality check and barcode-sequence tag sorting; 47 sequence reads (0.03% of filtered sequence reads) were detected as chimeras and were therefore removed. The sequence reads clustered into 81 OTUs, with an average of 13.9 ± 5.6 per sample. For most samples, the numbers of observed OTUs were close to the Chao 1 indices. In addition, Good’s coverage at the 97% similarity level for the samples, which provides an estimate of sampling completeness, was above 99.6%, with the average being 99.9% ± 0.1% (Table 1). These results indicate that the majority of bacterial phylotypes could be identified.

Alpha diversity scores

The average alpha diversity scores, observed OTUs, Chao 1 indices and PD whole trees were significantly higher (p<0.05) in the Airag and Khoormog samples than in the Tarag samples, whereas no significant differences were observed between the Airag and Khoormog samples. Moreover, the average Shannon index score was not significantly different between the Airag, Khoormog and Tarag samples (Table 2). Among the Tarag samples, we observed no differences between the alpha diversity scores of any of the combinations of animal species of milk.

Table 2.

Comparison of average alpha diversity scores

	Observed OTUs	Chao 1	PD whole tree	H’
Airag (n = 22)	17.4 ± 5.8a	21.9 ± 9.5a	0.9 ± 0.2a	0.9 ± 0.6
Khoormog (n = 5)	18.8 ± 2.4a	22.1 ± 3.2a	0.9 ± 0.1a	1.4 ± 0.5
Tarag (n = 26)	10.0 ± 3.8b	11.0 ± 4.8b	0.6 ± 0.2b	0.9 ± 0.4

a,b Different letters indicate statistically significant differences at p<0.05.

Bacterial diversity

After we had filtered out the OTUs with low relative abundance, we detected a total of 15 families, 21 genera and 41 species in 3 phyla (Actinobacteria, Firmicutes and Proteobacteria) in the Airag, Khoormog and Tarag samples. At the family level, the following 15 families were observed: Acetobacteraceae, Aeromonadaceae, Bacillaceae, Bifidobacteriaceae, Caulobacteraceae, Clostridiaceae, Enterobacteriaceae, Enterococcaceae, Lactobacillaceae, Leuconostocaceae, Micrococcaceae, Moraxellaceae, Planococcaceae, Staphylococcaceae and Streptococcaceae. Among these families, Lactobacillaceae (which comprises only genus Lactobacillus) predominated in all samples at a rate of 68.7%−99.6%. Streptococcaceae, which comprises 2 genera, Lactococcus and Streptococcus, was the second most common family (with 10%−29.6% in 7 samples and 1.1%−9.5% in 27 samples, respectively). The relative abundance (p<0.05) and detection rate (p<0.01) of Streptococcus in Tarag were significantly higher than those in Airag samples, whereas the relative abundance and detection rate of Lactococcus in Airag were significantly higher (p<0.05) than those in Tarag. Acetobacteraceae, which comprises the genus Acetobacter, was detected in 13 Airag, 5 Khoormog and 5 Tarag samples; the detection rates of Acetobacter in Airag (p=0.05) and Khoormog (p<0.01) were markedly higher than that in Tarag. Caulobacteraceae (which comprises genus Brevundimonas) was detected in almost all of the samples (22 Airag, 5 Khoormog and 22 Tarag) at a low abundance (0.1%–2.7%). Leuconostocaceae (which comprises genus Leuconostoc) was detected in only 8 Airag and 2 Khoormog samples. Micrococcaceae, which comprises the 3 genera Arthrobacter, Citricoccus and Kocuria, was detected in 17 Airag, 3 Khoormog and 10 Tarag samples; the detection rate of Arthrobacter in Airag (16 samples) was significantly higher (p<0.05) than that in Tarag. Staphylococcaceae (which comprises genus Macrococcus) was detected in only Airag (7 samples) and Khoormog (2 samples). Aeromonadaceae, Bacillaceae, Bifidobacteriaceae, Clostridiaceae, Enterobacteriaceae, Enterococcaceae, Moraxellaceae and Planococcaceae were detected in a few samples as minor microbial components (Table 3). Thus, at the family to genus level, significant diversity was not observed between the Airag and Khoormog samples.

Table 3.

Bacterial diversity of Airag, Khoormog and Tarag at the phylum to genus level

At the species level, a total of 41 OTUs were detected in all of the samples, and 5 genera, Acetobacter, Enterococcus, Lactobacillus, Lactococcus and Streptococcus, were classified into 4, 2, 13, 3 and 3 OTUs, respectively. Acetobacter pasteurianus was the predominant species in the genus Acetobacter, with a significantly higher detection rate (p<0.01) in Airag (11 samples) and Khoormog (5 samples) than in Tarag samples. The 13 OTUs in the genus Lactobacillus were assigned to 11 known species and 2 unknown species. In the 22 Airag samples, L. helveticus was predominant (at an average relative abundance of 85.9% across the 22 samples), and L. kefiranofaciens (2.4%, n = 20), Lactobacillus kefiri (1.2%, n = 21), Lactobacillus parakefiri (2.1%, n = 20) and Lactobacillus diolivorans (0.6%, n = 17) were dominant. In the 5 Khoormog samples, L. kefiranofaciens predominated (at an average relative abundance of 62.0%); L. helveticus (25.0%) and L. kefiri (4.5%) were dominant and detected in all 5 samples. Among the 26 Tarag samples, L. delbrueckii was the predominant species (74.6%, n = 24), whereas L. helveticus (19.2%, n = 24) and L. fermentum (2.7%, n = 11) were dominant. In the genus Lactococcus (Lc.), Lactococcus lactis was detected as the predominant species in Arag (1.3%, n = 13) and Khoormog (4.0%, n = 5), respectively, whereas S. thermophilus was predominant in Tarag samples (8.3%, n = 20).

Four OTUs showed less than 97% sequence similarity to any recognized species. These OTUs were assigned to Citricoccus sp., Lactobacillus sp. A, Lactobacillus sp. B and Lactococcus sp., respectively, because they showed similarity values of 96.4%, 92.2%, 93.1% and 94.8% with the type strains of their most closely related species―Citricoccus zhacaiensis FS24T (EU305672), Lactobacillus kalixensis Kx127A2T (AY253657), Lactobacillus kitasatonis JCM 1039T (AB107638) and Lactococcus chungangensis CAU 28T (EF694028), respectively (Fig. 1). Three of these four OTUs (the exception being Lactobacillus sp. B) were detected only in Airag samples.

Fig. 1.

Phylogenetic tree based on 16S rRNA gene sequences showing the taxonomic positions of the OTUs that showed less than 97% sequence similarities with their closest related species. The tree was constructed by using the neighbor-joining method on the basis of a comparison of approximately 300 bp; Escherichia coli ATCC 11775T served as an out-group. Bootstrap values (%) based on 1000 replicates are given at the nodes. Bar, 5% sequence divergence.

The OTUs that were not assigned to LAB species, such as Arthrobacter russicus, Bifidobacterium mongoliense, Brevundimonas nasdae, Clostridium perfringens, Enhydrobacter aerosaccus, Enterococcus durans, Enterococcus italicus, Macrococcus caseolyticus and Streptococcus parauberis, which are regarded as environmental-origin microbes, were detected in only a few samples and at low levels, with the exception of A. russicus and B. nasdae, which were detected in 14 Airag, 1 Khoormog and 9 Tarag samples and in 22 Airag, 5 Khoormog and 22 Tarag samples, respectively (Table 4).

Table 4.

Bacterial diversity of Airag, Khoormog and Tarag at the species level

Correlation between bacterial diversity and sample properties

The bacterial species diversity of each sample was analyzed based on species relative abundance, by using PCA, and compared with the differences among the animal species (mare, camel, cow, goat or yak) from which the milk was sourced and with the geographic differences of Airag, Khoormog and Tarag in 3 regions: the Gobi Desert region (provinces of Dundgobi and Umnugobi), the forest-steppe region (provinces of Arhangai and Uburhangai) and the steppe region (provinces of Bulgan and Tov) (Fig. 2). The relative abundance of L. delbrueckii, L. helveticus and L. kefiranofaciens had significant loadings to determine the sample plot location; relatively high loading was also observed for the relative abundance of S. thermophilus. Although the Tarag samples made from goat milk did not belong to any clusters, three clusters were found in the PCA plot on the basis of the animal species from which the milk was sourced (Fig. 2A). Cluster I comprised Airag samples made from mare’ s milk; L. helveticus was the significant factor in the formation of this cluster. Cluster II comprised Khoormog samples made from camel milk mainly, and L. kefiranofaciens was the main factor influencing the formation of this cluster. Cluster III comprised Tarag samples made from the milk of cows and yaks; L. delbrueckii and S. thermophilus were the main factors that contributed to the formation of this cluster. No correlations were found between sample bacterial diversities and the geographic sampling locations (Fig. 2B).

Fig. 2.

Two-dimensional scatter plot of PC1 and PC2 scores from the principal component analysis of bacterial species diversity. The sample plots were assigned on the basis of (A) the animal species from which the milk was sourced for the preparation of the fermented dairy products and (B) the geographic sampling locations. The four largest loads of bacterial species are shown with arrows. The contributions of PC1 (horizontal axis) and PC2 (vertical axis) are 76.1% and 21.1%, respectively.

DISCUSSION

In this study, we used a pyrosequencing method based on the 16S rRNA gene sequence to obtain detailed analyses of the bacterial diversity of traditional fermented dairy products of Mongolia. Based on the sequence of the V1−V2 hypervariable region, we were able to detect a total of 81 OTUs with an average of 13.9 ± 5.9 OTUs per sample (Table 1). An alpha diversity investigation suggested that the Airag and Khoormog samples were more diverse than the Tarag samples both in terms of OTU richness and phylogenetic relationships (Table 2). The Airag and Khoormog samples were made from non-pasteurized milk [3, 5], but the milk used for Tarag was boiled to remove milk fats [6], which would contribute to the differences in microbial diversity between Tarag and the other two products.

At the species level, a total of 22 LAB OTUs were observed, which showed high relative abundance in every sample, with the average being 98.5% ± 1.6%, suggesting that LAB were the core bacterial components involved in the fermentation of our samples. All of the LAB species that were isolated in our previous study that used the culture method [3] were also detected and identified with the pyrosequencing method in this study, except for Enterococcus faecium, Lactobacillus casei, Lactobacillus farciminis, Lactobacillus parafarraginis, Lactobacillus paraplantarum, Leuconostoc (Leuc.) pseudomesenteroides and Pediococcus pentosaceus. On the basis of sequence similarities, E. faecium, L. farciminis, L. paraplantarum and Leuc. pseudomesenteroides were not distinguishable from E. durans, Lactobacillus crustorum, Lactobacillus pentosus and Leuconostoc mesenteroides, respectively, in the V1−V2 hypervariable region with a 97% cutoff. The species name of L. casei was used based on the old taxonomy [28] in our previous study [3]; it should be corrected to Lactobacillus paracasei in accordance with the present taxonomy [29]. Thus, most species detected by using the culture method were also found with the pyrosequencing method. In addition, the pyrosequencing method detected 7 LAB species that were not detected by use of the culture method. These species could not be cultured in our previous study because of biases such as nutrient requirement, growth conditions or detection limit. Such biases may also explain the dissimilarity in LAB populations between the results obtained by culture and those obtained by pyrosequencing. The non-LAB OTUs were regarded as environmental contaminants from sources such as soils, animals and nomads [30,31,32,33,34]. Airag, Khoormog and Tarag were prepared by adding milk to a traditional container, such as the wooden barrel or cow’s skin bag, without the use of commercial starters. These containers were thought of as the main sources of the microbes involved in the fermentation of traditional fermented milks. The non-LAB OTUs, such as the species in the genera Clostridium, Enterococcus, Macrococcus and Moraxella, may have been transferred from the animals, because the milks for Airag and Khoormog were not heat-treated. Indeed, the average numbers of non-LAB OTUs in the Airag and Khoormog samples were higher than those in the Tarag samples. In this study, 4 OTUs did not show high 16S rRNA gene sequence similarities with known species, suggesting that these species were uncultured owing to specific culture requirements. In our previous study [35], we isolated two B. mongoliense species from two Airag samples at 5.5 × 107 cfu/ml and 4.0 × 107 cfu/ml, respectively, whereas in this study, B. mongoliense was found in only one Airag sample at a relative abundance of 0.2% (equal to about 106 cfu/ml). This result means that the pyrosequencing method only detected one-tenth (or less) of the bifidobacterial cells in the samples. This discrepancy between the present result and our previous result may be due to differences in the analytical methods used; in particular, the forward universal primer, 27F-mod, used in this study has mismatched bases for the 16S rRNA gene of Bifidobacterium [36], resulting in an underestimation of the bifidobacterial population by PCR.

PCA analysis based on the abundance of OTUs revealed an apparent correlation between the animal species providing the milk and the bacterial diversities of the Airag, Khoormog and Tarag samples. The abundance of the predominant LAB species made a major contribution to the formation of the clusters. Tarag samples made from goat milk did not show a clear PCA pattern, indicating that there are unknown factors influencing the microbial diversity of the LAB composition in these samples.

In conclusion, we used a pyrosequencing method to evaluate the traditional fermented dairy products of Mongolia (Airag, Khoormog and Tarag) and performed a detailed analysis of the microbial diversity of these products. We found a relationship between the predominant microbial component of these products and the animal species from which the milk was sourced for the preparation of the fermented dairy products. To our knowledge, this is the first study to apply a pyrosequencing method to these targets. Because pyrosequencing analysis was only able to reveal the microbial diversities at the time of sampling, further studies focused on manufacturing processes are needed to reveal the bacterial dynamics of these products.