All experimental procedures involved in this study were approved by the Institutional Animal Care and Use Committee of Sichuan Agricultural University (DKY-B20171906). The current study is a continuation of previous research, where the differentially expressed genes and pathways associated with ketosis were investigated using these animals (Wu et al., 2020). The animals and experimental design were fully described in the original article. In brief, a total of 74 multiparous Holstein cows at third parity with similar age, body condition score, and due dates (Supplementary Table 1) were enrolled at 3 weeks before parturition and raised in the same environment. All animals were fed regularly three times a day at 07:00, 13:00, and 19:00 with a total mixed ration, and the basal formulation may be found in Supplementary Table 2. Feed and water were offered ad libitum. Animals were managed by staff trained to identify medical problems and ketosis. Ketosis was determined at both prepartum (2 weeks before parturition) and postpartum (2 weeks after parturition) by testing plasma β-hydroxybutyrate (BHBA). A clinical ketosis case (case) was determined as a cow having BHBA concentration of ≥ 2.60 mM and a healthy control (control) with < 1.0 mM, respectively. Animals were excluded from the herd if diagnosed with clinical ketosis at prepartum or had other diseases during the whole experiment period. Among these cows, a subset of 19 animals were enrolled, with nine ultimately developing into clinical ketosis and 10 remaining healthy controls at postpartum.

Blood samples (10 mL) at both prepartum and postpartum were collected from coccygeal veins using vacutainer tubes before the morning feeding. After 30 min at room temperature for clot formation, samples were centrifuged at 3000 × g for 4°C for 15 min to obtain the corresponding blood serum sample. Serums were then stored in liquid nitrogen until used for further analysis. For these 19 finally enrolled animals, a total of 38 serum samples from two time points of both postpartum and prepartum were obtained and then divided into four groups: cases at postpartum (CK; BHBA = 2.79 ± 0.12 mM; n = 9), cases at prepartum (PCK; BHBA = 0.36 ± 0.05 mM; n = 9), controls at postpartum (HC; BHBA = 0.65 ± 0.22 mM; n = 10) and controls at prepartum (PHC; BHBA = 0.42 ± 0.08 mM; n = 10). For each group, all the serum samples were used for metabolomic analysis and five randomly selected samples were further subjected to proteomic analysis.

The LC-MS based metabolomic analysis has been described previously (Dunn et al., 2011; Luo et al., 2019). In brief, all 38 serum samples were slowly thawed at 4°C, and 100 μl of serum was mixed with 400 μl pre-cooled methanol. The mixed liquor was centrifugated at 12,000 rpm and 4°C for 10 min, and the supernatant was collected and blow-dried by vacuum concentration. Subsequently, the dried samples were dissolved with 150 μl of 2-chlorobenzylamine (4 ppm) methanol aqueous solution (4:1, 4°C). Finally, the supernatant was filtered through a 0.22 μm membrane and the prepared sample extracts were obtained for LC-MS analysis. For monitoring deviations of the analytical results and system stability over the entire experiment, 20 μl from each prepared sample were extracted and mixed for the preparation of quality control (QC) samples.

The liquid chromatographic separation was performed on a Thermo Ultimate 3000 system (Thermo Fisher Scientific Inc., Waltham, MA, United States) equipped with a Waters ACQUITY UPLC^® HSS T3 column (150 × 2.1 mm, 1.8 μm). The flow rate was 0.25 mL min^–1 and the column temperature was maintained at 40°C. The mobile phase consisted of 0.1% formic acid in water (A) and 0.1% formic acid in acetonitrile (B) or 5 mM ammonium formate in water (C) and acetonitrile (D). Injection of 2 μl of each sample was done after equilibration. An increasing linear gradient of solvent B (v/v) was used as follows: 0∼1 min, 2% B/D; 1∼9 min, 2%∼50% B/D; 9∼12 min, 50%∼98% B/D; 12∼13.5 min, 98% B/D; 13.5∼14 min, 98%∼2% B/D; 14∼17 min, 2%B/D. In addition, the QC sample was used to optimize the liquid chromatographic separation condition, as it contained the most information of the whole serum samples. The MS experiment was executed on the Thermo Q Exactive Focus mass spectrometer (Thermo Fisher Scientific Inc., Waltham, MA, United States) with the spray voltage of 3.8 kV and −2.5 kV in positive ion mode (ESI⁺) and negative ion mode (ESI^–), respectively. Data-dependent acquisition (DDA) MS/MS experiments were performed with HCD scan, and the normalized collision energy was 30 eV. Dynamic exclusion was implemented to remove some unnecessary information in MS/MS spectra.

