Metabolomic profiles in yak mammary gland tissue during the lactation cycle

The yak is one of the most important domestic animals in Tibetan life for providing basic resources such as milk, meat and transportation. Although yak milk production is not elevated, yak milk is superior to dairy cow milk in nutrient composition (protein and fat). However, the understanding of the metabolic mechanisms of yak mammary gland tissue during the lactation cycle remains elusive. In this study, GC-MS-based metabolomics was employed to study the metabolic variations in the yak mammary gland during the lactation cycle (pregnancy, lactation and dry period). Twenty-nine metabolites were up or downregulated during the lactation period. Compared to the dry period, during the lactation period the levels of oxalic acid were upregulated, while glycine and uridine were downregulated. Thirty-seven pathways were obtained when the 29 differential metabolites were imported into the KEGG pathway analysis. The most impacted pathways during the lactation cycle were glycine, serine and threonine metabolism; alanine, aspartate and glutamate metabolism; TCA cycle; glyoxylate and dicarboxylate metabolism; and pyrimidine metabolism. Our results provide important insights into the metabolic events involved in yak mammary gland development, lactogenesis and lactation, which can guide further research to improve milk yield and enhance the constituents of yak milk.


Introduction
The yak (Bos grunniens) lives in regions centered on the Chinese Qinghai-Tibet Plateau, including the highlands of Nepalese Himalayas, Indian Kasmir, Tibet, Mongolia and Bhutan [1] at high altitudes from 2000-5000 m where few other animals can survive. The number of yaks that live in China is approximately 13 million, comprising 92.8% of the gross number of yaks on this planet [2]. Yaks can use the pasture resources of this area and provide transportation, leather, meat and milk. Yak milk and milk products are the major dietary ingredients for 6.5 million Tibetan people as well as it is an important source of income to the local families. For these reasons, the yak is one of the most important domestic animals in Tibetan life [3], and the yak daily industry has rapidly grown in recent years [4]. a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 Although yaks have a limited production of milk (~3 kg/day), their milk is substantially richer in nutrients than dairy cows [5]. Yak milk contains 16.9%-17.7% solids, including 4.9%-5.3% protein, 5.5%-7.2% fat, 4.5%-5.0% lactose, and 0.8%-0.9% minerals [5]. A better understanding of the inherent metabolic status during different stages of lactation in yak mammary gland can enlighten important metabolic processes, which determines milk production and composition.
Metabolomics can be a powerful approach to investigate the low-molecular-weight metabolites and provide a global understanding of physiological alterations in specific organs/tissues [6]. In a previous study, the metabolite profile of yak milk was evaluated [7]; however, there is a lack of information regarding the metabolic profile of yak mammary gland. Considering the different physiological alterations of mammary gland during the lactation cycle, we hypothesized that the metabolite profile of yak mammary gland will also change over the lactation stage. The objective of the present study was to investigate the metabolite profile of yak mammary gland tissue during pregnancy, lactation and dry period by gas chromatography/mass spectrometry (GC-MS).

Ethics statement
The animal care and use were performed according to the regulations for the Administration of Affairs Concerning Experimental Animals (Ministry of Science and Technology, China, revised in June 2004) and approved by the Institutional Animal Care and Use Committee, Southwest Minzu University, Chengdu, Sichuan, China.

Yaks and sample collection
The surgical procedure for yak mammary biopsy was performed in accordance with the operations guide to ameliorate animal suffering. Seven healthy Maiwa yaks (4 years old) with unrelated background were obtained from Longri Breeding Farm of Sichuan Province, Hongyuan County, Sichuan Province, China. Individuals examined from pregnancy to the dry period were divided in to nine subgroups, class T1 (n = 5, -30±2 d), class T2 (n = 6, -15±2d), class T3 (n = 6, +1 d), class T4 (n = 7, +15 d), class T5 (n = 7, +30 d), class T6 (n = 6, +60 d), class T7 (n = 7, +120 d), class T8 (n = 7, +180 d), and class T9 (n = 6, +240 d) according to their parturition (d). Class T1 and T2 belonged to the pregnancy period, Class T3 to T8 belonged to the lactation period, Class T9 belonged to the dry period. A total of 57 mammary gland tissue samples were collected in the present study. Mammary gland tissue samples (approximately 500 mg) were collected by biopsy of the right or left rear quarters as previously described [8]. All samples were immediately frozen and stored in liquid nitrogen.

