Effect of Production Environment on Microbial Diversity in Raw Milk from Hami and Altay Regions Using Single-Molecule Real-Time Sequencing Technologye

1. Introduction

Raw milk harbors a complex and dynamic microbiome, thriving in a nutrient-rich aqueous environment that supports rapid microbial proliferation. These microbial consortia encompass spoilage agents, pathogens, beneficial probiotics, and technologically relevant starter cultures (Wang et al., 2019; Dufour et al., 2019; Moineau-Jean et al., 2019). Following milking and during refrigerated storage, microbial dynamics persist, ultimately influencing the microbiological quality of milk upon arrival at processing facilities (Castilho et al., 2018; Isele et al., 2019; Srisengfa et al., 2016). The raw milk microbiota is typically characterized by Gram-positive taxa, including Bacillus, Micrococcus, and Staphylococcus, alongside various lactic acid bacteria (LAB), such as Lactococcus and Streptococcus. Conversely, Gram-negative psychrotrophs, particularly Pseudomonas, Aeromonas, and Enterobacteriaceae members (e.g., Enterobacter, Hafnia, Klebsiella), pose significant spoilage risks due to their resilience. Yeasts, including Candida and Kazachstania, also contribute to the fungal diversity of raw milk (Krewinkel et al., 2015).

Psychrotrophic bacteria, such as Pseudomonas spp., Aeromonas spp., and Bacillus spp., dominate the cold-chain ecology (Qin et al., 2014). Of particular concern are the thermostable enzymes, specifically lipases and proteases, secreted by these microorganisms. These enzymes can withstand pasteurization and ultra-high temperature (UHT) treatments, representing the primary agents of spoilage that compromise the sensory and functional quality of dairy products. Specifically, lipases catalyze triglyceride hydrolysis, resulting in lipolytic rancidity and butyric acid off-flavors. Concurrently, proteases hydrolyze casein fractions, generating bitter peptides and inducing age-related gelation, thereby reducing shelf life (Xiang et al., 2024).

The taxonomic composition of raw milk is primarily dictated by direct-contact microbial reservoirs, notably the teat apex, udder skin, and milking apparatus. Secondary influences stem from indirect-contact environments, including bedding materials, feces, feed, water sources, and the broader barn and parlor environments (Stoeckel et al., 2016). Subsequent contamination can occur during milking, transportation, and storage, posing substantial risks to both the safety and processing quality of the final product. Consequently, characterizing the raw milk microbiome is critical for ensuring product safety and optimizing downstream processing protocols.

Traditional culture-independent methods, such as Single-Molecule Real-Time (SMRT) sequencing, bypass the limitations of microbial cultivation by analyzing total DNA extracted directly from samples. This technology offers superior resolution and throughput compared to conventional culture-dependent microbiological techniques (Fenelon et al., 2020). By applying SMRT sequencing to metagenomic DNA, researchers can achieve precise taxonomic classification down to the species level, significantly accelerating the diagnostic timeline (Hou et al., 2016). Furthermore, third-generation sequencing facilitates a holistic interrogation of microbial genomic architectures within environmental niches, thereby maximizing the exploitable potential of bioresources (Zhou et al., 2024). Specifically, targeting the full-length 16S rRNA gene allows for high-fidelity resolution of bacterial taxa, enabling localization to the species level (Liang et al., 2022).

Herein, we conducted a comprehensive investigation into the microbial ecology of raw camel milk and its associated production environments (i.e., feces, bedding, feed, and teat skin) to elucidate the overall microbial diversity. The primary objectives were to guide post-harvest processing strategies, optimize disease prevention protocols, and identify potential beneficial microbial resources for industrial exploitation. This work establishes a theoretical framework for understanding the dynamics of microbial contamination within camel milk production systems in Xinjiang, China. By mitigating potential contamination risks and deepening insights into raw milk microbiology, this study furnishes detailed contextual data to support dairy processing, enhance safety assurance, and facilitate novel strain discovery. Ultimately, these findings lay a robust scientific foundation for the systematic monitoring of raw milk quality and safety.

2. Materials and Methods

2.1. Sample Collection