The DIA large-scale proteomic method was described in detail previously (Bruderer et al., 2019). Briefly, all 20 serum samples were slowly thawed at 4°C, and then mixed with ammonium bicarbonate solution. Samples were reduced at 37°C for 1 h followed by alkylation in the dark for 1 h. Then, 100 μg of denatured serum were mixed with the ammonium bicarbonate buffer and 2.5 μg of trypsin and digested for 16 h at 37°C. Thereafter, the protein extracts were lyophilizated by freeze dryer according to the manufacture’s protocol. Desalting was performed using Sep-Pak C18 1CC Vac Cartridge (Waters, Milford, MA) following the manufacturer’s instructions. For library generation, pooled samples were fractionated using high pH reversed phase fraction chromatography (HPRP). 150 μl of digest was adjusted to pH 10 using pure ammonium formate, and then fractionated using HPRP separation on a H-Class UHPLC (Waters, Milford, MA) with a 2.1 × 150 mm BEH C18 1.7 μm column (Waters). Twelve fractions were collected and each fraction was dried in a vacuum concentrator for the next step. The fractions were resuspended with 40 μl solvent C containing 1× iRT kit, separated by nanoLC, and analyzed by on-line electrospray tandem mass spectrometry. Conditions for DDA analysis and DIA analysis were similar to those reported in Roland et al. (Bruderer et al., 2019).

The acquired MS spectra were analyzed by Mascot search engine (v.2.3.2; Matrix Science, London, United Kingdom) for protein identification by searching against the Bovine databases obtained from Uniprot¹. DTA files were generated from the raw data files and then converted to Mascot generic files using Proteome Discoverer software (v.1.4.0.288). Trypsin was specified as the proteolytic enzyme and two missed cleavage was allowed. Carbamidomethyl of cysteine was used as a fixed modification, methionine oxidation as a variable modification. The initial peptide mass tolerance was set at 10 ppm in the first search and 5 ppm in the main search, and fragment (MS/MS) mass deviation was set to 20 ppm; false discovery rate (FDR) for peptide and protein identification of all searches were less than 5%. Each protein identification involved at least one unique peptide. For protein quantification, a protein had to contain at least one unique spectra. The quantitative protein ratios were weighted and normalized by the median ratio in Mascot².

In the metabolomics analysis, all LC-MS data were extracted by ProteoWizard (v.3.0.878) and converted to mzXML format. All mass spectra were processed with peaks identification, peaks filtration, and alignment using the R package XCMS (Smith et al., 2006). The chromatographic peak data were normalized uniformly, and the multidimensional data were analyzed using SIMCA-P software (v.14.1). The metabolic peaks with relative standard deviations (RSDs) larger than 30% in QC samples were removed from the dataset (Dunn et al., 2011). Principal component analysis (PCA) was carried out to determine the global clustering and separation trends or possible outliers in an un-supervised manner. Orthogonal partial least squares discriminant analysis (OPLS-DA), a supervised model, was performed to obtain an overview of the complete data set and discriminate the inter-group differences. The model quality could be evaluated based on interpretation of variation for the X matrix (R²Y) and forecast ability of the model (Q²), which was discussed elsewhere (Yin et al., 2009). Generally, the model is believed to be reliable when Q² > 0.4. The differentially accumulated metabolites (DAMs) were screened out using variable importance projection threshold (VIP > 1.0) in the OPLS-DA model and p-value in student’s t-test (p < 0.05) (Tian et al., 2015). Identification of metabolites was carried out by searching the reference standard MS/MS spectral library or the HumanMetabolome Database (HMDB³), Metlin⁴, or mzcloud⁵ database. The functional enrichment analysis of DAMs was performed based on the Kyoto Encyclopedia of Genes and Genomes (KEGG) database using MetaboAnalyst 4.0 online tool (Chong et al., 2018).