Metabolite profiling analysis
Mammary gland tissue samples (~100 mg) was ground and transferred to 2-mL centrifuge tubes. A thousand μL of 80% methanol (pre-cooled at -20˚C) was then added and vortexed for 30 s, followed by the addition of 60 μL of 2-Chloro-L-phenylalanine (0.2 mg/mL stock in methanol) and 60 μL of heptadecanoic acid (0.2 mg/mL stock in methanol) as internal quantitative standards; the mixture was vortexed for 30 s. The tubes were placed into an ultrasound machine at room temperature for 30 min, incubated on the ice for 30 min, and centrifuged at 14,000 rpm at 4˚C for 10 min; subsequently, 0.8 mL of supernatant was transferred to a new centrifuge tube. The samples were blow-dried by moderate nitrogen, and 60 μL of 15 mg/mL methoxyamine pyridine solution was then added, followed by vortexing for 30 s and incubation for 120 min at 37˚C. Sixty μL of BSTFA reagent (containing 1% TMCS) was added to the mixture, followed by incubation for 90 min at 37˚C. Finally, the samples were examined for the metabolite contents by using GC-MS. A quality control (QC) sample was also prepared by mixing equal volumes of each sample; the samples were aliquoted for analysis prior to sample preparation. The QC samples were used to monitor deviations of the analytical results from the pooled mixtures and compare to the errors caused by the analytical instrument itself.
GC-MS data were acquired using Agilent 7890A GC system coupled to an Agilent 5975C inert XL EI/CI mass spectrometric detector (MSD) system (Agilent Technologies, Santa Clara, CA, USA). GC was performed on an HP-5MS capillary column (5% phenyl/95% methyl polysiloxane [30 m × 250 μm i.d., 0.25 μm film thickness, Agilent J & W Scientific, Folsom, CA, USA]) to separate the derivatives at a constant flow of 1 mL/min helium. A 1-μL aliquot of the sample was injected in split mode in a 20:1 split ratio by the auto-sampler. The injection temperature was 280˚C, the interface was set to 150˚C and the ion source was adjusted to 230˚C. The temperature-rise programs included an initial temperature of 80˚C for 5 min, followed by a 20˚C/min increased up to 300˚C and a hold at 300˚C for 6 min. Mass spectrometry was determined by the full-scan method with a range from 35 to 500 m/z.