Sampling was conducted across camel breeding farms in the Altay and Hami regions of Xinjiang, China. Environmental matrices were systematically collected, including bedding and feed from pens, feces and soil from the vicinity of the farms, and teat skin swabs from lactating camels. For each farm and sample type, three independent biological replicates were collected. Following quality control, a total of 30 samples were retained for downstream analysis. These comprised five categories (cushion, feces, teat skin, feed, and soil) from both regions, with three replicates per category. The coding scheme denoted samples from Altay as "07" and Hami as "09" , with suffixes -1 to -3 indicating replicates (See Supplementary Table S1 for detailed sample IDs).

2.2 DNA Extraction

Genomic DNA was extracted from all samples using the TIANamp Bacterial DNA Kit (Tiangen Biotech, Beijing, China), following the manufacturer’s protocol. The concentration and purity (A260/A280 ratio) of the extracted DNA were quantified using a NanoPhotometer® N60 (Implen GmbH, Munich, Germany). Subsequently, DNA integrity and fragment size were verified by electrophoresis on a 1% (w/v) agarose gel. Only samples with high purity (OD260/280 ratio between 1.8 and 2.0) and intact high-molecular-weight bands were selected for subsequent library construction.

2.3 PCR Amplification

Specific barcoded primers were synthesized to target the hypervariable V4 region of the 16S rRNA gene. To ensure data reproducibility and minimize PCR bias, amplification was performed using the minimum number of cycles required to obtain sufficient product, with strict consistency in cycle numbers across all samples. A preliminary trial was conducted on representative samples to determine the optimal cycle threshold for adequate yield.

PCR amplification was carried out using a ProFlex PCR System (GeneAmp PCR System 9700, Applied Biosystems, Foster City, CA, USA) with TransStart Fastpfu DNA Polymerase (TransGen Biotech, Beijing, China). Each sample was amplified in triplicate under standard cycling conditions. Replicate amplicons for each sample were pooled and visualized via 2% agarose gel electrophoresis. Target DNA fragments were excised from the gel and purified using the AxyPrep DNA Gel Recovery Kit (Axygen Biosciences, Union City, CA, USA). The purified products were eluted with Tris-HCl buffer (pH 8.8) and subsequently verified by 2% agarose gel electrophoresis.

Library quantification was performed using the QuantiFluor™ ST Blue Fluorescence System (Promega, Madison, WI, USA), ensuring a robust two-tiered quality control (QC) strategy prior to SMRT sequencing. Finally, libraries were normalized and pooled in equimolar ratios according to the sequencing platform's requirements.

2.4. Fluorescence Quantification

Library concentrations were accurately determined using the QuantiFluor™ ST Blue Fluorescence System (Promega, Madison, WI, USA). Subsequently, libraries were normalized and pooled in equimolar ratios based on the sequencing platform's requirements.

2.5. Library Construction and Sequencing

SMRTbell libraries were constructed using the SMRTbell Express Template Prep Kit 2.0 (Pacific Biosciences, Menlo Park, CA, USA) following the manufacturer’s instructions. Briefly, the PCR amplicons were purified using AMPure PB beads (Pacific Biosciences) to remove short fragments and primers. Size selection was performed via 2% agarose gel electrophoresis to ensure the integrity of the inserts. The purified DNA was eluted in TE buffer and quantified using a Qubit dsDNA HS Assay Kit (Invitrogen, Carlsbad, CA, USA). Finally, the double-stranded DNA templates were denatured and bound to the polymerase to form sequencing-ready complexes for Single-Molecule Real-Time (SMRT) sequencing.

2.6. SMRT Sequencing

Sequencing was performed on a PacBio Sequel II system (Pacific Biosciences, Menlo Park, CA, USA) following the manufacturer’s guidelines. The sequencing process relies on Single-Molecule Real-Time (SMRT) technology, wherein DNA polymerase complexes are immobilized at the bottom of zero-mode waveguides (ZMWs). Sequencing primers and the SMRTbell templates were bound to form sequencing-ready complexes. Upon loading into the SMRT Cell, the diffusion of nucleotides allowed for real-time monitoring of fluorescence pulses as phospholinked nucleotides were incorporated. This setup enabled the generation of high-fidelity reads, facilitating comprehensive profiling of the microbial communities.

2.7. Bioinformatics Analysis