In the proteomics analysis, Student’s t-test was used to compare protein differences between groups and to calculate p-values. Protein with a fold change of 1.5 and p < 0.05 was considered as differentially abundant protein (DAP) based on the published reference (Guo et al., 2019). We used PCA and partial least squares discriminant analysis (PLS-DA) to visualize the distribution of the samples between case and control groups and detect potential outliers. Functional enrichment analysis was performed using the differentially expressed proteins between different groups by functional categorization of Gene Ontology (GO) terms with agriGO toolkit (Tian et al., 2017), and KEGG pathway analysis was performed by KOBAS 3.0 (Xie et al., 2011).

Pearson correlation analysis was conducted for evaluating the metabolomics and proteomics integration. For this, the expression data of both DAMs and DAPs related to ketosis were calculated. Then, only the protein and metabolite with high correction (|R| ≥ 0.7) and p-value < 0.05 were considered. Finally, pathway over-representation analysis was conducted using Integrated Molecular Pathway Level Analysis (IMPaLA) (Kamburov et al., 2011).

The pooled QC sample was applied to ensure the reproducibility of the LC-MS system. The overlapped total ion chromatograms (TIC) of QC samples in positive and negative ion modes demonstrated the strong repeatability of the instruments, and more than 70% of main peaks had RSDs lower than 30% (Supplementary Figure 1). These results represented the robustness of the system. Therefore, the method was deemed acceptable for our subsequent metabolic analysis. Multivariate analysis of OPLS-DA was used to detect potential outliers and identify features potentially responsible for causing the variation between different groups of sera. Generally, the model was jointly assessed using R² (model fit) and Q² (predictive power) and the model is believed to be reliable when the R²Y and Q² values > 0.4. As shown in Figure 1, the OPLS-DA score plot could separately distinguish each of the four sera groups, indicating the differential metabolomic profiling of the four groups. Furthermore, in the positive ion mode, R²X = 0.297, R²Y = 0.999, and Q² = 0.936 (Figure 1A), whereas in the negative ion mode, R²X = 0.25, R²Y = 0.993, and Q² = 0.848 (Figure 1B). Both the R²Y and Q² values of the models were greater than 0.4, indicating that the models were predictable and reliable to discriminate among the four groups. After rigorous quality control and identification, we obtained 97 metabolites (Supplementary Table 3) among all samples.

We next examined the metabolomic profiles of sera from different groups. Firstly, we examined the variations in metabolomic profiles of cases during parturition. Multivariate analysis showed that when sera metabolites from CK were compared to those from PCK, both score plots from the PCA (Supplementary Figure 2) and OPLS-DA (Supplementary Figure 3) exhibited a clear separation without any overlap. It is noted that both the R²Y and Q² values of OPLS-DA models from positive ion mode and negative ion mode were greater than 0.4 (Supplementary Figure 3), demonstrating that the models were stable and reliable. We next examined dissimilarities in the abundance of identified metabolites between these two groups. Based on the VIP value in the OPLS-DA model > 1 and p-value in student’s t-test < 0.05, a total of 76 DAMs were identified. Of these, the level of 46 metabolites had increased, whereas those of 30 had decreased in CK group with respect to the level of PCK group (Supplementary Table 4). Sera belonging to a group (CK or PCK) clustered together perfectly (Figure 2A).

Secondly, in order to reveal the successful adaption changes in global metabolites during the transition phase, metabolomic profiles of controls during parturition were analyzed using multivariate analysis. Our results of PCA (Supplementary Figure 2) and OPLS-DA (Supplementary Figure 3) showed that sera metabolites from HC and those from PHC were clearly separated into two parts. We also observed that both the R²Y and Q² values of OPLS-DA models from both positive ion mode and negative ion mode were greater than 0.4 (Supplementary Figure 3), indicating a satisfactory effectiveness of the models which can be used to identify the difference between two groups. Based on the criteria of VIP > 1 and p < 0.05, 48 variables were screened out with 34 metabolites upregulated and 14 downregulated ones in the HC group when compared with PHC group (Supplementary Table 5). The hierarchical clustering illustrated that these metabolites clearly segregated the samples into two groups (Figure 2B).

