Advertisement

Characterization and association of bacterial communities and nonvolatile components in spontaneously fermented cow milk at different geographical distances

  • Chuantao Peng
    Affiliations
    Key Laboratory of Dairy Biotechnology and Engineering, Ministry of Education, Inner Mongolia Agricultural University, Hohhot, Inner Mongolia, 010018, China

    Key Laboratory of Dairy Products Processing, Ministry of Agriculture and Rural Affairs, Inner Mongolia Agricultural University, Hohhot, Inner Mongolia, 010018, China
    Search for articles by this author
  • Zhihong Sun
    Affiliations
    Key Laboratory of Dairy Biotechnology and Engineering, Ministry of Education, Inner Mongolia Agricultural University, Hohhot, Inner Mongolia, 010018, China

    Key Laboratory of Dairy Products Processing, Ministry of Agriculture and Rural Affairs, Inner Mongolia Agricultural University, Hohhot, Inner Mongolia, 010018, China
    Search for articles by this author
  • Yaru Sun
    Affiliations
    Key Laboratory of Dairy Biotechnology and Engineering, Ministry of Education, Inner Mongolia Agricultural University, Hohhot, Inner Mongolia, 010018, China

    Key Laboratory of Dairy Products Processing, Ministry of Agriculture and Rural Affairs, Inner Mongolia Agricultural University, Hohhot, Inner Mongolia, 010018, China
    Search for articles by this author
  • Teng Ma
    Affiliations
    Key Laboratory of Dairy Biotechnology and Engineering, Ministry of Education, Inner Mongolia Agricultural University, Hohhot, Inner Mongolia, 010018, China

    Key Laboratory of Dairy Products Processing, Ministry of Agriculture and Rural Affairs, Inner Mongolia Agricultural University, Hohhot, Inner Mongolia, 010018, China
    Search for articles by this author
  • Weicheng Li
    Affiliations
    Key Laboratory of Dairy Biotechnology and Engineering, Ministry of Education, Inner Mongolia Agricultural University, Hohhot, Inner Mongolia, 010018, China

    Key Laboratory of Dairy Products Processing, Ministry of Agriculture and Rural Affairs, Inner Mongolia Agricultural University, Hohhot, Inner Mongolia, 010018, China
    Search for articles by this author
  • Heping Zhang
    Correspondence
    Corresponding author
    Affiliations
    Key Laboratory of Dairy Biotechnology and Engineering, Ministry of Education, Inner Mongolia Agricultural University, Hohhot, Inner Mongolia, 010018, China

    Key Laboratory of Dairy Products Processing, Ministry of Agriculture and Rural Affairs, Inner Mongolia Agricultural University, Hohhot, Inner Mongolia, 010018, China
    Search for articles by this author