Raw subreads generated from SMRT sequencing were processed using the ccs (Circular Consensus Sequencing) algorithm to generate High-Fidelity (HiFi) reads. Following rigorous quality filtering and primer trimming, demultiplexed reads were denoised and clustered into Amplicon Sequence Variants (ASVs) or Operational Taxonomic Units (OTUs) using the QIIME2 pipeline (version 2023.2). Taxonomic assignment was performed by querying the sequences against the SILVA reference database (release 138).

Alpha diversity metrics (including Shannon, Simpson, and Chao1 indices) and beta diversity analyses (such as Principal Coordinates Analysis, PCoA) were calculated based on the ASV/OTU table to evaluate microbial richness and structural variations. Taxonomic composition was statistically analyzed across all hierarchical levels (from phylum to genus) to characterize the microbial landscape. Furthermore, multivariate statistical analyses (e.g., ANOSIM, PERMANOVA) and differential abundance testing (e.g., DESeq2, LEfSe) were conducted to identify significant variations in microbial composition between sample groups.

3. Results

3.1. Overall Microbiome Analysis

Taxonomic profiling revealed a total of 35 bacterial phyla and 100 genera across all samples. The relative abundance of the dominant phyla is illustrated in Figure 1A, while the genus-level distribution is presented in Figure 1B (taxa with an abundance < 1.0% or those remaining unclassified were aggregated as "Others").

The majority of OTUs were assigned to three predominant phyla: Firmicutes (43.20%), Proteobacteria (20.12%), and Actinobacteriota (12.02%). In the Altay region, Firmicutes was the predominant phylum (43.04%). At the genus level, Rothia (20.63%) and sequences identified as Chloroplast (18.54%) represented the most abundant taxa. Similarly, Firmicutes dominated the Hami region (46.80%), although a high abundance of Chloroplast (18.23%) and Mitochondria (9.22%) sequences was also observed (Fig. 1A, B). It is noteworthy that the high abundance of chloroplast and mitochondrial sequences likely reflects the presence of plant-derived material and eukaryotic cells in the samples, which is a common occurrence in 16S rRNA amplicon sequencing of environmental samples.

3.2. Microbiota Diversity Based on Alpha and Beta Analysis

Alpha diversity was evaluated using the Chao1 and ACE indices (richness), alongside the Shannon and Simpson indices (diversity and evenness). As illustrated in Figure 2, the ACE index revealed a significantly higher community richness in soil samples from Altay compared to those from Hami (P < 0.05). In contrast, no significant differences were observed in Shannon or Simpson indices across most sample types between the two regions (P > 0.05).

Specifically, fecal samples exhibited significantly higher Chao1 and ACE values in the Altay region than in the Hami region (P < 0.05, Table 1). Similarly, feed samples showed significantly greater richness in Altay based on the ACE index, suggesting regional disparities in the abundance of microbial species.

To visualize structural variations, Principal Coordinates Analysis (PCoA) based on Bray-Curtis dissimilarity was employed. The PCoA plot (Fig. 2B) revealed distinct clustering patterns between the Altay and Hami groups. The ANOSIM test confirmed that the inter-group differences were significantly larger than the intra-group variations (R = 0.7503, P = 0.001), indicating distinct microbial profiles between the two regions. Furthermore, the NMDS ordination yielded a low stress value of 0.0094 (Fig. 2C), indicating an excellent fit of the data in the reduced dimensional space.

The heatmap visualization of genus-level abundance (Fig. 3) revealed distinct clustering patterns between the Altay and Hami regions, reflecting significant differences in microbial composition. Achromobacter was consistently detected across nipple skin samples from both regions. In contrast, sequences identified as mitochondria (likely representing eukaryotic contaminants or degraded host DNA) were exclusively found in feed samples from Hami.

Anaerobic taxa, including Clostridium sensu stricto, Romboutsia, and Paeniclostridium, exhibited a broad distribution across soil and fecal samples in both regions, with an extended presence in feed and bedding samples specifically within the Hami region. Conversely, specific genera, such as Bifidobacterium, Treponema, and Bacteroides, were predominantly restricted to soil, bedding, and fecal samples in the Hami region, with sporadic detection in Altay feces.

Certain taxa demonstrated regional specificity; for instance, Bacillus was exclusively identified in fecal, soil, and bedding samples from the Altay region. Similarly, Salinimicrobium, Flavobacterium, and Halomonas were unique to the soil and bedding of Altay. In Hami, Christensenellaceae R-7 group and Rikenellaceae RC9 group were prevalent in soil and bedding. The genus Jeotgalicoccus was primarily associated with nipple skin and bedding in Hami, whereas Salinicoccus showed a ubiquitous presence across bedding and skin in both regions. These distinct distribution profiles underscore the feasibility of differentiating bacterial communities between the two regions based on environmental reservoirs.

