Procedures for collection of samples were approved by the Animal Welfare Committee of China Agricultural University, Beijing, China and were conducted strictly following the regulations and guidelines established by this committee (approval number: AW10111202-1-3). All cows sampled in this research were guaranteed to have no records of clinical mastitis for three consecutive months, and were reared in the same dairy farm in Beijing, Northern China. Milk samples collection was performed after a regular disinfection procedure i.e. teat cleaning with 70% ethanol and the first three streams of milk were discarded [8]. Milk from each of the four quarters was collected in a 50 mL-sterile falcon tube, evenly mixed together and kept on dry ice during the whole collection process. For each cow, two milk samples were taken with one put in −80°C for later DNA extraction and the other sent to Beijing Dairy Cattle Center for SCC measurement. Other phenotype information was acquired according to records provided by the dairy farm. A total of 23 peripheral blood samples were aseptically taken from caudal vein and stored in 10-mL pre-added ethylenediaminetetra-acetic acid-K2 vacutainers on the same day the milk samples were collected.

Total genomic DNA was extracted from milk using cetyltrimethyl ammonium bromide. The integrity and purity of DNA was tested using 1% agarose gel electrophoresis. Total RNA was isolated from peripheral white blood cells using Trizol reagent (Ambion, Waltham, MA, USA) following the manufacturer’s instructions. Purity of RNA was measured by Nanodrop (Thermo Fish, Waltham, MA, USA) and RNA integrity was assessed using electrophoresis and Aglient 2100 BioAnalyzer. Qubit 2.0 Fluorimeter was used for accurate quantification of DNA and RNA. Both metagenomic shotgun sequencing and RNA-seq were performed on Illumina NovaSeq 6000 platform (Illumina, San Diego, CA, USA), preceded by a construction of 150-bp paired-end library. The average percent of bases scoring Q30 or above were over 90% for reads generated from either metagenomic shotgun sequencing or RNA-seq techniques.

Based on a threshold of 100,000 SCC/mL, the samples were divided into MS (medium SCC) and LS (low SCC) group [8]. Low quality reads, adaptors and host reads filtering were accomplished using KneadData. The quality of all clean data was assessed utilizing FastQC first before subsequent analyses. Non-host reads that passed quality control were taken for taxonomic classification using Kraken2 and abundance re-estimation using Bracken (v 2.6). Eukaryotes and viruses were discarded from the analysis. De novo assembly of reads into contigs was conducted by Megahit (v 1.2.8) and protein-coding gene prediction was performed using Prodigal (v 2.6.3). Redundant sequences of coding sequence predicted were then removed by CD-HIT (v 4.8.1) in order to facilitate downstream analysis. The remaining non-redundant representative sequences then served as an index for Salmon alignment-free quantification and also the input to eggnog-mapper (v 2.0.1) for functional annotation. Alpha diversity indices (Shannon, Chao1, Simpson) at different sampling depths were calculated using function diversity using R package vegan. Difference of alpha diversity between the two groups was examined using Wilcoxon rank-sum test. Bray-Curtis distance matrices based on relative abundance table at different taxonomic levels were built up using function vegdist in vegan and permutational multivariate analysis of variance was done by adonis to examine the effect of various experimental factors (SCC, parity, lactation stage) on differences among samples (Bray-Curtis distance ~ SCC + parity + lactation stage). Dissimilarity among samples (beta-diversity) was visualized with function pcoa in ape package. To eliminate the effect of library size, edgeR was utilized to determine differentially abundant taxa (DAT) and metagenomic function.

Raw reads generated from RNA-seq were subjected to strict quality control (Trimmomatic v0.39) and were then aligned to Bos taurus ARS-UCD1.2 using STAR (v 2.7.9a) with an average mapping rate over 90%. Reads mapped to the reference genome were quantified by feature Counts in subread (v 2.0.2) package followed by TPM normalization using StringTie (v 2.1.5). Subsequent analyses were carried out using R package. Differentially expressed genes (DEG) were identified by edgeR. Gene ontology (GO) and Kyoto encyclopedia of genes and genomes (KEGG) pathway enrichment analysis were conducted using online website KOBAS (kobas.cbi.pku.edu.cn/kobas3).

To find the correlation between co-expression gene modules and DAT, weighted correlation network analysis and module-trait relationship assessment were implemented by WGCNA using top 5,000 genes with the largest median absolute deviation (MAD). Pearson correlation coefficients of the relative abundance of DAT and genes both identified as DEGs and found in modules associated with DAT were calculated using function cor.test in R.

A total of 31 (MS = 13, LS = 18) milk samples were obtained and used for sole metagenomics analysis. The parity of each selected cow was less than three and days in lactation ranged from 54 to 164 d (97±30.5 d). Distribution of SCC, parity and lactation are presented in Supplementary Figure S1. In general, 23 (MS = 8, LS = 15) samples that had both transcriptome data and metagenome data were contained in integrative analysis.

A total of 35 phyla were observed, of which 4 phyla (Proteobacteria 54.73±14.35; Actinobacteria 16.47±6.56; Firmicutes 19.71±12.83; and Bacteroidetes 9.40±6.70) were shared across all samples. Top 6 most abundant phyla selected based on the average relative abundance were Proteobacteria, Actinobacteria, Firmicutes, Bacteroidetes, Tenericutes, and Spirochaetes (Figure 1). Additionally, 70 classes, 166 orders, 375 families and 1,413 genera were identified. Species level were excluded from subsequent analyses due to inadequate sequencing depth for accurate classification at this level.