Thirdly, we used PCA and OPLS-DA score plots to detect potential outliers and identify features potentially responsible for causing the variation between cases and controls at postpartum. Our results of PCA (Supplementary Figure 2) and OPLS-DA (Supplementary Figure 3) also showed excellent separation between CK and HC. The parameters of the OPLS-DA models from positive ion mode and negative ion mode showed that both the R²Y and Q² values were greater than 0.4 (Supplementary Figure 3). This indicated that these are reliable and predictable models to discriminate between the two groups. To identify the DAMs between the CK and HC groups, we compared the abundance of identified metabolites between these two groups. A total of 31 DAMs were obtained from the comparison, 14 of which had a higher relative abundance in the CK than HC group; the other 17 metabolites significantly decreased in the CK group (Supplementary Table 6). Cluster hierarchization showed that the clusters of these two groups were obviously separated (Figure 2C).

In addition, we found 21 shared metabolites between the comparisons of CK vs. HC and CK vs. PCK (Figure 3A). Strikingly, of these shared metabolites, 10 were consistently upregulated and the other 10 consistently downregulated in the CK group; only one metabolite of L-Pipecolic acid did not show agreement in the direction of the fold change between comparisons of CK vs. HC and CK vs. PCK (Figure 3B), indicating that these 20 metabolites were consistent with the clinical determination of the cases. The most upregulated metabolites in the CK group included the metabolites of 4-Hydroxy-6-Methylpyran-2-one, BHBA, and cinnamoylglycine (Figure 3).

We evaluated the interactions of the 76 DAMs identified between CK vs. PCK, which would help to unravel potential metabolic changes contributed to ketosis from prepartum to postpartum. The enrichment analysis revealed that these DAMs were enriched in 32 KEGG pathways. Among these pathways, carbohydrate metabolism (35%), amino acid metabolism (22%), and metabolism of cofactors and vitamins (13%) accounted for a large proportion (Figure 4A). Furthermore, a total of five pathways of “D-Glutamine and D-glutamate metabolism”, “Alanine, aspartate, and glutamate metabolism”, “Arginine and proline metabolism”, “Histidine metabolism”, and “Citrate cycle (TCA cycle)” were the most enriched pathways (p < 0.05, impact value > 0.10) (Figure 4B). For the 48 DAMs identified between HC vs. PHC, these would help to unravel the metabolic pathways involved in successful adaption from prepartum to postpartum. The enrichment analysis revealed that these DAMs were enriched in 20 KEGG pathways. Among these pathways, amino acid metabolism (30%), carbohydrate metabolism (25%), and metabolism of cofactors and vitamins (15%) accounted for a large proportion (Figure 4C). The functional impact pathways were shown in Figure 4D, by which we found only the pathway of “Arginine and proline metabolism” was significantly changed in lactating cows.

To further investigate metabolic variations of cases and controls at postpartum, we mapped the 31 DAMs identified between CK vs. HC to the KEGG database. The results of enrichment demonstrated that there were 22 pathways enriched. Among these pathways, amino acid metabolism (36%), carbohydrate metabolism (23%), and metabolism of cofactors and vitamins (14%) accounted for a large proportion (Figure 4E). The five significantly enriched pathways included “D-Glutamine and D-glutamate metabolism”, “Valine, leucine, and isoleucine biosynthesis”, “Alanine, aspartate, and glutamate metabolism”, “Histidine metabolism”, and “Arginine and proline metabolism” (Figure 4F).

To further understand the systematic changes associated with metabolic adaptions to transition phase stress as well as metabolic variations associated with ketosis during the transition period, the proteomic profiles among cases and controls at both prepartum and postpartum were carried out. Of the total sequenced spectra, 7,036 were mapped to the bovine reference protein database. Of these, 4,072 were uniquely mapped to specific peptides. In total, 540 proteins were identified under the 5% false discovery rate threshold at both the peptide and protein levels. Next, all the 540 proteins were subjected to PCA analysis for clustering all the samples. The first two principle components accounted for 30% of total variance (Figure 5A). This panel did not show a clear separation between sera from HC and PHC groups. However, it can be easily observed that sera within each group tended to cluster together (Figure 5A). Subsequently, the PLS-DA model was used for further multivariate analysis, which revealed that the proteomics of each group could be clearly distinguished from the others (Figure 5B), indicating the differential proteomic profiling of the four groups of sera. In summary, clear separation in proteomic profiles was found among sera from cases and controls.