Figure 3. Heat map of clustering of bacterial genera in all samples from both regions.

Figure 3. Heat map of clustering of bacterial genera in all samples from both regions.

3.3. Common and Difference Analysis of Bacterial Compositions

Venn diagram analysis revealed a core microbiome comprising 21 shared genera and 23 shared species across both regions, representing 84.02% and 82.95% of the total taxonomic units identified, respectively (Figure 4A).

Fecal samples were dominated by Romboutsia (10.77%), followed by Paeniclostridium (8.10%) and Jeotgalicoccus (4.85%) (Figure 4B). Nipple skin samples exhibited a distinct profile, with Achromobacter spp. (22.23%) and Moraxella (6.29%) being predominant, alongside Clostridium sensu stricto (5.00%) (Figure 4C). Feed samples were characterized by a high abundance of chloroplast sequences (54.47%) and mitochondrial sequences (17.08%), likely reflecting plant-derived material and eukaryotic DNA, along with members of the Planococcaceae family (3.63%) (Figure 4D). Soil samples were primarily composed of Planococcaceae (7.05%), Clostridium sensu stricto (4.16%), and Romboutsia ilealis (2.95%) (Figure 4E). Finally, bedding samples were dominated by Corynebacterium maris (5.98%), Jeotgalicoccus spp. (5.61%), and Salinicoccus spp. (4.62%) (Figure 4F).

To identify the specific bacterial taxa driving the significant differences in community composition between the Altay and Hami regions, biomarker analysis was performed using the LEfSe method. Taxa labeled as "uncultured," "unranked," or "unclassified" were excluded from the dataset prior to analysis. Biomarker discovery was then conducted using Linear Discriminant Analysis (LDA) Effect Size (LEfSe), which effectively discriminates between biologically relevant groups by identifying features with significant effect sizes.

As illustrated in Figure 5A–C, Firmicutes was the dominant phylum in fecal samples from Altay, whereas Actinobacteriota predominated in those from Hami. A similar pattern was observed in teat skin samples (Firmicutes in Altay vs. Actinobacteriota in Hami). In soil samples, Actinobacteriota dominated in Altay, while Firmicutes was predominant in Hami. Feed samples lacked a dominant phylum in Altay but were dominated by Firmicutes in Hami. Additionally, bedding samples in Altay were characterized by Actinobacteriota, whereas those in Hami were dominated by Firmicutes.

These findings are consistent with the overall beta-diversity trends observed in previous sections. Furthermore, LEfSe identified 10 bacterial phyla with LDA scores exceeding 4.0, highlighting distinct taxonomic signatures between the regions.

3.4. Model Prediction and Predictive Functions

Functional profiling based on the Clusters of Orthologous Groups (COG) database revealed that the majority of genes were annotated to metabolic functions, including global and overview maps, as well as secondary metabolite biosynthesis (Fig. 6A). Predominant pathways specifically encompassed amino acid transport and metabolism, alongside the transport and metabolism of inorganic ions. Furthermore, secondary metabolite biosynthesis was linked to energy production and conversion, as well as translation, ribosomal structure, and biogenesis.

These results highlight the interconnected nature of the cellular metabolic network, emphasizing systemic functionality over isolated biochemical reactions (Fig. 6B). Comparative analysis indicated highly consistent COG functional profiles across both regions, suggesting conserved core metabolic potentials despite variations in taxonomic composition

Figure 6. Predictive functional analysis of microbial communities based on the COG database. (A) Similar COG functional classification profiles across all samples; (B) Key metabolic pathways associated with COG-encoded proteins and the structure of the overall cellular metabolic network.

Figure 6. Predictive functional analysis of microbial communities based on the COG database. (A) Similar COG functional classification profiles across all samples; (B) Key metabolic pathways associated with COG-encoded proteins and the structure of the overall cellular metabolic network.

4. Discussion