Open ArchivePublished:January 14, 2021DOI:https://doi.org/10.3168/jds.2020-19303

      ABSTRACT

      In the ecosystem of spontaneously fermented cow milk, the characteristics and relationship of bacterial communities and nonvolatile components at different scales of geographical distances (provincial, county, and village levels) are unclear. Here, 25 sampling sites from Xin Jiang and Tibet, 2 provinces of China, were selected based on the distribution of spontaneously fermented cow milk and used for metagenomic and metabolomic analysis. At the provincial geographical distance, the same predominant species, Lactobacillus delbrueckii ssp. bulgaricus and Streptococcus thermophilus, were detected in Xin Jiang and Tibet. Further, the richness of the bacterial composition of samples from Tibet was higher than those from Xin Jiang; specifically, at the species level, 28 species were identified in Tibet samples but only 7 species in Xin Jiang samples. At the provincial geographical level, we detected significant differences in bacterial structure, shown in principal coordinate analysis plots, and significant differences (Simpson index) in bacterial diversity were also detected. However, at the county and village levels, no significant differences were detected in bacterial communities and diversity, but a difference in bacterial compositions was detectable. This indicates that bacterial communities and diversity of spontaneously fermented milk dissimilarity significantly increased with geographic distance. For the nonvolatile component profiles, the partial least squares discriminant analysis plot (R2Y > 0.5 and Q2 > 0.5 for the goodness-of-fit and predictive ability parameter, respectively) showed that samples from different geographical distances (provincial, county, and village) were all separated, which indicated that all the discriminations in nonvolatile components profiles were from different geographical distances. Investigating relationships between lactic acid bacteria and discriminatory nonvolatile components at the county level showed that 9 species were positively correlated with 16 discriminatory nonvolatile components, all species with low abundance rather than the predominant species L. delbrueckii ssp. bulgaricus and Strep. thermophilus, which indicates the importance of the selection of autochthonous nonpredominant bacteria.

      Key words

      INTRODUCTION

      Spontaneous milk fermentation has a long history, and the traditional methods and processes have been handed down through generations (
      • Yu J.
      • Wang W.H.
      • Menghe B.L.
      • Jiri M.T.
      • Wang H.M.
      • Liu W.J.
      • Bao Q.H.
      • Lu Q.
      • Zhang J.C.
      • Wang F.
      • Xu H.Y.
      • Sun T.S.
      • Zhang H.P.
      Diversity of lactic acid bacteria associated with traditional fermented dairy products in Mongolia.
      ). Many traditional fermented dairy products are regarded as nutraceuticals and have typical organoleptic features that depend on local and regional traditions, environmental factors, and indigenous microorganisms (
      • Corsetti A.
      • Settanni L.
      • Lopez C.C.
      • Felis G.E.
      • Mastrangelo M.
      • Suzzi G.
      A taxonomic survey of lactic acid bacteria isolates from wheat (Triticum durum) kernels and non-conventional flours.
      ;
      • Turbes G.
      • Linscott T.D.
      • Tomasino E.
      • Waite-Cusic J.
      • Lim J.
      • Meunier-Goddik L.
      Evidence of terroir in milk sourcing and its influence on Cheddar cheese.
      ). Spontaneous milk fermentation practices encourage or rely entirely on “native” microbiota to carry out fermentations, and the indigenous microorganisms are considered to play a major role in determining regional typicity (
      • Mathara J.M.
      • Schillinger U.
      • Kutima P.M.
      • Mbugua S.K.
      • Holzapfel W.H.
      Isolation, identification and characterisation of the dominant microorganisms of kule naoto: The Maasai traditional fermented milk in Kenya.
      ).
      Milk is rapidly colonized following milking by a variety of other microbes coming from the teat canal, udder skin, milking machines, and containers used for milk storage, which can reflect the farm and pasture environment (
      • Coorevits A.
      • De Jonghe V.
      • Vandroemme J.
      • Reekmans R.
      • Heyrman J.
      • Messens W.
      • De Vos P.
      • Heyndrickx M.
      Comparative analysis of the diversity of aerobic spore-forming bacteria in raw milk from organic and conventional dairy farms.
      ;
      • Vacheyrou M.
      • Normand A.C.
      • Guyot P.
      • Cassagne C.
      • Piarroux R.
      • Bouton Y.
      Cultivable microbial communities in raw cow milk and potential transfers from stables of sixteen French farms.
      ;
      • Quigley L.
      • O'Sullivan O.
      • Stanton C.
      • Beresford T.P.
      • Ross R.P.
      • Fitzgerald G.F.
      • Cotter P.D.
      The complex microbiota of raw milk.
      ;
      • Addis M.F.
      • Tanca A.
      • Uzzau S.
      • Oikonomou G.
      • Bicalho R.C.
      • Moroni P.
      The bovine milk microbiota: Insights and perspectives from -omics studies.
      ). Also, the bacterial communities of spontaneously fermented milk vary geographically, and regional environmental factors and other abiotic factors are important in shaping the bacterial communities (
      • Parker M.
      • Zobrist S.
      • Donahue C.
      • Edick C.
      • Mansen K.
      • Nadjari M.H.Z.
      • Heerikhuisen M.
      • Sybesma W.
      • Molenaar D.
      • Diallo A.M.
      • Milani P.
      • Kort R.
      Naturally fermented milk from northern Senegal: Bacterial community composition and probiotic enrichment with lactobacillus rhamnosus..
      ). A large number of studies have reported the diversity of the bacterial communities in spontaneously fermented milk from different regions, such as Russia, Mongolia, Inner Mongolia, and Xin Jiang (China), and differences have been found in bacterial compositions based on different regions (
      • Sun Z.
      • Liu W.
      • Bao Q.
      • Zhang J.
      • Hou Q.
      • Kwok L.
      • Sun T.
      • Zhang H.
      Investigation of bacterial and fungal diversity in tarag using high-throughput sequencing.
      ;
      • Liu W.
      • Zheng Y.
      • Kwok L.Y.
      • Sun Z.
      • Zhang J.
      • Guo Z.
      • Hou Q.
      • Menhe B.
      • Zhang H.
      High-throughput sequencing for the detection of the bacterial and fungal diversity in Mongolian naturally fermented cow's milk in Russia.
      ;
      • Gesudu Q.
      • Zheng Y.
      • Xi X.X.
      • Hou Q.C.
      • Xu H.Y.
      • Huang W.Q.
      • Zhang H.P.
      • Menghe B.
      • Liu W.J.
      Investigating bacterial population structure and dynamics in traditional koumiss from Inner Mongolia using single molecule real-time sequencing.
      ). However, these studies may not meticulously explore microbial communities or metabolic profiles of one region at different geographical distances.
      Recently, some researches have demonstrated metabolic characteristics of milk and its derivative products also linked to the territory of production, and have provided evidence of terroir in milk (
      • Parker M.
      • Zobrist S.
      • Donahue C.
      • Edick C.
      • Mansen K.
      • Nadjari M.H.Z.
      • Heerikhuisen M.
      • Sybesma W.
      • Molenaar D.
      • Diallo A.M.
      • Milani P.
      • Kort R.
      Naturally fermented milk from northern Senegal: Bacterial community composition and probiotic enrichment with lactobacillus rhamnosus..
      ). Microbial activity is an inner driving force of spontaneous milk fermentation; however, few reports exist on the relationship between regional bacterial communities and metabolic profiles of spontaneously fermented milk. Furthermore, it is presently unknown what type of bacterial species, predominant species, or bacterial species with low abundance determine regional metabolic characteristics. Therefore, it is of interest to further characterize bacterial communities and metabolic profiles at different geographical distances and then build associations between them.
      The geography and traditional way of life of Xin Jiang and Tibet offer advantages for livestock farming, and local nomadic families have a rich history of making traditional dairy products. In this study, 25 sampling sites from Xin Jiang and Tibet at different scales of geographical distances (provincial, county, and village levels) were selected. To understand bacterial communities and nonvolatile component characteristics in spontaneously fermented cow milk at different geographical distances, we conducted an exploratory study to assess (1) whether the bacterial communities and nonvolatile components exhibit difference at different scales of geographical distances: province (distance above 1,000 km), county (distance above 100 km), and village (distance within 10 km); and (2) whether predominant species or bacterial species with low abundance determine regional nonvolatile component characteristics. A simplified schematic of the experimental approach used to study characterization and associations of bacterial communities and nonvolatile components in spontaneously fermented cow milk from different geographical distances is shown in Figure 1.
      Figure thumbnail gr1
      Figure 1Simplified schematic of the experimental approach used to study characterization and association of bacterial community and metabolic profiles in spontaneously fermented cow milks at different scales of geographical distances: province (distance above 1,000 km), county (distance above 100 km), and village (distance within 10 km). Spontaneously fermented cow milk samples were collected from 25 sampling sites across 6 counties (TK = Tashenkuergan, YC = Yecheng, HT = Hetian, BM = Bomin, JD = Jiangda, LZ = Linzhi) in Xin Jiang and Tibet.

      MATERIALS AND METHODS

      Sampling of Spontaneously Fermented Cow Milk

      In total, 25 sampling sites from 6 counties in China (TK = Tashenkuergan, YC = Yecheng, HT = Hetian, BM = Bomin, JD = Jiangda, LZ = Linzhi) were selected based on the distribution of spontaneously fermented cow milk in Xin Jiang and Tibet, covering an area of approximately 2.9 million km2, and the sampling sites from each county correspond to several adjacent farms. Specifically, 9 sampling sites were from YC, 3 sites from HT, 3 sites from TK, 3 sites from JD, 4 sites from BM, and 3 sites from LZ. The geographical distances between sampling sites from TK, YC, and HT in Xin Jiang and those from JD, LZ, and BM in Tibet were above 1,000 km; the geographical distances between sampling sites from TK, YC, and HT in Xin Jiang and between the sampling sites from JD, LZ, and BM in Tibet were above 100 km. The geographical distances of sampling sites within one county were below 10 km. Farmers adopted similar procedures for making spontaneously fermented cow milk in Xin Jiang and Tibet. Briefly, the milk used for spontaneous fermentation was locally produced in the sampling sites. Fresh unpasteurized milk is poured into a large custom-made, specially treated jar made of tung wood and is fermented spontaneously without any inoculation of starter cultures at room temperature of around 15 to 30°C. After 2 to 3 d, the desirable curdy texture and flavor will be achieved, and then the products are sampled. A total of 75 samples (3 replicate samples from each site) were collected from the 25 sampling sites across 6 counties in Xin Jiang (n = 45 sites) and Tibet (n = 30) between July and August in 2019. Cows had access to pasture during the sampling period. Approximately 50 mL of samples from each sampling site were collected aseptically and transported in a tank of liquid nitrogen. The samples were then transported to the laboratory under refrigerated conditions and were stored at −80°C until further analysis.

      Single Molecule, Real-Time Sequencing and Bioinformatic Analysis

      The samples were thawed in a water bath at 37°C until they were completely melted and then uniformly shaken for 3 min to mix well. For each sample, 3 g of fermented milk were centrifuged at 10,000 × g for 20 min in an Eppendorf 5810R centrifuge (Eppendorf AG, Hamburg, Germany). The supernatant was then removed, and the remaining pellet was resuspended in 5 mL of sterile water; 500 mg of lysozyme (Sigma-Aldrich, Gallarate, Italy) was added to improve the extraction of bacterial DNA and then incubated at 37°C for 12 h. Genomic DNA was extracted from each sample using the DNeasy Mericon Food Kit (69514, Qiagen, Hilden, Germany) according to the manufacturer's instructions. The quality of extracted DNA was evaluated by 0.8% agarose gel electrophoresis and a NanoDrop ND-1000 spectrophotometer (Thermo Fisher Scientific, Wilmington, DE).
      The full-length sequences of 16S rRNA gene were amplified by PCR for single molecule, real-time (SMRT) sequencing. The primers and reaction conditions were performed as described previously (
      • Hou Q.
      • Xu H.
      • Zheng Y.
      • Xi X.
      • Kwok L.
      • Sun Z.
      • Zhang H.
      • Zhang W.
      Evaluation of bacterial contamination in raw milk, ultra-high temperature milk and infant formula using single molecule, real-time sequencing technology.
      ), and the quality control for PCR amplification and sequence preprocessing were performed as described previously (
      • Quigley L.
      • O'Sullivan O.
      • Beresford T.P.
      • Ross R.P.
      • Fitzgerald G.F.
      • Cotter P.D.
      High-throughput sequencing for detection of subpopulations of bacteria not previously associated with artisanal chesses.
      ;
      • Yu J.
      • Mo L.
      • Pan L.
      • Yao C.
      • Ren D.
      • An X.
      • Tsogtgerel T.
      • Zhang H.
      • Liu W.
      Bacterial microbiota and metabolic character of traditional sour cream and butter in Buryatia, Russia.
      ). The amplicons were sequenced using the P6-C4 chemistry of the PacBio RS II instrument (Pacific Biosciences, Menlo Park, CA), and purified PCR products were used to construct DNA libraries with the Pacific Biosciences Template Prep Kit 2.0. The RS_ReadsOfinsert.1 protocol available on the SMRT Portal (v. 2.7, https://smrt-analysis.readthedocs.io/en/latest/SMRT-Pipe-Reference-Guide-v2.2.0/) was used to process raw sequence data, and, briefly, the restrictive filtering parameters of the raw data were based on the following criteria: (1) minimum full passes of up to 5; (2) minimum predicted accuracy of 90; (3) minimum read length of inserts of 1,400; and (4) maximum read length of inserts of 1,800 (
      • Hou Q.
      • Xu H.
      • Zheng Y.
      • Xi X.
      • Kwok L.
      • Sun Z.
      • Zhang H.
      • Zhang W.
      Evaluation of bacterial contamination in raw milk, ultra-high temperature milk and infant formula using single molecule, real-time sequencing technology.
      ;
      • Gesudu Q.
      • Zheng Y.
      • Xi X.X.
      • Hou Q.C.
      • Xu H.Y.
      • Huang W.Q.
      • Zhang H.P.
      • Menghe B.
      • Liu W.J.
      Investigating bacterial population structure and dynamics in traditional koumiss from Inner Mongolia using single molecule real-time sequencing.
      ).
      After removal of barcodes and primers, high-quality sequences were extracted by the Quantitative Insights into Microbial Ecology (QIIME) package (v. 1.7, https://github.com/biocore/qiime). Briefly, the most abundant sequence from each cluster was selected as representative sequence and aligned and obtained using the PyNAST and UCLUST systems with 100% clustering of sequence identity (
      • Caporaso J.G.
      • Bittinger K.
      • Bushman F.D.
      • DeSantis T.Z.
      • Andersen G.L.
      • Knight R.
      PyNAST: A flexible tool for aligning sequences to a template alignment.
      ;
      • Edgar R.C.
      Search and clustering orders of magnitude faster than BLAST.
      ). After selection of representative sequences, sequences were classified into operational taxonomic units (OTU) below a threshold of 97% using UCLUST software (
      • Lozupone C.
      • Hamady M.
      • Knight R.
      UniFrac—An online tool for comparing microbial community diversity in a phylogenetic context.
      ). Chimer-aSlayer (
      • Haas B.J.
      • Gevers D.
      • Earl A.M.
      • Feldgarden M.
      • Ward D.V.
      • Giannoukos G.
      • Ciulla D.
      • Tabbaa D.
      • Highlander S.K.
      • Sodergren E.
      Chimeric 16S rRNA sequence formation and detection in Sanger and 454-pyrosequenced PCR amplicons.
      ) was applied to remove potential chimeric sequences in the representative OTU set. Subsequently, the taxonomy of each OTU representative sequence was assigned using the Ribosomal Database Project II database at a minimum bootstrap threshold of 80% (
      • Cole J.R.
      • Chai B.
      • Farris R.J.
      • Wang Q.
      • Kulam-Syed-Mohideen A.S.
      • McGarrell D.M.
      • Bandela A.M.
      • Cardenas E.
      • Garrity G.M.
      • Tiedje J.M.
      The ribosomal database project (RDP-II): Introducing myRDP space and quality controlled public data.
      ). The principal coordinate analysis (PCoA) using the Bray-Curtis distance of species abundance was performed. Permutational multivariate analysis of variance (PERMANOVA) was used to evaluate differences in the structure and composition of bacteria between different samples. The permutations helped avoid possible biases. Shannon-Wiener and Simpson indices were calculated to evaluate bacterial α-diversity. Bacteria with lower than 0.1% of total abundance and unclassified bacteria were not discussed.

      Determination of Nonvolatile Components by Ultra-Performance Liquid Chromatography with Quadrupole Time-of-Flight MS

      In preparing samples for analyzing nonvolatile component profiles, for each sample, 2 mL of the spontaneously fermented milk was centrifuged at 4,000 × g for 10 min. The fat layer was discarded, and 1 mL of the supernatant was collected and transferred to a fresh tube. Seven milliliters of acetonitrile solution was added to the supernatant, and the mixture was vortexed for 2 min at room temperature. Then, the mixture was centrifuged at 12,000 × g for 10 min again, and the supernatant was collected and transferred to a fresh tube. The supernatant was concentrated to near-dryness in a vacuum concentrator and resuspended with 500 μL of 40% (vol/vol) acetonitrile solution. The resuspended samples were filtered through a 0.22-μm water-insoluble microporous membrane and then stored at −20°C for subsequent analysis by liquid chromatography-MS.
      The pretreated milk samples were loaded onto an Acquity ultra-performance liquid chromatography system coupled to a Waters quadrupole time-of-flight mass spectrometer (Waters Corp., Milford, MA). Injection volume was 4 μL per sample, and each treatment was performed in triplicate. Chromatographic separations were performed on an Acquity HSS T3 column (2.1 mm × 100 mm × 1.8 μm; Waters Corp.), the flow rate was 0.3 mL/min, and the column heater was set to 40°C throughout the analysis. An electrospray ionization source that operated in both positive and negative ion modes was used to generate ions for analysis. The ion scan mode, liquid phase, and mass spectrometry conditions were selected according to a previous study (
      • Xi X.M.
      • Kwok L.Y.
      • Wang Y.N.
      • Ma C.
      • Mi Z.H.
      • Zhang H.P.
      Ultra-performance liquid chromatography-quadrupole-time of flight mass spectrometry MSE-based untargeted milk metabolomics in dairy cows with subclinical or clinical mastitis.
      ). To ensure mass reproducibility and accuracy, leucine-enkephalin was used as the lock mass for positive ion mode [556.2771(M + H)+] and negative ion mode [554.2615(M − H)].
      Raw data obtained via ultra-performance liquid chromatography with quadrupole time-of-flight MSE, where E represents elevated energy mode of acquisition, were collected in continuum mode and processed using the Progenesis QI software (Waters Corp.). Data preprocessing procedures included peak alignment, peak identification, and deconvolution. Subsequently, data were further processed, and compounds showing significant difference between different regional samples were identified as subsequently described in detail.

      Multivariate Statistical Analysis

      All statistical analysis was conducted using R v. 3.2.4 (2016, R Foundation for Statistical Computing, Vienna, Austria) or MetaboAnalyst (https://www.metaboanalyst.ca). For the preprocessed liquid chromatography-MS data, multivariate statistical analysis such as principal coordinate analysis, and partial least squares discriminant analysis (PLS-DA) was then performed on the data set using MetaboAnalyst. The PLS-DA was used to visualize discrimination and differences of nonvolatile component profiles between samples. Before multivariate statistical analysis, all data were mean-centered and Pareto-scaled. One-tail Student's t-test (P < 0.05) was used to assess statistical significance between fermented milk samples, and a variable importance in projection score is a measure of a variable's importance in the PLS-DA model. For variable selections, the nonvolatile components with P < 0.05 and variable importance in projection >1 were selected as significantly different nonvolatile components.
      According to retention time, mass-to-charge ratio, and fragment information, the selected differentiating components were identified using the Progenesis QI software, which was set to automatically search the selected databases such as Human Metabolome Database (https://www.hmdb.ca), Kyoto Encyclopedia of Genes and Genomes (https://www.genome.jp/kegg/), ChemSpider (https://www.chemspider.com), METLIN (https://metlin.scripps.edu/), and Yeast Metabolome Database (http://www.ymdb.ca/) databases.

      RESULTS

      Characterization and Difference of Bacterial Communities and Nonvolatile Component Profiles at Different Geographical Distances

      A source map of bacterial communities driving spontaneous cow milk fermentation at the provincial geographical level (distance above 1,000 km) for Xin Jiang and Tibet is displayed in Figure 2. Overall, we identified 1 order, 4 families, 6 genera, and 28 species in Xin Jiang and Tibet samples (Supplemental File S1, https://doi.org/10.6084/m9.figshare.13297712). In both provinces, we found the same predominant species: Lactobacillus delbrueckii ssp. bulgaricus and Streptococcus thermophilus. In Xin Jiang, these 2 species constitute, on average, 55 and 39%, respectively, of the total bacterial communities; in Tibet, 68 and 12%, respectively. However, the richness of the samples from Tibet was higher than that from Xin Jiang, and Simpson index (P = 0.012) for Tibet was significantly higher than that for Xin Jiang (Supplemental Figure S1, https://doi.org/10.6084/m9.figshare.13297709); specifically, at the species level, 28 species were identified in Tibetan samples, and only 7 species in samples from Xin Jiang. Also, we performed PCoA using Bray-Curtis distances based on the OTU abundance, shown in Figure 3a, and the results revealed significant difference (P = 0.008) in bacterial communities. In the present study, it should be noted that some non-lactic acid bacteria (LAB) and unclassified bacteria were identified, such as Staphylococcus chromogenes, but the abundances of these non-LAB were extremely low (less than 0.1% of total abundance). As previously described, the bacteria with abundance less than 0.1% of total abundance, as well as unclassified bacteria, were not discussed. We next sought to test whether differentiable nonvolatile component profiles exist at the provincial geographical level. We analyzed the nonvolatile component profiles of 30 samples from 10 Xin Jiang sampling sites and 30 samples from 10 Tibet sampling sites in triplicate (Supplemental File S2, https://doi.org/10.6084/m9.figshare.13297715) using PLS-DA. The PLS-DA plot showed that nonvolatile component profiles of fermented milk samples from Xin Jiang and Tibet were separated clearly, which indicated discrimination in nonvolatile components profiles between Xin Jiang and Tibet (Figure 3b). The R2Y and Q2 values, representing the goodness-of-fit and the predictive ability parameter, were 0.88 and 0.81, respectively, which were relatively high, indicating that the model had good accuracy and predictability.
      Figure thumbnail gr2
      Figure 2Source map of bacterial communities at the order, family, genus, and species levels, respectively, derived from spontaneous cow milk fermentation at the provincial level (distance above 1,000 km) in Xin Jiang and Tibet. The diameter of each circle is proportional to relative abundance.
      Figure thumbnail gr3
      Figure 3Characteristics and differences of bacterial community performed by principal coordinate analysis (PCoA) using Bray-Curtis distances and permutational multivariate analysis of variance (PERMANOVA; P = 0.012; a) and metabolic profiles performed via partial least squares discriminant analysis (b) at the provincial level (distance greater than 1,000 km) in Xin Jiang and Tibet. t[1] = weight of regression coefficient of the predictive principal component; t[2] = weight of the regression coefficient of the orthogonal principal component.
      At relatively smaller geographical distance, the county geographical distance (above 100 km), we found no significant differences (P > 0.05) in bacterial structure shown in the PCoA plots, and also no significant differences were detectable [Shannon index (P > 0.05); Simpson index (P > 0.05)] in bacterial diversity (Supplemental Figure S2, https://doi.org/10.6084/m9.figshare.13297709). However, the clusters of bacterial communities were not irregular, and, to some extent, the distributions were conditioned by the geographical distances between the 6 counties (Figure 4a). For the nonvolatile component profiles, the PLS-DA plot showed that samples from the 6 counties were separated, which indicated discrimination in nonvolatile components profiles among the counties (Figure 4b). The R2Y and Q2 values were 0.96 and 0.86, respectively. Also, the distributions of nonvolatile component profiles of fermented milk samples in the PLS-DA plot were conditioned by the geographical distances between the 6 counties (Figure 4b).
      Figure thumbnail gr4
      Figure 4Characteristics and differences of bacterial community performed by principal coordinate analysis (PCoA) using Bray-Curtis distances and permutational multivariate analysis of variance (PERMANOVA; P > 0.05; a) and metabolic profiles performed via partial least squares discriminant analysis (b) at the county level (distance greater than 100 km) in Bomin (BM), Jiangda (JD), and Linzhi (LZ) from Tibet, and Tashenkuergan (TK), Yecheng (YC), and Hetian (HT) from Xin Jiang. t[1] = weight of regression coefficient of the predictive principal component; t[2] = weight of the regression coefficient of the orthogonal principal component.
      At the smallest geographical distance, we sought to test whether bacterial communities and nonvolatile component profiles can be distinguished between 2 villages (distance within 10 km) within 1 county (YC). No significant differences (P > 0.05) in bacterial structure were shown in the PCoA plots (Supplemental Figure S3, https://doi.org/10.6084/m9.figshare.13297709), and also we detected no significant differences [Shannon index (P = 0.18); Simpson index (P = 0.43)] in bacterial diversity (Supplemental Figure S4, https://doi.org/10.6084/m9.figshare.13297709). However, a difference in bacterial composition was detectable. In the village of Kekeya, we found 2 predominant species, L. delbrueckii ssp. bulgaricus and Strep. thermophilus, which constitute, on average, 48.5 and 49.6%, respectively, of the total bacterial communities. In the village of Xinbashen, we detected 3 predominant species: L. delbrueckii ssp. bulgaricus, Strep. thermophilus, and Lactobacillus helveticus, which constitute, on average, 45, 32.9, and 20.7% of the total bacterial communities, respectively. The richness of the samples from Xinbashen was higher than that of the samples from Kekeya; specifically, at the species level, 14 species were identified in Xinbashen samples and only 7 species in Kekeya samples (Figure 5a). The PLS-DA plot showed that nonvolatile component profiles between samples from the 2 villages were separated clearly, which indicated discrimination in nonvolatile components profiles between samples from the 2 villages (Figure 5b). The R2Y and Q2 values were 0.98 and 0.93, respectively. The results indicate that even in adjacent villages (within 10 km) differences in bacterial composition and nonvolatile component profiles still exist.
      Figure thumbnail gr5
      Figure 5Characteristics and differences of bacterial community and percentages of bacterial species detected in samples (shown in pie charts, a) and metabolic profiles performed by partial least squares discriminant analysis (b) at the village level (distance within 10 km) in Kekeya and Xinbashen, 2 villages of Yecheng (YC), China. t[1] = weight of regression coefficient of the predictive principal component; t[2] = weight of the regression coefficient of the orthogonal principal component.

      Correlations Between Discriminatory Nonvolatile Components and Bacterial Species at the County Level

      The discriminatory nonvolatile components that can distinguish different counties were analyzed and identified by variable importance in projection and P-value, and 22 discriminatory nonvolatile components were identified (Supplemental File S2, https://doi.org/10.6084/m9.figshare.13297715). We next generated a heat-map of discriminatory nonvolatile components between 6 different counties and analyzed them via hierarchal clustering to visualize county-specific discriminatory nonvolatile components (Figure 6a).
      Figure thumbnail gr6
      Figure 6Heat-map and hierarchal clustering of discriminatory metabolites between 6 different counties (a). Correlation network analysis was used to investigate underlying relationships between bacterial community and discriminatory metabolites at the county level (distance above 100 km; b). TK = Tashenkuergan, YC = Yecheng, HT = Hetian, BM = Bomin, JD = Jiangda, LZ = Linzhi. GADGM = galactopyranosyl-(1->4)-2-amino-2-deoxy-beta-d-glucopyranosyl-(1->6)-d-mannose.
      To further dissect the relationship between LAB and nonvolatile components of spontaneously fermented cow milk, correlation network analysis was used to investigate underlying relationships between LAB and 22 discriminatory nonvolatile components at the county level. The false discovery rate method (false discovery rate < 0.05) was used to define statistical significance. Nine species were positively correlated with 16 discriminatory nonvolatile components, and no negative correlations were detected, as demonstrated through Spearman correlation coefficients (Figure 6b). Surprisingly, we detected no correlations between discriminatory nonvolatile components and the 2 predominant species [L. delbrueckii ssp. bulgaricus (relative abundance >40%) and Strep. thermophilus (relative abundance >20%)].
      Among these 9 species, Limosilactobacillus reuteri, Streptococcus pasteurianus, Streptococcus anginosus, and Vagococcus fluvialis were the top 4 species, indicating that these 4 species with low abundance had a great positive influence on the nonvolatile component characteristics of spontaneously fermented cow milk from 6 different counties. For example, L. reuteri was positively correlated with 5 discriminatory nonvolatile components: isoleucylphenylalanine, l-isoleucine, 1-piperidinecarboxaldehyde, prolyltyrosine, and phenylalanylproline. Among samples from the 6 counties, the abundance of L. reuteri was the highest in fermented milk samples from BM, and, correspondingly, the relative concentrations of the 5 discriminatory nonvolatile components were also higher in BM fermented milk samples than in those from other counties. This indicates high correlations between these 9 bacterial species and the 16 discriminatory nonvolatile components.

      DISCUSSION

      In previous research, the main focus has been on microbial communities or nonvolatile component profiles of spontaneously fermented dairy products from a single region, and the collected samples may be from one sampling site or random sampling sites within a region (
      • Mathara J.M.
      • Schillinger U.
      • Kutima P.M.
      • Mbugua S.K.
      • Holzapfel W.H.
      Isolation, identification and characterisation of the dominant microorganisms of kule naoto: The Maasai traditional fermented milk in Kenya.
      ;
      • Liu W.
      • Zheng Y.
      • Kwok L.Y.
      • Sun Z.
      • Zhang J.
      • Guo Z.
      • Hou Q.
      • Menhe B.
      • Zhang H.
      High-throughput sequencing for the detection of the bacterial and fungal diversity in Mongolian naturally fermented cow's milk in Russia.
      ;
      • Gesudu Q.
      • Zheng Y.
      • Xi X.X.
      • Hou Q.C.
      • Xu H.Y.
      • Huang W.Q.
      • Zhang H.P.
      • Menghe B.
      • Liu W.J.
      Investigating bacterial population structure and dynamics in traditional koumiss from Inner Mongolia using single molecule real-time sequencing.
      ;
      • Yu J.
      • Mo L.
      • Pan L.
      • Yao C.
      • Ren D.
      • An X.
      • Tsogtgerel T.
      • Zhang H.
      • Liu W.
      Bacterial microbiota and metabolic character of traditional sour cream and butter in Buryatia, Russia.
      ). However, these studies may not exactly reflect microbial communities or truly nonvolatile components of one region at different geographical distances. In the present study, from different scales of geographical distances [province (distance above 1,000 km), county (distance above 100 km), and village (distance within 10 km)], the characteristics and association of bacterial community and nonvolatile components in spontaneously fermented cow milk were analyzed. Also, it should be noted that this research mainly focused on bacterial composition analyzed by 16S rRNA sequencing and nonvolatile components using liquid chromatography-MS. Because of resistance to low pH and antimicrobial compound production, spontaneously fermented cow milk is mainly dominated by LAB.
      At the provincial geographical distance (greater than 1,000 km) in Xin Jiang and Tibet, we found the same predominant species: L. delbrueckii ssp. bulgaricus and Strep. thermophilus. In Mongolia, the predominant species of spontaneously fermented milk samples from 13 different regions were Strep. thermophilus and L. helveticus (
      • Yu J.
      • Wang W.H.
      • Menghe B.L.
      • Jiri M.T.
      • Wang H.M.
      • Liu W.J.
      • Bao Q.H.
      • Lu Q.
      • Zhang J.C.
      • Wang F.
      • Xu H.Y.
      • Sun T.S.
      • Zhang H.P.
      Diversity of lactic acid bacteria associated with traditional fermented dairy products in Mongolia.
      ). Lactobacillus delbrueckii ssp. bulgaricus and Strep. thermophilus were reported as the most frequently isolated bacteria in Columbian dairy products (
      • Perea Vélez M.
      • Hermans K.
      • Verhoeven T.L.A.
      • Lebeer S.E.
      • Vanderleyden J.
      • De Keersmaecker S.C.J.
      Identification and characterization of starter lactic acid bacteria and probiotics from Columbian dairy products.
      ). This indicates that different regions may have different predominant species in spontaneously fermented milk products. In addition, the bacterial composition of the spontaneously fermented cow milk samples from Tibet was much more diverse than that of the samples from Xin Jiang. The microbial community of cow teat skin depends on the dairy farm environment, such as forage plants and housing conditions (
      • Zdanowicz M.
      • Shelford J.A.
      • Tucker C.B.
      • Weary D.M.
      • von Keyserlingk M.A.G.
      Bacterial populations on teat ends of dairy cows housed in free stalls and bedded with either sand or sawdust.
      ;
      • Verdier-Metz I.
      • Gagne G.
      • Bornes S.
      • Monsallier F.
      • Veisseire P.
      • Delbès-Paus C.
      • Montel M.C.
      Cow teat skin, a potential source of diverse microbial populations for cheese production.
      ). The teat skin of cows is considered a major reservoir of microbial diversity for raw milk and dairy products, and the microbiota of raw milk has a pivotal role in dairy product fermentation (
      • Verdier-Metz I.
      • Gagne G.
      • Bornes S.
      • Monsallier F.
      • Veisseire P.
      • Delbès-Paus C.
      • Montel M.C.
      Cow teat skin, a potential source of diverse microbial populations for cheese production.
      ;
      • Doyle C.J.
      • Gleeson D.
      • O'Toole P.W.
      • Cotter P.D.
      Impacts of seasonal housing and teat preparation on raw milk microbiota: A high-throughput sequencing study.
      ). Some previous studies have reported that extremes of temperature, UV exposure, and low humidity might decrease microbial diversity of plants (
      • Whipps J.M.
      • Hand P.
      • Pink D.
      • Bending G.D.
      Phyllosphere microbiology with special reference to diversity and plant genotype.
      ;
      • Redford A.J.
      • Bowers R.M.
      • Knight R.
      • Linhart Y.
      • Fierer N.
      The ecology of the phyllosphere: Geographic and phylogenetic variability in the distribution of bacteria on tree leaves.
      ). Therefore we speculate that the large differences in temperature between day and night, low humidity, long daylight hours, and strong solar radiation of Xin Jiang may reduce microbial diversity of plants and thus decrease the microbial diversity of fermented milk samples from Xin Jiang.
      At the provincial geographical level, bacterial communities and diversity of fermented milk samples in Xin Jiang were significantly different from those of samples from Tibet. However, at either the county or the village geographical level, we detected no significant differences in bacterial communities and diversity, despite differences in bacterial compositions. This indicates that bacterial communities and diversity of spontaneously fermented milk dissimilarity significantly increased with geographic distance. This phenomenon has also been found in plant and wine microbial research (
      • Bokulich N.A.
      • Thorngate J.H.
      • Richardson P.M.
      • Mills D.A.
      Microbial biogeography of wine grapes is conditioned by cultivar, vintage, and climate.
      ;
      • Wang Y.L.
      • Gao C.
      • Chen L.
      • Ji N.N.
      • Wu B.W.
      • Li X.C.
      • Lü P.P.
      • Zheng Y.
      • Guo L.D.
      Host plant phylogeny and geographic distance strongly structure Betulaceae-associated ectomycorrhizal fungal communities in Chinese secondary forest ecosystems.
      ). Such studies have concluded that this could be because geographic isolation generates dispersal barriers and reduces migration rates of bacterial and fungi propagules (
      • Peay K.G.
      • Garbelotto M.
      • Bruns T.D.
      Evidence of dispersal limitation in soil microorganisms: Isolation reduces species richness on mycorrhizal tree islands.
      ;
      • Peay K.G.
      • Schubert M.G.
      • Nguyen N.H.
      • Bruns T.D.
      Measuring ectomycorrhizal fungal dispersal: Macroecological patterns driven by microscopic propagules.
      ;
      • Wang Y.L.
      • Gao C.
      • Chen L.
      • Ji N.N.
      • Wu B.W.
      • Li X.C.
      • Lü P.P.
      • Zheng Y.
      • Guo L.D.
      Host plant phylogeny and geographic distance strongly structure Betulaceae-associated ectomycorrhizal fungal communities in Chinese secondary forest ecosystems.
      ). Whether such a conclusion can also be applied to the microbial communities of spontaneously fermented milk requires further research.
      Whether at the provincial, county, or village level, spontaneous milk fermentation is a complex microecosystem, and many different factors can influence bacteria composition and further metabolic profiles. On different farms, farmers may employ different fermentation processes and different traditional fermentation conditions, such as temperature, pH, and time, which are known to influence bacterial growth and diversity (
      • Gatti M.
      • Bottari B.
      • Lazzi C.
      • Neviani E.
      • Mucchetti G.
      Microbial evolution in raw-milk, long-ripened cheeses produced using undefined natural whey starters.
      ). Similarly, milking machines can contain a reservoir of microorganisms, and thus, unsurprisingly, differences between machines and related practices can influence the microbial population of spontaneously fermented milk (
      • Verdier-Metz I.
      • Michel V.
      • Delbès C.
      • Montel M.C.
      Do milking practices influence the bacterial diversity of raw milk?.
      ).
      • Zhang R.
      • Huo W.
      • Zhu W.
      • Mao S.
      Characterization of bacterial community of raw milk from dairy cows during subacute ruminal acidosis challenge by high-throughput sequencing.
      described the effects of different dairy cattle diets (high-concentrate versus low-concentrate diets) on milk microbial communities, and the results indicated diet-associated differences in milk microbial communities. The bovine teat surface can contain a high diversity of bacteria, and it was also noted that the composition of the microbial community on the teat surface varied qualitatively and quantitatively from one farm to another (
      • Verdier-Metz I.
      • Gagne G.
      • Bornes S.
      • Monsallier F.
      • Veisseire P.
      • Delbès-Paus C.
      • Montel M.C.
      Cow teat skin, a potential source of diverse microbial populations for cheese production.
      ). The microbial community of cow teat skin mainly depends on the dairy farm environment, such as type of forage plants and housing conditions (
      • Zdanowicz M.
      • Shelford J.A.
      • Tucker C.B.
      • Weary D.M.
      • von Keyserlingk M.A.G.
      Bacterial populations on teat ends of dairy cows housed in free stalls and bedded with either sand or sawdust.
      ;
      • Verdier-Metz I.
      • Gagne G.
      • Bornes S.
      • Monsallier F.
      • Veisseire P.
      • Delbès-Paus C.
      • Montel M.C.
      Cow teat skin, a potential source of diverse microbial populations for cheese production.
      ). The cow's diet can modulate metabolites in raw milk (
      • Villeneuve M.P.
      • Lebeuf Y.
      • Gervais R.
      • Tremblay G.F.
      • Vuillemard J.C.
      • Fortin J.
      • Chouinard P.Y.
      Milk volatile organic compounds and fatty acid profile in cows fed timothy as hay, pasture, or silage.
      ), and differences in metabolic and nutritional composition of raw milk may also account for variation of bacterial communities (
      • Wang Z.
      • Jiang S.
      • Ma C.
      • Huo D.
      • Peng Q.
      • Shao Y.
      • Zhang J.
      Evaluation of the nutrition and function of cow and goat milk based on intestinal microbiota by metagenomic analysis.
      ), because effects of milk nutrients on bacterial composition of fermented milk are due to different responses of different bacterial species to variations in milk nutrient availability, such as carbohydrates and protein. However, the influences of these factors on bacterial composition and nonvolatile components remains unclear, and their specific contributions need to be calculated and elucidated in future research.
      In addition to their contribution to milk fermentation by transforming lactose into lactic acid, bacterial communities play important roles in the nonvolatile components and influence the organoleptic properties of fermented milk, contributing to a regional terroir (
      • Wouters J.T.M.
      • Ayad E.H.E.
      • Hugenholtz J.
      • Smit G.
      Microbes from raw milk for fermented dairy products.
      ). By identifying discriminatory nonvolatile components of different regions, it may be possible to define potentially predictive regional features. Furthermore, investigating underlying relationships between bacterial communities and discriminatory nonvolatile components contributes to understanding how bacterial communities influence nonvolatile components and what type of species influence the regional nonvolatile component characteristics of spontaneously fermented cow milk. In the present study, L. delbrueckii ssp. bulgaricus and Strep. thermophilus are predominant species; however, we found no correlations between discriminatory nonvolatile components and these 2 predominant species. Our results indicate that specific species with low abundance, rather than predominant species, might shape the regional metabolic characteristics of spontaneously fermented cow milk. Spontaneously fermented milk is a complex ecosystem, and species with high or low abundances make their own unique contributions to metabolic profiles. Thus, specific contributions of species with high or low abundances to regional nonvolatile component characteristics need further research. Meanwhile, such data give us some suggestions that we need to pay more attention to the selection and study of some nonpredominant bacterial species if we develop fermented milk products with regional characteristics. In addition, previous research has reported that growth and metabolic interactions occur between predominant species L. delbrueckii ssp. bulgaricus and Strep. thermophilus, which can result in higher acidification rates, stimulation of metabolite production, and increase of exopolysaccharide production (
      • Moon N.J.
      • Reinbold G.W.
      Commensalism and competition in mixed cultures of Lactobacillus bulgaricus and Streptococcus thermophilus..
      ;
      • Bouzar F.
      • Cerning J.
      • Desmazeaud M.
      Exopolysaccharide production and texture-promoting abilities of mixed-strain starter cultures in yogurt production.
      ;
      • Herve-Jimenez L.
      • Guillouard I.
      • Guedon E.
      • Boudebbouze S.
      • Hols P.
      • Monnet V.
      • Maguin E.
      • Rul F.
      Postgenomic analysis of Streptococcus thermophilus cocultivated in milk with Lactobacillus delbrueckii ssp. bulgaricus: Involvement of nitrogen, purine, and iron metabolism.
      ). In spontaneous milk fermentation, these interactions among microbial communities are more complex, which can directly influence the metabolic profiles of fermented milk products. Thus, further research is needed to determine whether interactions between predominant species and nonpredominant species influence regional metabolic characteristics of spontaneously fermented milk.
      Limosilactobacillus reuteri is a well-studied probiotic bacterium, and several possible beneficial effects on host health have been reported for specific strains, such as producing antimicrobial metabolites, anti-inflammatory effects, and inhibition of colonization of pathogenic microbes (
      • Mu Q.
      • Tavella V.J.
      • Luo X.M.
      Role of Lactobacillus reuteri in human health and diseases. 9:757.
      ). In the present study, L. reuteri was mainly positively correlated with amino acids and small peptide molecules, such as isoleucylphenylalanine, l-isoleucine, prolyltyrosine, and phenylalanylproline. However, it remains unclear whether these metabolites originated from raw cow milk or from fermentation, and this requires further study. Lactobacillus delbrueckii ssp. lactis was mainly positively correlated with indoleacrylic acid, indole-3-carboxaldehyde, 2-octenedioic acid, and indolelactic acid, and interestingly, these are all metabolites of tryptophan, which are synthesized by bacteria, particularly species of the genus Lactobacillus (
      • Zhang L.S.
      • Davies S.S.
      Microbial metabolism of dietary components to bioactive metabolites: Opportunities for new therapeutic interventions.
      ). This indicates that, in spontaneously fermented cow milk, L. delbrueckii ssp. lactis may synthesize these 4 metabolites, among which indole-3-carboxaldehyde has antifungal properties (
      • Brucker R.M.
      • Harris R.N.
      • Schwantes C.R.
      • Gallaher T.N.
      • Flaherty D.C.
      • Lam B.A.
      • Minbiole K.P.C.
      Amphibian chemical defense: Antifungal metabolites of the microsymbiont Janthinobacterium lividum on the salamander Plethodon cinereus.
      ). Limosilactobacillus fermentum is considered probiotic, or friendly bacteria, in humans and animals, and some strains have been applied to treat urogenital infections in women (
      • Reque E.D.F.
      • Pandey A.
      • Franco S.G.
      • Soccol C.R.
      Isolation, identification and physiological study of Lactobacillus fermentum LPB for use as probiotic in chickens.
      ;
      • Gardiner G.E.
      • Heinemann C.
      • Bruce A.W.
      • Beuerman D.
      • Reid G.
      Persistence of Lactobacillus fermentum RC-14 and Lactobacillus rhamnosus GR-1 but not L. rhamnosus GG in the human vagina as demonstrated by randomly amplified polymorphic DNA.
      ). Limosilactobacillus fermentum was mainly positively correlated with oryzalic acid B, which has antibacterial activity (
      • Kono Y.
      • Uzawa J.
      • Kobayashi K.
      • Suzuki Y.
      • Uramoto M.
      • Sakurai A.
      • Watanabe M.
      • Teraoka T.
      • Hosokawa D.
      • Watanabe M.
      • Kondo H.
      Structures of oryzalides A and B, and oryzalic acid A, a group of novel antimicrobial diterpenes, isolated from healthy leaves of a bacterial leaf blight-resistant cultivar of rice plant.
      ). The results showed that fermented cow milk from different sampling sites contains different potential probiotic LAB and bioactive metabolites.
      In the present study, we focused on the nonvolatile component characteristics and bacteria communities at the endpoint of spontaneously fermented cow milk rather than during fermentation. That is, these nonvolatile component characteristics and bacteria communities are the endpoints of the spontaneous fermentation and do not necessarily reflect those occurring during fermentation. Thus, it is also important to assess how bacterial composition evolves during milking, transport, storage, and dairy processing in different regions, and how different bacteria communities influence nonvolatile component profiles during fermentation. Such knowledge would provide valuable physiological insight into the dynamics of microbial growth and metabolic interactions in different regions, and this knowledge needs to be further investigated and supplemented.

      CONCLUSIONS

      In the present study, 25 sampling sites from Xin Jiang and Tibet were selected based on the distribution of spontaneously fermented cow milks and used for metagenomic and metabolomic analysis. The characteristics and relationships of bacterial communities and nonvolatile components at different scales of geographical distance (provincial, county, and village) were studied. At the provincial geographical level, we found significant differences (P = 0.008) in bacterial structure, shown in the PCoA plots, as well as significant differences [Simpson index (P = 0.012)] in bacterial diversity. However, at the county and village levels, we detected no significant differences (P > 0.05) in bacterial communities and diversity, but a difference in bacterial compositions was detectable. This indicates that bacterial communities and diversity of spontaneously fermented milk dissimilarity significantly increased with geographic distance. For the nonvolatile component profiles, the PLS-DA plot (R2Y > 0.5 and Q2 > 0.5) showed that samples from different geographical distances (provincial, county, and village levels) were all separated, indicating discriminations in nonvolatile component profiles. Investigating relationships between LAB and discriminatory nonvolatile components at the county level, we found that 9 species were positively correlated with 16 discriminatory nonvolatile components, and these were all species with low abundance, rather than predominant species.

      ACKNOWLEDGMENTS

      This study was funded by the China Agriculture Research System (Beijing, China; grant no. CARS-36) and Inner Mongolia Science and Technology Major Projects (Inner Mongolia; grant no. zdzx2018018). The authors declare no competing financial or other conflicts of interest.

      REFERENCES

        • Addis M.F.
        • Tanca A.
        • Uzzau S.
        • Oikonomou G.
        • Bicalho R.C.
        • Moroni P.
        The bovine milk microbiota: Insights and perspectives from -omics studies.
        Mol. Biosyst. 2016; 12 (27216801): 2359-2372
        • Bokulich N.A.
        • Thorngate J.H.
        • Richardson P.M.
        • Mills D.A.
        Microbial biogeography of wine grapes is conditioned by cultivar, vintage, and climate.
        Proc. Natl. Acad. Sci. USA. 2014; 111 (24277822): E139-E148
        • Bouzar F.
        • Cerning J.
        • Desmazeaud M.
        Exopolysaccharide production and texture-promoting abilities of mixed-strain starter cultures in yogurt production.
        J. Dairy Sci. 1997; 80: 2310-2317
        • Brucker R.M.
        • Harris R.N.
        • Schwantes C.R.
        • Gallaher T.N.
        • Flaherty D.C.
        • Lam B.A.
        • Minbiole K.P.C.
        Amphibian chemical defense: Antifungal metabolites of the microsymbiont Janthinobacterium lividum on the salamander Plethodon cinereus.
        J. Chem. Ecol. 2008; 34 (18949519): 1422-1429
        • Caporaso J.G.
        • Bittinger K.
        • Bushman F.D.
        • DeSantis T.Z.
        • Andersen G.L.
        • Knight R.
        PyNAST: A flexible tool for aligning sequences to a template alignment.
        Bioinformatics. 2010; 26 (19914921): 266-267
        • Cole J.R.
        • Chai B.
        • Farris R.J.
        • Wang Q.
        • Kulam-Syed-Mohideen A.S.
        • McGarrell D.M.
        • Bandela A.M.
        • Cardenas E.
        • Garrity G.M.
        • Tiedje J.M.
        The ribosomal database project (RDP-II): Introducing myRDP space and quality controlled public data.
        Nucleic Acids Res. 2007; 35 (17090583): D169-D172
        • Coorevits A.
        • De Jonghe V.
        • Vandroemme J.
        • Reekmans R.
        • Heyrman J.
        • Messens W.
        • De Vos P.
        • Heyndrickx M.
        Comparative analysis of the diversity of aerobic spore-forming bacteria in raw milk from organic and conventional dairy farms.
        Syst. Appl. Microbiol. 2008; 31 (18406093): 126-140
        • Corsetti A.
        • Settanni L.
        • Lopez C.C.
        • Felis G.E.
        • Mastrangelo M.
        • Suzzi G.
        A taxonomic survey of lactic acid bacteria isolates from wheat (Triticum durum) kernels and non-conventional flours.
        Syst. Appl. Microbiol. 2007; 30: 561-571
        • Doyle C.J.
        • Gleeson D.
        • O'Toole P.W.
        • Cotter P.D.
        Impacts of seasonal housing and teat preparation on raw milk microbiota: A high-throughput sequencing study.
        Appl. Environ. Microbiol. 2016; 83 (27815277): e02616-e02694
        • Edgar R.C.
        Search and clustering orders of magnitude faster than BLAST.
        Bioinformatics. 2010; 26 (20709691): 2460-2461
        • Gardiner G.E.
        • Heinemann C.
        • Bruce A.W.
        • Beuerman D.
        • Reid G.
        Persistence of Lactobacillus fermentum RC-14 and Lactobacillus rhamnosus GR-1 but not L. rhamnosus GG in the human vagina as demonstrated by randomly amplified polymorphic DNA.
        Clin. Vaccine Immunol. 2002; 9 (Clin. Diagn. Lab. Immunol.) (11777835): 92-96
        • Gatti M.
        • Bottari B.
        • Lazzi C.
        • Neviani E.
        • Mucchetti G.
        Microbial evolution in raw-milk, long-ripened cheeses produced using undefined natural whey starters.
        J. Dairy Sci. 2014; 97 (24290824): 573-591
        • Gesudu Q.
        • Zheng Y.
        • Xi X.X.
        • Hou Q.C.
        • Xu H.Y.
        • Huang W.Q.
        • Zhang H.P.
        • Menghe B.
        • Liu W.J.
        Investigating bacterial population structure and dynamics in traditional koumiss from Inner Mongolia using single molecule real-time sequencing.
        J. Dairy Sci. 2016; 99 (27522429): 7852-7863
        • Haas B.J.
        • Gevers D.
        • Earl A.M.
        • Feldgarden M.
        • Ward D.V.
        • Giannoukos G.
        • Ciulla D.
        • Tabbaa D.
        • Highlander S.K.
        • Sodergren E.
        Chimeric 16S rRNA sequence formation and detection in Sanger and 454-pyrosequenced PCR amplicons.
        Genome Res. 2011; 21: 494-503
        • Herve-Jimenez L.
        • Guillouard I.
        • Guedon E.
        • Boudebbouze S.
        • Hols P.
        • Monnet V.
        • Maguin E.
        • Rul F.
        Postgenomic analysis of Streptococcus thermophilus cocultivated in milk with Lactobacillus delbrueckii ssp. bulgaricus: Involvement of nitrogen, purine, and iron metabolism.
        Appl. Environ. Microbiol. 2009; 75 (19114510): 2062-2073
        • Hou Q.
        • Xu H.
        • Zheng Y.
        • Xi X.
        • Kwok L.
        • Sun Z.
        • Zhang H.
        • Zhang W.
        Evaluation of bacterial contamination in raw milk, ultra-high temperature milk and infant formula using single molecule, real-time sequencing technology.
        J. Dairy Sci. 2015; 98 (26476945): 8464-8472
        • Kono Y.
        • Uzawa J.
        • Kobayashi K.
        • Suzuki Y.
        • Uramoto M.
        • Sakurai A.
        • Watanabe M.
        • Teraoka T.
        • Hosokawa D.
        • Watanabe M.
        • Kondo H.
        Structures of oryzalides A and B, and oryzalic acid A, a group of novel antimicrobial diterpenes, isolated from healthy leaves of a bacterial leaf blight-resistant cultivar of rice plant.
        Agric. Biol. Chem. 1991; 55: 803-811
        • Liu W.
        • Zheng Y.
        • Kwok L.Y.
        • Sun Z.
        • Zhang J.
        • Guo Z.
        • Hou Q.
        • Menhe B.
        • Zhang H.
        High-throughput sequencing for the detection of the bacterial and fungal diversity in Mongolian naturally fermented cow's milk in Russia.
        BMC Microbiol. 2015; 15 (25887414): 45
        • Lozupone C.
        • Hamady M.
        • Knight R.
        UniFrac—An online tool for comparing microbial community diversity in a phylogenetic context.
        BMC Bioinformatics. 2006; 7 (16893466): 371
        • Mathara J.M.
        • Schillinger U.
        • Kutima P.M.
        • Mbugua S.K.
        • Holzapfel W.H.
        Isolation, identification and characterisation of the dominant microorganisms of kule naoto: The Maasai traditional fermented milk in Kenya.
        Int. J. Food Microbiol. 2004; 94 (15246238): 269-278
        • Moon N.J.
        • Reinbold G.W.
        Commensalism and competition in mixed cultures of Lactobacillus bulgaricus and Streptococcus thermophilus..
        J. Milk Food Technol. 1976; 39: 337-341
        • Mu Q.
        • Tavella V.J.
        • Luo X.M.
        Role of Lactobacillus reuteri in human health and diseases. 9:757.
        Front. Microbiol. 2018; 9 (29725324): 757
        • Parker M.
        • Zobrist S.
        • Donahue C.
        • Edick C.
        • Mansen K.
        • Nadjari M.H.Z.
        • Heerikhuisen M.
        • Sybesma W.
        • Molenaar D.
        • Diallo A.M.
        • Milani P.
        • Kort R.
        Naturally fermented milk from northern Senegal: Bacterial community composition and probiotic enrichment with lactobacillus rhamnosus..
        Front. Microbiol. 2018; 9 (30298060)2218
        • Peay K.G.
        • Garbelotto M.
        • Bruns T.D.
        Evidence of dispersal limitation in soil microorganisms: Isolation reduces species richness on mycorrhizal tree islands.
        Ecology. 2010; 91 (21302834): 3631-3640
        • Peay K.G.
        • Schubert M.G.
        • Nguyen N.H.
        • Bruns T.D.
        Measuring ectomycorrhizal fungal dispersal: Macroecological patterns driven by microscopic propagules.
        Mol. Ecol. 2012; 21 (22703050): 4122-4136
        • Perea Vélez M.
        • Hermans K.
        • Verhoeven T.L.A.
        • Lebeer S.E.
        • Vanderleyden J.
        • De Keersmaecker S.C.J.
        Identification and characterization of starter lactic acid bacteria and probiotics from Columbian dairy products.
        J. Appl. Microbiol. 2007; 103 (17714400): 666-674
        • Quigley L.
        • O'Sullivan O.
        • Beresford T.P.
        • Ross R.P.
        • Fitzgerald G.F.
        • Cotter P.D.
        High-throughput sequencing for detection of subpopulations of bacteria not previously associated with artisanal chesses.
        Appl. Environ. Microbiol. 2012; 78 (22685131): 5717-5723
        • Quigley L.
        • O'Sullivan O.
        • Stanton C.
        • Beresford T.P.
        • Ross R.P.
        • Fitzgerald G.F.
        • Cotter P.D.
        The complex microbiota of raw milk.
        FEMS Microbiol. Rev. 2013; 37 (23808865): 664-698
        • Redford A.J.
        • Bowers R.M.
        • Knight R.
        • Linhart Y.
        • Fierer N.
        The ecology of the phyllosphere: Geographic and phylogenetic variability in the distribution of bacteria on tree leaves.
        Environ. Microbiol. 2010; 12 (20545741): 2885-2893
        • Reque E.D.F.
        • Pandey A.
        • Franco S.G.
        • Soccol C.R.
        Isolation, identification and physiological study of Lactobacillus fermentum LPB for use as probiotic in chickens.
        Braz. J. Microbiol. 2000; 31: 303-307
        • Sun Z.
        • Liu W.
        • Bao Q.
        • Zhang J.
        • Hou Q.
        • Kwok L.
        • Sun T.
        • Zhang H.
        Investigation of bacterial and fungal diversity in tarag using high-throughput sequencing.
        J. Dairy Sci. 2014; 97 (25129502): 6085-6096
        • Turbes G.
        • Linscott T.D.
        • Tomasino E.
        • Waite-Cusic J.
        • Lim J.
        • Meunier-Goddik L.
        Evidence of terroir in milk sourcing and its influence on Cheddar cheese.
        J. Dairy Sci. 2016; 99 (27085404): 5093-5103
        • Vacheyrou M.
        • Normand A.C.
        • Guyot P.
        • Cassagne C.
        • Piarroux R.
        • Bouton Y.
        Cultivable microbial communities in raw cow milk and potential transfers from stables of sixteen French farms.
        Int. J. Food Microbiol. 2011; 146 (21429612): 253-262
        • Verdier-Metz I.
        • Gagne G.
        • Bornes S.
        • Monsallier F.
        • Veisseire P.
        • Delbès-Paus C.
        • Montel M.C.
        Cow teat skin, a potential source of diverse microbial populations for cheese production.
        Appl. Environ. Microbiol. 2012; 78 (22081572): 326-333
        • Verdier-Metz I.
        • Michel V.
        • Delbès C.
        • Montel M.C.
        Do milking practices influence the bacterial diversity of raw milk?.
        Food Microbiol. 2009; 26 (19269573): 305-310
        • Villeneuve M.P.
        • Lebeuf Y.
        • Gervais R.
        • Tremblay G.F.
        • Vuillemard J.C.
        • Fortin J.
        • Chouinard P.Y.
        Milk volatile organic compounds and fatty acid profile in cows fed timothy as hay, pasture, or silage.
        J. Dairy Sci. 2013; 96 (24035021): 7181-7194
        • Wang Y.L.
        • Gao C.
        • Chen L.
        • Ji N.N.
        • Wu B.W.
        • Li X.C.
        • Lü P.P.
        • Zheng Y.
        • Guo L.D.
        Host plant phylogeny and geographic distance strongly structure Betulaceae-associated ectomycorrhizal fungal communities in Chinese secondary forest ecosystems.
        FEMS Microbiol. Ecol. 2019; 95 (30889238)
        • Wang Z.
        • Jiang S.
        • Ma C.
        • Huo D.
        • Peng Q.
        • Shao Y.
        • Zhang J.
        Evaluation of the nutrition and function of cow and goat milk based on intestinal microbiota by metagenomic analysis.
        Food Funct. 2018; 9 (29577121): 2320-2327
        • Whipps J.M.
        • Hand P.
        • Pink D.
        • Bending G.D.
        Phyllosphere microbiology with special reference to diversity and plant genotype.
        J. Appl. Microbiol. 2008; 105 (19120625): 1744-1755
        • Wouters J.T.M.
        • Ayad E.H.E.
        • Hugenholtz J.
        • Smit G.
        Microbes from raw milk for fermented dairy products.
        Int. Dairy J. 2002; 12: 91-109
        • Xi X.M.
        • Kwok L.Y.
        • Wang Y.N.
        • Ma C.
        • Mi Z.H.
        • Zhang H.P.
        Ultra-performance liquid chromatography-quadrupole-time of flight mass spectrometry MSE-based untargeted milk metabolomics in dairy cows with subclinical or clinical mastitis.
        J. Dairy Sci. 2017; 100 (28342601): 4884-4896
        • Yu J.
        • Mo L.
        • Pan L.
        • Yao C.
        • Ren D.
        • An X.
        • Tsogtgerel T.
        • Zhang H.
        • Liu W.
        Bacterial microbiota and metabolic character of traditional sour cream and butter in Buryatia, Russia.
        Front. Microbiol. 2018; 9 (30459729)2496
        • Yu J.
        • Wang W.H.
        • Menghe B.L.
        • Jiri M.T.
        • Wang H.M.
        • Liu W.J.
        • Bao Q.H.
        • Lu Q.
        • Zhang J.C.
        • Wang F.
        • Xu H.Y.
        • Sun T.S.
        • Zhang H.P.
        Diversity of lactic acid bacteria associated with traditional fermented dairy products in Mongolia.
        J. Dairy Sci. 2011; 94 (21700007): 3229-3241
        • Zdanowicz M.
        • Shelford J.A.
        • Tucker C.B.
        • Weary D.M.
        • von Keyserlingk M.A.G.
        Bacterial populations on teat ends of dairy cows housed in free stalls and bedded with either sand or sawdust.
        J. Dairy Sci. 2004; 87 (15453481): 1694-1701
        • Zhang L.S.
        • Davies S.S.
        Microbial metabolism of dietary components to bioactive metabolites: Opportunities for new therapeutic interventions.
        Genome Med. 2016; 8 (27102537): 46
        • Zhang R.
        • Huo W.
        • Zhu W.
        • Mao S.
        Characterization of bacterial community of raw milk from dairy cows during subacute ruminal acidosis challenge by high-throughput sequencing.
        J. Sci. Food Agric. 2015; 95 (24961605): 1072-1079