To survey the DAPs associated with transition phase stress as well as ketosis, we compared the abundance of identified proteins among the following three comparisons. Firstly, we identified 37 DAPs in sera between the CK vs. PCK (fold change > 1.5 and p < 0.05), of which 16 were upregulated and 21 were downregulated in the CK group (Supplementary Table 7). A heatmap of these DAPs was generated to visualize expression patterns across all 10 sera samples, and expression patterns in the heatmap were accompanied by hierarchical clustering of proteins (horizontal axis) and samples (vertical axis) (Figure 6A). Secondly, we identified 30 DAPs, with 10 significantly higher and 20 significantly lower relative concentration proteins in the HC group when compared with PHC group (Supplementary Table 8). Cluster hierarchization using the expression data of these 30 DAPs confirmed the presence of two distinct groups (Figure 6B). Thirdly, a total of 30 DAPs was identified between CK vs. HC, with 18 upregulated and 12 downregulated in CK group compared with HC group (Supplementary Table 9). The protein abundance data of these 30 DAPs revealed an obvious separation of two parts (Figure 6C), which confirmed the presence of discriminating features between CK and HC groups.

The functional enrichment analysis of DAPs based on GO categories was performed and the significantly (p-adjusted < 0.05) enriched terms were shown in Figures 7A–C. By which, the 37 DAPs obtained between CK vs. PCK were enriched in 40 GO terms. Three of these terms corresponded to molecular function (MF), namely protein binding, enzyme regulator activity, and enzyme inhibitor activity, and most of these proteins were concentrated in protein binding (Figure 7A). For biological process (BP) ontology, 22 terms were enriched, and most of these proteins were enriched in the terms of biological regulation, regulation of biological quality, and regulation of molecular function (Figure 7A). The cellular component (CC) ontology presented 15 enriched terms, and those of extracellular region, extracellular region part, membrane-bounded organelle, and extracellular region were ranked at the top of the category (Figure 7A). In addition, a total of 36 GO terms were obtained by DAPs between HC vs. PHC. For biological processes, most proteins were enriched in multicellular organismal process, localization, developmental process, and response to stimulus (Figure 7B). For molecular function, the top term was protein binding (Figure 7B), while no MF terms were enriched by DAPs between CK vs. HC (Figure 7C). Meanwhile, the proteins participate in several CC terms, such as extracellular region, extracellular region part, organelle, and membrane-bounded organelle, and these terms had high ratios among the DAPs that were identified among all three comparisons within each comparison.

For cellular components, most proteins were enriched in the extracellular region, extracellular region part, and membrane-bounded organelle. For molecular function, the top term was protein binding (Figure 7B), while no MF terms were enriched by DAPs between CK vs. HC (Figure 7C).

The KEGG analysis of DAPs allowed us to better understand the key proteins and pathways affected by transition phase stress and ketosis. By which, the KEGG pathways enriched by DAPs between different comparisons were displayed in Figures 7D–F. A total of four pathways, including amebiasis, vitamin digestion and absorption, phagosome, and legionellosis, were significantly (p-adjusted < 0.05) enriched by DAPs between CK vs. PCK (Figure 7D); those between HC vs. PHC were cholesterol metabolism and phagosome (Figure 7E), and those between CK vs. HC were ECM-receptor interaction (Figure 7F).

To investigate the protein and metabolite regulatory network of ketosis implicated in transition cows, we performed a pathway over-representation analysis using IMPaLA tool. By which, the DAMs and DAPs obtained from both comparisons of CK vs. PCK and CK vs. HC were used to identify significantly perturbed pathways. In total, 85 metabolites and 53 annotated proteins were subjected to Pearson correlation analysis (Supplementary Table 10). The results showed that 38 metabolites had 78 strong corrections (R ≥ 0.7) with 25 proteins (Figure 8).

Based on the results, interaction networks between the 38 metabolites and 25 proteins were organized. The list of pathways with joint p-values < 0.05 is shown in Supplementary Table 11. Ketosis altered a wide range of metabolic pathways, among them mainly metabolism, metabolism of proteins, and metabolism of angiotensinogen to angiotensins. The cellular processes involved were ones such as gene expression (transcription), RNA polymerase II transcription, and post-translational protein modification. Also involved were vitamin metabolism (e.g., vitamin digestion and absorption, metabolism of vitamins and cofactors, metabolism of water-soluble vitamins and cofactors) and signaling (e.g., G alpha (q) signaling events, GPCR ligand binding, GPCR downstream signaling, signaling by GPCR, and signal transduction). Other pathways associated with the vascular system (hemostasis, ion channels, and their functional role in vascular endothelium) and disease (amebiasis) were identified.