This study provides a comprehensive characterization of the microbial ecology across the raw camel milk production chain in the Altay and Hami regions of Xinjiang. The dominant phyla identified—Firmicutes (43.04%), Proteobacteria (21.98%), and Actinobacteriota (16.40%)—differ significantly from those reported in pasteurized milk studies. For instance, Ding et al. (2020) observed that Proteobacteria and Firmicutes each exceeded 50% in HTST-treated milk, with Pseudomonas dominating the genera. The lower relative abundance of these taxa in our raw samples suggests that pasteurization selectively eliminates sensitive strains, thereby enriching thermotolerant genera.

When compared to other non-bovine milks, our results align with the trends reported by Luoyizha et al. (2020) for donkey milk, particularly regarding the abundance of Bacteroidota (~10%). However, unlike Zhang et al. (2023), who found Firmicutes to constitute 96% of mare's milk, our camel milk samples exhibited a more balanced distribution among the top phyla. This variation likely stems from species-specific factors, regional climates, and dietary differences. Similarly, Ouyang et al. (2022) reported Proteobacteria dominance (43.3%–84.0%) in seawater, noting Cyanobacteria (0.3%–10.4%); our observed Cyanobacteria (7.91%) falls within this range, although in the context of ruminant feces/feed, these signals likely represent plant-derived chloroplast sequences rather than free-living cyanobacteria (Sánchez-Baracaldo, 2025).

Environmental sourcing analysis revealed distinct distribution patterns. In Altay, Firmicutes was highly prevalent in teat skin (49.58%) and soil (48.65%), suggesting these as primary reservoirs. Conversely, in Hami, fecal samples showed high abundances of Cyanobacteria (68.34%) and Firmicutes (29.07%). These findings corroborate Dong et al. (2021), who identified feces as a critical vector for Firmicutes transmission to milking equipment and storage tanks. Furthermore, Zhang et al. (2023) highlighted that Proteobacteria often originate from raw milk and water sources, which aligns with our observation of high Proteobacteria in Hami feed (21.80%) and soil (20.46%). The high prevalence of Cyanobacteria (often indicative of chloroplast sequences from plant matter) in feces and feed highlights the influence of roughage composition on the microbial landscape.

Functional profiling based on the COG database indicated that the microbiome is heavily geared towards metabolic activities, including amino acid transport and energy conversion. The consistent functional profiles across both regions, despite taxonomic variations, suggest a degree of functional redundancy within these ecosystems. Notably, the enrichment of Firmicutes and Bacteroidota is relevant to the production of short-chain fatty acids (SCFAs) like butyrate and propionate, which can suppress pro-inflammatory cytokines and strengthen epithelial barriers (Zhi et al., 2022). However, the presence of psychrotrophic taxa implies risks associated with heat-stable lipases and proteases (Xiang et al., 2024), which can lead to lipolytic rancidity and age-related gelation even after pasteurization (Fanelli et al., 2020).

Ensuring raw milk safety is paramount for public health and industry sustainability. Microbial contamination not only increases spoilage risks—often driven by Pseudomonas spp. and other psychrotrophs—but also undermines consumer confidence (Palii et al., 2020; Martin et al., 2021). Our study identifies critical control points, particularly the management of fecal contamination and bedding hygiene, supporting the notion that centrifugation and heat treatment can reduce microbial loads prior to processing (Peruzi et al., 2019; Jia et al., 2021). By mapping the environmental sources of key taxa, this work provides actionable evidence for upgrading quality surveillance systems and optimizing production practices in Xinjiang's dairy sector (Pierezan et al., 2022). Ultimately, the integration of SMRT sequencing data with traditional microbiological parameters offers a robust scientific foundation for regulatory bodies to implement precision monitoring across the entire dairy supply chain (Yuan et al., 2022).

5. Conclusions

This study systematically mapped the microbial landscapes of raw camel milk production chains in Xinjiang's Altay and Hami regions, resolving key taxa to the species level during the same seasonal window. We elucidated the primary environmental reservoirs of these microbes and assessed their implications for raw milk safety and processing stability, proposing targeted intervention strategies.

Furthermore, by characterizing bacterial diversity across critical nodes of the production continuum, this work deepens our understanding of how environmental factors shape the raw milk microbiome. These findings underscore the dual value of this research: advancing dairy safety protocols while identifying potential probiotic candidates or starter cultures. Ultimately, this integrated analysis provides a robust scientific foundation for regulatory bodies to implement precision monitoring and optimize quality assurance systems across the entire dairy supply chain.