To compare richness of the microbial community in the two groups, three different indices were used to measure the overall diversity in each sample (Figure 2). The significance was calculated using paired Wilcoxon rank-sum test. The results showed that samples, subsampled at a depth of 12,000, 24,000, and 30,000 had no significant differences between MS and LS group using any of the three alpha diversity indices, indicating the richness difference was quite subtle. We further explored the impact of SCC, parity and lactation stage on dissimilarities of microbial composition among samples using PERMANOVA. The result suggested that SCC has little impact on Bray-Curtis dissimilarity at almost all classification levels (p>0.05). Moreover, parity and lactation stage also accounted for little of the variation on overall microbial composition (Table 1). Principal coordinate analysis also showed no clustering tendency within either group and that the discrepancy between MS and LS samples remained unclear at any classification level (Supplementary Figure S2).

To evaluate the difference of relative abundance of specific taxon between the two groups, we utilized edgeR to normalize the counts and then identify SCC-associated taxa. Genus level was mainly focused because the difference in other lower-ranked classification levels was already captured at genus level. To further avoid the bias caused by sequencing depth, only taxa existent in at least 30% of samples in either group were kept for differential abundance analysis. Differentially abundant taxa were selected based on |log2 (fold change)| ≥1 and p<0.05 (Figure 3). A total of 7 genera were identified using this threshold. Two Gram-negative bacteria - Ralsotina and Sphingomonas were enriched in MS group while the abundance of Xanthomonas, Psychrobacter, Moraxella, Pantoea, Staphylococcus was significantly higher in LS group.

In order to find out whether certain host genes could potentially get involved in the interplay between bacteria and host, a combined analysis of DEG and weighted gene co-expression network was carried out to evaluate the correlation between host phenotype (SCC), gene expression pattern and the relative abundance of DAT. We first performed a DEGs analysis in MS and LS group. DEGs were identified with |log2 (fold change)| ≥1 and p<0.05. In general, 355 genes were differentially expressed in the two groups, of which 297 genes were up-regulated and 58 were down-regulated in MS compared with LS. Next we ranked genes according to MAD and the top 5000 genes were selected to be clustered into modules by WGCNA based on the correlation.

To further explore which module was potentially related to DAT, relative abundance of each seven differentially abundant genera identified above was used as a trait and the correlation between it and each module eigengene was evaluated. Four (M10, M11, M13, M14) out of 14 modules exhibited significant relevance to the relative abundance of Ralsotina and Sphingomonas, which were enriched in MS group (Figure 4a). A total of 297 up-regulated DEGs were mapped to genes in these four modules to find genes that were both related to SCC and MS-enriched bacteria. 0, 21, 2, 2 genes in M10, M11, M13, M14 overlapped with those in DEGs, respectively (Supplementary Figure S3). Relatively high overlap between genes in M11 and DEGs provided evidence that genes in M11 might play a role in host-microbiota interaction.

To get functional interpretation of these genes and pathways they were involved in, we conducted GO and KEGG pathways enrichment analysis. Because of the inadequate number of overlapping genes to complete enrichment analyses, an overlap between GO and KEGG pathways enriched by genes in M11 and DEGs were investigated instead. As shown in Figure 4b, approximately nine KEGG pathways were shared by both sets. Several pathways related to inflammation and immune response were identified, such as phagosome, cytokine-cytokine receptor interaction, interleukin-17 (IL-17) signaling pathway, chemokine signaling pathway, apoptosis. The overlapping GO terms that fell in biological process categorical included chemokine-mediated signaling pathway, antimicrobial humoral immune response, neutrophil chemotaxis, which were also associated with immune signals, suggesting that there was an increase in the immune system activities in MS group compared with LS group and were consistent with the physiological process during subclinical mastitis, because subclinical mastitis was usually characterized by active immune responses (Supplementary Figure S4). Equally important is that the result also demonstrated a higher level of kinase and phosphatases activities in MS which was typically involved in initiation of immune cells activation in response to antigens. Despite the fact that none of the 21 genes shared by DEGs and M11 was enriched in these pathways or terms, the result infers that the intensity of immune responses varied from subclinical mastitis cows to healthy cows, and also from cows with higher abundance of Ralstonia and Sphingomonas to cows with lower abundance of those, though the variation might not be obvious. However, one gene, F11R was found to be involved in two different enriched KEGG pathways. F11R was significantly enriched in tight junction using DEGs for pathway enrichment analysis and leukocyte transendothelial migration using M11 genes, indicating its role in initiation of inflammatory immune response and barrier function maintenance. To further determine a relevance between the change in expression level of some specific genes and the relative abundance of Ralstonia and Sphingomonas, Pearson correlation was used to find genes directly linked to these two kinds of bacteria (Figure 4d). The results of correlation test showed that transcobalamin 1 (TCN1) and uridine phosphorylase 1 (UPP1) had a strong correlation with either Ralstonia (0.78, 0.69, respectively) or Sphingomonas (0.85, 0.76, respectively).

Based on the functional annotation of metagenome, we sought to detect the biological processes that either intensified or weakened in MS group. Differentially enriched GO terms determined by edgeR are shown in Figure 5, of which four were less abundant in MS compared with LS, and 18 showed the opposite tendency. The four down-regulated GO in MS were organic hydroxy compound metabolic process, transmembrane transporter activity, kinase activity and response to organic substance. Notably, lipopolysaccharide (LPS) biosynthetic and metabolic process, response to oxygen levels and hypoxia were up-regulated. It is noteworthy that the activities of some enzymes involved in DNA replication and nucleotide synthesis also increased in MS, suggesting active bacterial reproduction was underway.