Data analysis
Raw GC-MS data were converted into CDF format (NetCDF) files by Agilent GC/MS 5975 Data Analysis software and subsequently processed by XCMS (www.bioconductor.org) [9,10]. The signal integration area of each metabolite was normalized with the total peak area for each sample. The identification of metabolites using the Automatic Mass Spectral Deconvolution and Identification System (AMIDS) was searched against commercially available databases such as the National Institute of Standards and Technology (NIST) and Wiley libraries. Metabolites were confirmed by comparison of mass spectra and retention indices to the spectra library using a cutoff of 70% identity [11].
For multivariate statistical analysis, the XCMS output was further processed using Microsoft Excel (Microsoft, Redmond, WA, USA). Finally, the normalized data were imported into Simca-P software (version 13.0, http://www.umetrics.com/simca) for multivariate statistical analyses, including principal component analysis (PCA) and partial least-squares discriminant analysis (PLS-DA). All data were mean-centered and unit variance (UV)-scaled before PCA and PLS-DA applied to avoid overfitting. In the present study, a default 7-fold (Leave-1/7th samples-Out) cross validation procedure and 100 random permutations tests were performed to avoid overfitting of the supervised PLS-DA models.
Hierarchical cluster analysis (HCA) was performed and visualized by the "pheatmap" package available in R software (www.r-project.org). Pearson's product-moment correlation was performed to calculate the correlation. Corresponding P-values and false discovery rate (FDR) of each correlation were also calculated using "cor.test function" [12] in R software. Identified metabolites were mapped onto general biochemical pathways according to annotation in Kyoto Encyclopedia of Genes and Genomes (KEGG) [13].

Mammary gland tissue metabolic profiles
In total, 80 retention time−exact mass pairs remained in each sample profile using the GC-MS analysis protocol and subsequent processes. As shown in S1 Fig. no drift was observed in the PCA scores plot representation of QC samples. Thus, the metabolic features demonstrated acceptable reproducibility and stability in GC−MS profiling analyses. Changes in the metabolite profile over the lactation cycle were observed. In addition, a plot of scores at different time phases on the axis of the first principal component (PC1, R 2 X 1 = 0.30) accounted for a sizable portion of variance during the pregnancy, lactation and dry periods ( Fig 1A). Subsequently, the PLS-DA models were established and demonstrated satisfactory modeling (R 2 X = 0.45, R 2 Y = 0.15, Q 2 = 0.09), which showed differences of the nine classes during lactation in the direction of the first and second components (Fig 1B). The correlation of 80 metabolites was calculated by "Pearson" correlation coefficient, and the results are showed in Fig 2.

Metabolic variations during different periods
Based on the mammary gland tissue metabolic profiles, the scores plots and permutation tests of PLS-DA model discriminating the pregnancy period from the lactation period are presented in Fig 3, the results show a clear separation between the two groups. Mammary gland tissue metabolites passing the VIP threshold (VIP > 2) in this model (R 2 X = 0.41, R 2 Y = 0.77, Q 2 = 0.67) were selected as significantly different between the pregnancy and lactation periods (Fig 3A and 3B). Similarly, another two-component PLS-DA model (R 2 X = 0.51, R 2 Y = 0.69, Q 2 = 0.34) was also constructed to discriminate the dry period from the lactation period, achieving a clear discrimination between the two groups (Fig 3C and 3D). Although the Q 2 value between the dry period and lactation period was slightly lower, this value was acceptable, considering the many uncontrollable factors. Model cross-validation through permutation tests (100 times) generated intercepts of R 2 and Q 2 (0.30 and -0.21 for pregnancy and lactation period, 0.48 and -0.24 for dry and lactation period, respectively). The low values of the intercepts indicate that the model was not over-fitted.
A list of 29 differential expressed metabolites (Table 1) including beta-alanine, aspartic acid, fructose, xylitol, arabitol, glycine and myo-inositol were identified by commercially available reference standards with MS/MS fragments and retention time (P < 0.05). As shown in Table 1, the relative concentration of 19 metabolites were significantly lower while 5 were significantly higher in the lactation period compared to pregnancy period. The relative concentration of 4 metabolites were significantly lower and 4 were significantly higher in the dry period compared to lactation period. These metabolites were carbohydrates, amino acids, lipids and their metabolites, suggesting that these metabolic pathways were altered through the lactation cycle. Compared to the lactation period, oxalic acid concentrations were higher during both pregnancy (P < 0.001, VIP = 1.17, fold-change = 0.51) and dry periods (P < 0.05, VIP = 1.59, fold-change = 1.78) (Fig 4). Glycine and uridine showed decreased levels during both the lactation and dry periods of the lactation cycle (lactation period VS. pregnancy period, P < 0.01, fold change: 0.59, 0.60, respectively; dry period VS lactation period, P < 0.05, fold change: 0.60, 0.54, respectively) (Fig 4).

Metabolic KEGG pathway analysis
Overall, 37 pathways were obtained when the 29 differential metabolites were imported into the KEGG analysis.  most enriched functional pathways between lactation period and pregnancy period included glycine, serine and threonine metabolism (FDR < 0.001); alanine, aspartate and glutamate metabolism (FDR < 0.001); and TCA cycle (FDR = 0.008). The most enriched functional pathways between lactation and dry periods were glycine, serine and threonine metabolism (FDR < 0.001); glyoxylate and dicarboxylate metabolism (FDR = 0.01); and pyrimidine metabolism (FDR = 0.02). Among the 29 metabolites affected by the treatments, 5 are involved in the glycine, serine and threonine metabolism; 6 in the alanine, aspartate and glutamate metabolism; 4 in the TCA cycle; 9 in the glyoxylate and dicarboxylate metabolism; and 6 in the pyrimidine metabolism.

Discussion
The defining characteristic of the mammalians is the provision of milk. Milk composition differs greatly across species, although species in the same taxonomic order tend to produce milk with a somewhat similar composition. The mammary gland is one of the few tissues in the body that undergoes repeated cycles of structural development, functional differentiation, and metabolism. However, these changes, particularly those related to the metabolism of the mammary gland, are still relatively unexplored. The inherent metabolic status during the lactation cycle play a vital role in the lactation physiology, which determines the milk production, udder health and dairy sustainability [14]. To our knowledge, this is the first paper addressing the lactation-related metabolic profiles using mammary gland tissue of yaks. This data provides novel insights and a greater understanding of the metabolic mechanism during the lactation cycle of yaks, and can guide further research to improve milk yield and enhance the constituents.
As shown in Fig 3A and 3B, PLS-DA scores plot revealed a clear and statistically significant separation of the lactation period versus the pregnancy period and the dry period versus the lactation period. A total of 29 differential metabolites were identified throughout the lactation cycle, and these differential metabolites were involved in pathways including amino acid, biosynthesis of other secondary metabolites, carbohydrate, energy, lipid, cofactors and vitamins and nucleotide metabolism, indicating that these metabolic pathways were influenced by different periods of the lactation cycle. The significant changes observed in mammary gland tissue composition throughout lactation cycle are primarily related to the identified compounds: oxalic acid, glycine, and uridine. Oxalic acid has usually been seen as an inert end product of the metabolism in animals because mammalian cells cannot metabolize it [15]. However, recent studies indicate that the levels of oxalate are too high for the substance to only be an end-product of the metabolism in animals [16]. Therefore, it has been suggested that there could be an oxalate oxidase pathway in animals which uses oxalate to promote phagocytes [17]. Oxamide is a diamide that is derived from oxalic acid and participates in the TCA cycle, and the TCA cycle is one of the most important pathways in the mammary grand for lactation initiation [18,19]. The expression of oxalic acid was showing significantly lower in the mammary grand tissues of lactating than non-lactating cows [18] which is consistent with ours. Citric acid (fold-change = 5.50) and other secondary metabolites participating in the TCA cycle were significantly higher in the lactating yaks, indicating that the TCA cycle is more active during lactation initiation. The changes of the expression of oxalic acid was presumed to be related to the up-regulated TCA cycle.
Glycine is the simplest amino acid and a basic nutrient. As part of endogenous antioxidant glutathione, glycine is a semiessential amino acid. Glycine is primarily involved in the proliferation of immune cells and synthesis of proteins related to inflammatory and immune responses [20][21][22][23]. Additionally, glycine can also regulate the proliferation and apoptosis of mammary epithelial cells of dairy cows and promote the conversion of mammary epithelial cells to the S phase and cell division [24]. Two studies have reported that glycine was the most abundant metabolite in the goat mammary gland [25] and lactating dairy cow mammary grand [18]. Glycine metabolism pathway is important in the mammary gland for lactation initiation, because glycine can be used to synthesize creatine with arginine and methionine. It is reported that more creatine is needed to ease the serious negative energy balance in dairy cows [26]. We speculate that the decreased level of glycine was related to the negative energy balance in yaks.
Uridine is one of the five standard nucleosides that make up nucleic acids. The relative concentration of uridine was significantly lower in Holstein than yak and other dairy animals [7]. Besides, the level of uridine was identified significantly higher in the lactation than non-lactation cows and that was consistent with our result [18]. Uridine metabolism can be regulated by the prolactin stimulation, the uridine uptake and its incorporation in to RNA are concomitantly stimulated by prolactin in mammary gland [27]. In mammals, prolactin stimulates the mammary gland to produce milk and has important cell cycle-related functions such as a growth factor and differentiating factor [28]. The prolactin levels increase concentration during pregnancy cause enlargement of the mammary glands and preparing for milk production and fall when the process of breastfeeding is terminated. Therefore, we speculate that the decreased level of uridine was affected by the decreased prolactin level.
In our study, all of the shared metabolites in the mammary grand from the pregnancy, lactation and dry periods were evaluated together. For the lactation yaks, the most active and important pathways included glycine, serine and threonine metabolism, alanine, aspartate and glutamate metabolism, TCA cycle, glyoxylate and dicarboxylate metabolism and pyrimidine metabolism. Among them, TCA cycle is the key pathways of carbohydrate metabolism and plays a central role in cellular respiration and the supply of energy to all living cells [29], which is of paramount significance to the cell's metabolic efficiency and, therefore, to the yak's metabolism and production. The differential expression level of citric acid between the lactating and non-lactating cows (fold-change = 5.50) [18] was consistent with our results (foldchange = 3.43). Citric acid is formed in the TCA cycle or from the diet and participates in the intermediate metabolism of carbohydrate oxidation in animal tissues [30]. The level of citric acid is significantly higher in the lactating yaks and may enhance energy by participating in the TCA cycle. TCA cycle and glycine, serine and threonine metabolism pathways are the most important and work together in the mammary grand for lactation initiation [18,31]. As a prime metabolic source of glutathione, creatine, purines and serine and a protein constituent in the lactating mammary epithelial cells, glycine plays a significant role in various biological processes such as ischemia reperfusion injury, oxygen stress, cell membrane injury, tumor metastasis and so on [32,33]. Glyoxylate and dicarboxylate metabolism pathway were determined to be pathways that were significantly impacted in the lactation cows [12]. In the alanine, aspartate and glutamate metabolism pathway, aspartate is a precursor of many compounds that are involved in cellular signaling, such as N-acetyl-aspartate. It is also a metabolite in the urea cycle and participates in gluconeogenesis in dairy cows [34]. Pyrimidines are involved in growth via DNA and RNA formation and they also have a special role in mammary grand because of the essential role of uridine in the synthesis of lactose and other uridinesugar complexes, and the vital role of cytidine in the synthesis of complex lipids resulting from the mechanism of milk secretion [35]. Thymine is broken down into β-aminoisobutyrate which can be further broken down into intermediates eventually leading into the TCA cycle. Based on the overview map of the KEGG pathway, these 5 key metabolic pathways were closely integrated together, suggesting that the carbohydrate pathways, amino acid pathways and nucleotide pathway play vital roles in regulating lactation maintenance for the yak.

Conclusion
The present study investigates the metabolite profile of yak mammary gland tissue during pregnancy, lactation and dry periods by gas chromatography/mass spectrometry. The analysis revealed substantial metabolic alteration through the lactation cycle; 29 metabolites were up or downregulated during the lactation cycle. Thirty-seven pathways were obtained when the 29 differential metabolites were imported into the KEGG analysis; the main pathways affected by the lactation cycle were glycine, serine and threonine metabolism; aspartate and glutamate metabolism; TCA cycle; glyoxylate and dicarboxylate metabolism; and pyrimidine metabolism. Overall, our results provide a better physiological understanding of the lactation metabolism of yaks, which can help elucidate the regulated metabolic strategies for the lactation yaks in the future.