Comparative Analysis of Milk Triglycerides Profile between Jaffarabadi Buffalo and Holstein Friesian Cow

Milk lipids are known for a variety of biological functions, however; little is known about compositional variation across breeds, especially for Jaffarabadi buffalo, an indigenous Indian breed. Systematic profiling of extracted milk lipids was performed by mass spectrometry across summer and winter in Holstein Friesian cow and Jaffarabadi buffalo. Extensive MS/MS spectral analysis for the identification (ID) of probable lipid species using software followed by manual verification and grading of each assigned lipid species enabled ID based on (a) parent ion, (b) head group, and (c) partial/full acyl characteristic ions for comparative profiling of triacylglycerols between the breeds. Additionally, new triacylglycerol species with short-chain fatty acids were reported by manual interpretation of MS/MS spectra and comparison with curated repositories. Collectively, 1093 triacylglycerol species belonging to 141 unique sum compositions between the replicates of both the animal groups were identified. Relative quantitation at sum composition level followed by statistical analyses revealed changes in relative abundances of triacylglycerol species due to breed, season, and interaction effect of the two. Significant changes in triacylglycerols were observed between breeds (81%) and seasons (59%). When the interaction effect is statistically significant, a higher number of triacylglycerols species in Jaffarabadi has lesser seasonal variation than Holstein Friesian.


Introduction
Milk is a complete nutritious drink composed of proteins, carbohydrates, fats, vitamins, and a range of bioactivities influencing the human health of all age groups in a positive way [1,2]. Milk composition is influenced by several factors such as breed, environmental conditions, lactation stage, and animal's nutritional status [3,4]. To date, the majority of lipidomics studies have focused on improving dairy nutritional management strategies to enhance milk quality and production by improving pasteurization and preservation techniques to extend milk shelf-life and to remove bacteria, spores, and somatic cells from milk [5]. However, the role of different exogenous (e.g., season) and endogenous (e.g., breed) factors influencing milk composition, particularly lipids, a potential source of functional food component of milk and dairy products, has been limited in both cows and buffaloes [6,7]. breed-and season-dependent, and consideration of this information is critical for estimating milk's nutritional value as well as to identify potential markers of milk fraud impacting human health.

Cumulative Lipid Identifications
The UPLC/MS-based lipidomics relies on compressive data analysis for success. A variety of commercial packages and freeware approaches have been developed and implemented with vendor-neutral software allowing analysis of data acquired from any mass spectrometer software [31]. Here, we use SimLipid's inbuilt functionality for automated assignment of MS/MS spectra with possible lipid species (detailed description on the database search parameters and filters provided in Appendix A.1) together with a manual interpretation of MS and MS/MS spectra for verification of the SimLipid reported lipid ion species as: i.
Group-ID: lipid IDs at sum composition level with only class-specific head group identification (e.g., PC 32:0). iii.
PA-ID (partial-acyl identification): For example, PC 14:0*_18:0 where the suffix * indicates that the MS/MS spectrum does not feature characteristic ion/s corresponding to the fatty acid chain 14:0. iv.
FA-ID (lipid IDs with full-acyl identification) [32]: For example, PC 14:0_18:0, indicating that the MS/MS spectrum features all the characteristic ion/s corresponding to the head group/neutral loss of head group along with the fatty acid chains.
A detailed description of the process of manual verification and grading of the lipid ion species reported by SimLipid Database search along with example spectrum lists annotated with fragment ion species is provided (Appendix A.2). Although our goal of the study is to achieve in-depth profiling of triacylglycerols (TG), we have reported phosphatidylcholine (PC), phosphoethanolamine (PE), sphingomyelins (SM), sterol esters, cholesterol and derivatives (Cho Der), and diacylglycerols (DG) species that ionized well in positive mode of MS experiment.
We present the summary of the SimLipid MS/MS database search results in Table 1. SimLipid database search using LIPIDMAPS (LM) databases reported 1132 lipids ion species that consist of 9 Mass-ID, 25 Group-ID, 251 PA-ID, and 848 FA-ID (Table 2 for detailed breakdown). Of note, only the TG and DG lipid ions at the FA-ID level make it to the final list as an effort to retain only the highly confident IDs for our further studies (we removed 243/7 TG/DG species identified at PA-ID level from the final list, respectively). Hence, we have 882 non-redundant unique lipid species that we consider for further comparative and relative quantitative analysis.
The likelihood of MS/MS spectra not assigned with any possible lipid species by SimLipid software is high, because SimLipid relies on a database containing lipid species migrated from LIPIDMAPS (LM) (https://lipidmaps.org/) that does not have many of the glycerolipids with short-chain fatty acids, e.g., 2-8 carbons, that have been reported in the published literature on fatty acid profiles of bovine milks. Hence, those unassigned MS/MS spectra were manually interpreted for the identification and characterization of novel glycerolipids (Appendix A.3).

Characterization of Novel DG and TG Lipid Species that Are Not Yet Cataloged in the LIPIDMAPS Database
The process of manually interpreting MS/MS spectra (Appendix A.3) resulted in the identification of 260 TG and one DG species that have not yet been cataloged by the LM database. Besides, we have reported nine phospholipids by manually investigating the presence of LC-compounds corresponding to ions of PC/PI/PS with their fatty acyl compositions made up of palmitic acid, stearic acid, and oleic acid between the replicates (detailed processed explained in the next section). Finally, we obtained a total of 1145 unique graded IDs comprising of 882 lipid species assigned by the initial SimLipid database search, 260 TG species and one DG species by manual interpretation of MS/MS spectra, and nine phospholipids detected using LC-compounds. Figure 1A shows the breakdown of the number of lipid molecular species belonging to each subclass. Figure 1B shows the overlay of total ion chromatograms (TIC) of all the runs annotated with the periods in which lipid ions from different classes eluted in the LC domain. Table 1. Summary of the SimLipid MS/MS database search identified using 15/10 ppm error tolerances for precursor and product ions and grading the lipid molecular ions based on observed diagnostic ions in the MS/MS spectra using MS excel.  Besides, we have reported nine phospholipids by manually investigating the presence of LCcompounds corresponding to ions of PC/PI/PS with their fatty acyl compositions made up of palmitic acid, stearic acid, and oleic acid between the replicates (detailed processed explained in the next section). Finally, we obtained a total of 1145 unique graded IDs comprising of 882 lipid species assigned by the initial SimLipid database search, 260 TG species and one DG species by manual interpretation of MS/MS spectra, and nine phospholipids detected using LC-compounds. Figure 1A shows the breakdown of the number of lipid molecular species belonging to each subclass. Figure 1B shows the overlay of total ion chromatograms (TIC) of all the runs annotated with the periods in which lipid ions from different classes eluted in the LC domain.

Comparative Analysis of Lipids between Holstein Friesian Cow and Jaffarabadi Buffalo across Season
Bovine milk contains a diverse set of lipids that perform multiple functions and that have evolved due to sustained pressure on livestock to increase milk production by adopting and implementing appropriate animal husbandry practices. Nevertheless, there is limited information on the influence of endogenous (e.g., cow vs. buffalo) and exogenous (e.g., season) factors on change in the composition of lipid species.  1  0  0  0  0  TG  39:0  1  0  0  0  0  TG  48:1  1  1  1  1  1  TG  54:5  1  1  1  1  0   TG  24:0  1  1  0  0  0  TG  39:1  1  0  0  0  0  TG  48:2  1  1  1  1  0  TG  54:6  1  0  1  0  0   TG  24:1  1  0  0  0  0  TG  39:2  1  0  0  0  0  TG  48:3  1  1  1  1  0  TG  55:0  1  0  0  0    As a group, 1145 lipids graded IDs were identified in HF cow and J buffalo. Of these lipids, 187 and 170 were detected exclusively in each group, while the remaining 788 lipids were common to both animal groups (Figure 2A), although a lipid species may not be present in a particular season. A total of 496 lipid species were detected across breed as well as season. For instance, lipid species, such as Cho Der-[C33H58O6+Na]+, were detected only in the JW, but not detected in JS, HFS, and HFW; similarly, DG 12:0_16:0_0:0 and DG 16:0_20:1_0:0, DG 18:1_18:2_0:0 were identified in the winter season of HF and J, respectively. TG 2:0_6:0_20:1, TG 3:0_9:0_13:0 lipid species present in JW were found to be made up of short fatty acyl chain while long fatty acyl chains were present in JS, HFW. Extensive breed and season associated changes were observed in different lipid subclasses (e.g., triacylglycerols and phospholipids) in both animal groups. Additionally, new triacylglycerol lipid species have been reported for the first time by manually interpreting the MS/MS spectra and comparing the detected lipid species against a curated repository of published lipids and HFS. The variation and higher numbers of TG species were detected in the winter season compared to summer season, possibly due to availability and feeding of fresh seasonal grass and breed genetic factors [7]. Similarly, sphingomyelin (SM d43:1) was detected exclusively in the summer season. The lipid classes of TG, DG, PL, identified in HF cow were found to be consistent with the previous reports of [20], while lipid IDs from J buffalo have been reported for the first time.

Variation in Unsaturation between Holstein Friesian Cow and Jaffarabadi Buffalo
The profile of unsaturated lipids obtained from this study followed the same trend that has been reported in HF. For example, C14:1, C16:1, C18:1, C18:2, C22:1 unsaturated lipids were observed in HF during the winter season. The pattern of increasing the unsaturation level in lipids during the winter season was consistent with previous studies of lipid profiling across the season in HF [35]. The DG and TG species were detected across both the seasons in HF, while phospholipids such as PE (e.g., 18:1_18:1) and PC (e.g., 34: 1) were exclusively present in HF during the winter season. The observation of lipids exclusively in the winter season could be due to the nutritional availability of greener pasture and feed, which contribute to the availability of lipid precursors [33]. The results of our study showed a similar expression pattern to that observed for saturated and unsaturated fatty acids across both seasons, as reported previously in HF cow milk [36,37]. However, an important factor that has not been adequately controlled and is not clear is the influence of physiological and It is interesting to note that 653 lipid species were common in both animal groups across the winter season, while 236 and 185 were exclusively detected in HF and J, respectively. Similarly, comparison across summer season in both groups resulted in the detection of 590 lipid species common in both groups with 105 and 163 species observed exclusively in HF and J, respectively. Similarly, multiple comparisons between animal groups across both seasons lead to the detection of 20 unique lipid species in HF in summer and 143 in winter, while 39 and 93 lipid species were observed in JS and JW, respectively ( Figure 2B). In contrast, 496 lipid species were common to both groups across the summer and winter seasons. Interestingly, 38 unique lipids were found to be common between J winter and summer season while only 24 unique lipids (e.g., TG 12:0_14:1_14:1, TG 20:0_20:0_20:1, TG 15:1_15:1_19:0, SM d32:1, SM d38:2, SM d42:2) were common across seasons in HF. The identified lipids for HF and J are listed in (See Supplementary Table S2). Note that the cell value corresponding to a lipid species and a study group is the sum of the relative intensity of all the matched characteristic ions of the reported lipid species. For example, we have Cho Der-[C 33 H 54 O 6 +H]+ with 100 as the value in all the cells corresponding to HFS, HFW, JS, and JW because this lipid species, graded with Mass ID level, was identified with the base peaks in all the MS/MS spectra. Hence, cell value shall not be used as a means of relative quantitation of the lipid species between the experimental runs.

Variation in Unsaturation between Holstein Friesian Cow and Jaffarabadi Buffalo
The profile of unsaturated lipids obtained from this study followed the same trend that has been reported in HF. For example, C14:1, C16:1, C18:1, C18:2, C22:1 unsaturated lipids were observed in HF during the winter season. The pattern of increasing the unsaturation level in lipids during the winter season was consistent with previous studies of lipid profiling across the season in HF [35]. The DG and TG species were detected across both the seasons in HF, while phospholipids such as PE (e.g., 18:1_18:1) and PC (e.g., 34:1) were exclusively present in HF during the winter season. The observation of lipids exclusively in the winter season could be due to the nutritional availability of greener pasture and feed, which contribute to the availability of lipid precursors [33]. The results of our study showed a similar expression pattern to that observed for saturated and unsaturated fatty acids across both seasons, as reported previously in HF cow milk [36,37]. However, an important factor that has not been adequately controlled and is not clear is the influence of physiological and endocrinological factors on lactation stages and its correlation with the change in lipid composition. The majority of lipids identified in J during winter season were unsaturated lipid species (e.g., C14:1, C18:1, C18:2, C20:4) compared to the summer season. Additionally, phospholipid species like PC 34:1, PE 18:0_18:1 and PE 18:1_18:1 including cholesterol species such as Cho Der-[C 27 H 48 O 6 +Na]1+, Cho Der-20:0, Cho Der-22:0 were detected exclusively in winter, while GlcCer 34:2 was exclusively detected in the summer season. In the present study, higher numbers of saturated fatty acids (SFA) were detected during summer in both animal groups. However, trends in HF were found to be consistent with previous reports of the higher mean value of SFA during summer, and a higher level of unsaturation mean value was detected during the winter season [38], while it has been reported for the first time in J buffalo. The extent of unsaturation plays a vital role in the processing and physical properties of dairy products, for instance, the melting point and spread ability of fat are dependent on the acyl chain length and the degree of unsaturation [39,40].
One of the primary functions of milk lipids is to provide essential nutrients for the development of immune systems in newborns and protect against diseases. For example, a significant number of lipids of Omega family such as oleic acid (18:1), conjugated linoleic acid (18:2) and α-linolenic acid (18:3) previously reported with anti-diabetic and anticarcinogenic effects, arthritis, depression and Alzheimer's disease [41] were confidently identified and constitutively expressed across summer and winter season in both animal groups.
Similarly, the milk fat C18:1 cis-9-to-C15:0 ratio previously known to plays a critical role as a discriminating factor in the diagnosis of hyperketonemia in dairy cows [42] was also constitutively expressed across summer and winter seasons in both groups. Additionally, inflammatory mediators such as oxylipids (e.g., 20-hydroxyeicosatetraenoic acid) derived from polyunsaturated fatty acids, including linoleic acid and arachidonic acid, a known marker of mastitis [43] were confidently identified in our study.

Collating Lipid Ion Species as Sum Compositions and Measuring Ion Abundances for Relative Quantitation
Many of the 1145 MS/MS identified lipid molecular species have the same sum composition or are structural isomers, i.e., lipid species that have same head-group as well as the same total number of carbons and double bonds on the fatty acyl chains. Since these structural isomers have the same observed parent ion mass, they can be grouped into 187 unique sum compositions or Group-IDs. For example, we identified seven structural isomers, namely TG 4:0_6:0_14:0, TG 4:0_8:0_12:0, TG 4:0_4:0_16:0, TG 4:0_10:0_10:0, TG 2:0_6:0_16:0, TG 2:0_4:0_18:0, and TG 6:0_8:0_10:0 from the MS/MS spectra for the lipid molecular composition TG 24:0, we call Group -ID (Supplementary Figure S9).
Of note, many of the lipid structural isomers belonging to a Group-ID or a sum composition comigrated in the LC domain [44,45], thereby posing a significant challenge in deconvoluting these peaks. Hence, we performed relative quantitation at the Group-ID level as done in the literature [32]. For each Group-ID, we also consider the sum of peak areas under-extracted ion chromatograms (XIC) belonging to all the molecular ion species of a Group-ID [46].

Data Normalization and Data Transformation
The abundance of each lipid species was first normalized the measured against the total fat content of the milk. Subsequently, the relative ion abundance was calculated for each lipid molecular species as the ratio of the measured ion abundance of the molecular species to the cumulative sum of ion abundances of all the molecular species belonging to the same class. For example, the relative ion abundance of TG 25:0 species is the ratio of its measured peak area to the cumulative sum of the peak areas of all the 141 TG species. Finally, we use the log2 (log with base 2) of the ratio for analysis and visualization of fold changes describing the n-fold decrease of the ion abundance of an analyte to the total ion abundance of all the species belonging to the same class. Hence, all the fold change values are negative in our study.

Two-Way Analysis of Variance Analysis (ANOVA)
A two-way ANOVA with replication for a "balanced design" was conducted-each subgroup has an equal number of observations-to investigate how the main effects (independent variables), namely breed and season and their interaction, affect the relative abundance of each analyte for the following hypotheses.
Null Hypotheses: bH 0 : There is no difference in the average ion abundances of the analyte between the J and HF milk. sH 0 : There is no difference in the average ion abundance of the analyte between the summer and winter.
iH 0 : There is no interaction effect.
Alternative Hypotheses: bHa: There is a difference in the average ion abundances of the analyte between J and HF. sHa: There is a difference in the average ion abundances of the analyte between summer and winter. iHa: There is an interaction effect between season and breed type on average ion abundances of the analyte.
Detailed proofs for the verification of requisite assumptions of Two-Way ANOVA with Replicates are provided (Appendix A.6).
Finally, Tukey's HSD (honestly significant difference) test was conducted as "post hoc analysis" for analytes that showed significant interaction effect to check if there is significant mean difference between the seasons of each breed (Appendix A.7).

Quantitative Comparison of TG Profiles between Holstein Friesian Cow and Jaffarabadi Buffalo
We summarize the results of the Two-Way ANOVA with replication for balanced design for the 141 TG species as follows: (i) Interaction effect is significant for 91 TG species ( Figure 3A). (ii) Out of the 91 species, 71 TG species have higher seasonal variation-defined as the absolute difference of relative abundances between the summer season and the winter season-for the HF samples than the J samples. On performing Tukey's HSD test to check whether there is a significant variation between the seasons, 47 TG species have statistically significant seasonal variations ( Figure 3B). (iii) Twenty TG species have higher seasonal variation in J samples than in HF samples. However, only four TG species, namely TG 48:4, TG 50:4, TG 51:4, and TG 54:6, vary significantly based on Tukey's HSD Test.
(a) The relative abundances of 81% of the total reported 114 TG species have significant variation between the breeds. Figure 3C showed the 27 TG species that do not have significant variation between the breeds. (b) Relative abundances of 83 TG species, i.e., 59% vary significantly between the seasons. Figure 3D shows the relative abundances of the 58 TG species that do not vary significantly between the season. (c) A total of 53 TG species have not only significant interaction effect but also significant main effects.
Out of these 53 TG species, 48 TG species have higher seasonal variation in HF samples than the J samples of which 40 species are found to be statistically significant ( Figure 3E). Only 5 TG species have higher seasonal variation in J samples than the HF samples. However, Tukey's HSD test showed that none of the TG species is statistically significant.
Metabolites 2020, 10, x 10 of 28 Metabolites 2020, 10, x; doi: www.mdpi.com/journal/metabolites 2.7.3. Data Normalization and Data Transformation The abundance of each lipid species was first normalized the measured against the total fat content of the milk. Subsequently, the relative ion abundance was calculated for each lipid molecular species as the ratio of the measured ion abundance of the molecular species to the cumulative sum of ion abundances of all the molecular species belonging to the same class. For example, the relative ion abundance of TG 25:0 species is the ratio of its measured peak area to the cumulative sum of the peak areas of all the 141 TG species. Finally, we use the log2 (log with base 2) of the ratio for analysis and visualization of fold changes describing the n-fold decrease of the ion abundance of an analyte to the total ion abundance of all the species belonging to the same class. Hence, all the fold change values are negative in our study.

Two-Way Analysis of Variance Analysis (ANOVA)
A two-way ANOVA with replication for a "balanced design" was conducted-each subgroup has an equal number of observations-to investigate how the main effects (independent variables), namely breed and season and their interaction, affect the relative abundance of each analyte for the following hypotheses.
Null hypotheses: bH0: There is no difference in the average ion abundances of the analyte between the J and HF milk.
sH0: There is no difference in the average ion abundance of the analyte between the summer and winter.
iH0: There is no interaction effect. Alternative Hypotheses bHa: There is a difference in the average ion abundances of the analyte between J and HF. sHa: There is a difference in the average ion abundances of the analyte between summer and winter.
iHa: There is an interaction effect between season and breed type on average ion abundances of the analyte.
Detailed proofs for the verification of requisite assumptions of Two-Way ANOVA with Replicates are provided (Appendix A.6).
Finally, Tukey's HSD (honestly significant difference) test was conducted as "post hoc analysis" for analytes that showed significant interaction effect to check if there is significant mean difference between the seasons of each breed (Appendix A.7).

Quantitative Comparison of TG Profiles between Holstein Friesian Cow and Jaffarabadi Buffalo
We summarize the results of the Two-Way ANOVA with replication for balanced design for the 141 TG species as follows: (a) Interaction effect is significant for 91 TG species ( Figure 3A) (A)  (i) Out of the 91 species, 71 TG species have higher seasonal variation-defined as the absolute difference of relative abundances between the summer season and the winter season-for the HF samples than the J samples. On performing Tukey's HSD test to check whether there is a significant variation between the seasons, 47 TG species have statistically significant seasonal variations ( Figure  3B).
(a). The relative abundances of 81% of the total reported 114 TG species have significant variation between the breeds. Figure 3C showed the 27 TG species that do not have significant variation between the breeds.
(b). Relative abundances of 83 TG species, i.e., 59% vary significantly between the seasons. Figure  3D shows the relative abundances of the 58 TG species that do not vary significantly between the season.
(c). A total of 53 TG species have not only significant interaction effect but also significant main effects. Out of these 53 TG species, 48 TG species have higher seasonal variation in HF samples than the J samples of which 40 species are found to be statistically significant ( Figure 3E). Only 5 TG species have higher seasonal variation in J samples than the HF samples. However, Tukey's HSD test showed that none of the TG species is statistically significant.
Some observations from Figure 3E are as follows: i. Different effect of season on different breeds: J samples showed lesser seasonal variation than HF samples. On conducting Tukey's HSD test for Post Hoc analysis of the 53 TG species with Some observations from Figure 3E are as follows: i. Different effect of season on different breeds: J samples showed lesser seasonal variation than HF samples. On conducting Tukey's HSD test for Post Hoc analysis of the 53 TG species with significant interaction effect, only one species, TG 23:0, has significant variation between the summer and winter seasons whereas for the HF samples, 40 TG species showed significant variation between the seasons. ii.
Potential seasonal effect based on the length of the fatty acids: For the HF milk samples, 34 TG species with TC-where TC is total number of carbons in Sn1, Sn2, and Sn3 chains-ranging from 23 to 49 have higher relative abundance in winter season than the summer season while 19 TG species with TC ranging between 51 to 59 have higher relative abundances in summer season than the winter season (rectangle in Figure 3B shows these 19 TG species). Of note, TG 62:2 has higher relative abundance in winter than in summer.
Quantitative comparison of TG lipid profiles of bovine milk [21,32,47] was also performed between HF and J. Of the 78 lipid molecular species (Figure 4) reported in the recently published literature [21,32,34], 45 lipid molecular compositions showed significant variation in their relative ion abundances between the HF and J groups (On performing two-way ANOVA, **: p-value < 0.01, and **: p-value < 0.05, (n = 78, i*,** = 49, b*,** = 61, s*,** = 44, bsi*,** = 29 where n is the total number of TG species, b/s/i/bsi with suffix *,** denotes the number of TG species for which the effect of breed/season/interaction/each main effect as well as the interaction is significant). Each marker in Figure 4 is annotated with lipid composition with suffixes b/s/i followed by */** indicating statistical significance, then followed by reference number of the paper that reported the lipid species.
Metabolites 2020, 10, x; doi: www.mdpi.com/journal/metabolites identified between the breeds. Quantitative comparison of TG lipid profiles of bovine milk [21,32,47] was also performed between HF and J. Of the 78 lipid molecular species (Figure 4) reported in the recently published literature [21,32,34], 45 lipid molecular compositions showed significant variation in their relative ion abundances between the HF and J groups (On performing two-way ANOVA, **: p-value < 0.01, and **: p-value < 0.05, (n = 78, i*,** = 49, b*,** = 61, s*,** = 44, bsi*,** = 29 where n is the total number of TG species, b/s/i/bsi with suffix *,** denotes the number of TG species for which the effect of breed/season/interaction/each main effect as well as the interaction is significant). Each marker in Figure 4 is annotated with lipid composition with suffixes b/s/i followed by */** indicating statistical significance, then followed by reference number of the paper that reported the lipid species. We used A previous study [6] reported the identification of 400 lipids from milk fat, while over 100 triacylglycerols (TG) groups were reported from HF cows [46], however, in our study 57 lipids in HF and 54 lipids in J did not overlap with the previously catalogued lipids, suggesting possible breedspecific differences. Taken together, our results significantly expand the number of identified lipids in the bovine milk and emphasize that it is essential to consider these normal changes before looking for useful lipids for incorporating into dairy foods to benefit human health and markers for disease diagnosis. For instance, odd fatty acids such as 15:0 and 17:0 have been reported as markers for dairy intake as plasma LDL-cholesterol has been related to the dietary intake of saturated fatty acids [48]. A previous study [6] reported the identification of 400 lipids from milk fat, while over 100 triacylglycerols (TG) groups were reported from HF cows [46], however, in our study 57 lipids in HF and 54 lipids in J did not overlap with the previously catalogued lipids, suggesting possible breed-specific differences. Taken together, our results significantly expand the number of identified lipids in the bovine milk and emphasize that it is essential to consider these normal changes before looking for useful lipids for incorporating into dairy foods to benefit human health and markers for disease diagnosis. For instance, odd fatty acids such as 15:0 and 17:0 have been reported as markers for dairy intake as plasma LDL-cholesterol has been related to the dietary intake of saturated fatty acids [48].

Multivariate Analysis of Lipid Profiles between Holstein Friesian Cow and Jaffarabadi Buffalo
Two datasets-(1) Quantitative TG profiles, only the 141 TG species, and (2) quantitative lipid profiles, all the 193 lipid species containing TG, DG, PC, PE, SM, PI, PS, and Chol Der-obtained above were subjected to unsupervised analysis with principal component analysis (PCA) and hierarchical cluster analysis (HCA) looking for natural clustering behavior due to the ion abundances of the lipid molecular species among the samples studied. The analysis was conducted using ClustVis software [49], an open source web tool developed by BIIT (a grouped between the Institute of Computer Science, University of Tartu, Tartu; competence center STACC, and Quretec https://biit.cs.ut.ee/), using the following settings-selected data transformation: none, and data scaling: Pareto scaling (mean-centered and divided by the square root of standard deviation of each variable). Note that, ClustVis software does not display loadings plot, and thus, MS excel chart functionality was used to generate the loadings plot.
The first two principal components-PC1 and PC2-of the dataset 1 accounted for 89% of the explained variance. The scores plot ( Figure 5A) and the corresponding loadings plot ( Figure 5B) [50] show a clear separation between J samples and the HF samples. Furthermore, from the scores plot, we observed that HF samples between the summer and the winter seasons are separated farther apart in both the PC1 and PC2 axes than the J samples indicating higher seasonal variation in the quantitative TG profile of HF samples than the J samples. This observation is in line with our previous finding from Figure 3E that J samples showed lesser seasonal variation than HF samples.
Metabolites 2020, 10, x; doi: www.mdpi.com/journal/metabolites HCA was performed on both the datasets using the following settings: each row, i.e., quantitative profile of a lipid species is centered; Pareto scaling is applied to rows. Each column represents a sample of J/HF collected in a particular season. Both rows and columns are clustered using Euclidean distance and average linkage. Branches of the clustering trees are ordered with highest mean values first for easy visualization of the dendrograms.
HCA of the columns on both the datasets (column dendrogram of Figure 6A,B) shows similar groupings were observed on PCA score plots-JS and JW samples are closer to each other than HFW and HFS samples-indicating lesser seasonal variation on the J samples. Figure 6A,B also shows the row clusters and the heat maps based on datasets 1 and 2 respectively.  PCA analysis of the dataset 2 also produced a similar result in terms of clear separation between the samples based on breed and season (scores plot and loadings plot shown in Figure 5C,D) respectively. However, PCA scores plot based on dataset 1, i.e., TG profiles ( Figure 5A) provides a clearer separation between the J samples collected in summer and winter seasons than the separation between these samples in Figure 5D.
HCA was performed on both the datasets using the following settings: each row, i.e., quantitative profile of a lipid species is centered; Pareto scaling is applied to rows. Each column represents a sample of J/HF collected in a particular season. Both rows and columns are clustered using Euclidean distance and average linkage. Branches of the clustering trees are ordered with highest mean values first for easy visualization of the dendrograms.
HCA of the columns on both the datasets (column dendrogram of Figure 6A,B) shows similar groupings were observed on PCA score plots-JS and JW samples are closer to each other than HFW and HFS samples-indicating lesser seasonal variation on the J samples. Figure 6A,B also shows the row clusters and the heat maps based on datasets 1 and 2 respectively.

Sample Collection
The Director of Research, Junagadh Agricultural University approved (protocol number is 18/1/2016/503) the collection of milk samples. A total of 60 animals, Jaffarabadi buffaloes (n = 30) and Holstein Friesian cows (n = 30) were used to collect milk samples for the study. In summer (April-May), milk samples were collected from 15 Jaffarabadi buffaloes and 15 Holstein Friesian cows post-calving, i.e., the first between 30-40 days while the second between 50-60 days. Similarly, in the winter season (December-January), milk samples were collected from 15 animals from each animal group between 30-40 days and 50-60-days post-calving. The average fat percentage (AFP) for J, and HF were 6.1%, and 3.02%, respectively. The milk fat content was measured by Gerber's method using a butyrometer. Feeding of J buffalo comprised of green fodder, maize (Zea mays), jawar (Sorghum bicolor), dry fodder and concentrate mixture like Amul-dan, Sabar-dan and Banas-dan, a balanced ration consisting of 20-22% protein and 6-7% oil in summer season, while fresh Jinjwa grass (Dichanthium annulatum) was fed in addition to the standard feed in winter season. Feed for Holstein Friesian consisted of doob grass (Cynodon dactylon), concentrate mixture (Amul-dan, Sabar-dan and Banas-dan), maize, jawar in summer season while Berseem (Trifolium alexandrinum) was added to the feed in winter season. The average milk yield of Jaffarabadi buffalo was 18-20 L/day while for HF cows yield is 10-15 L/day.
Raw milk samples were collected and stored at −20 • C and processed for lipid isolation. Thus, for each animal group, two study groups were formed: HFS and HFW for the HF cow, and JS and JW groups for the J buffalo. The suffixes S and W represent summer and winter seasons, respectively.

Isolation of Lipids from Milk Samples
Lipids were extracted from milk, as reported previously [24,51]. In brief, an equal volume of milk (20 µL) from each of the fifteen animals within first lactation stage was pooled (300 µL) and diluted by mixing with 800 µL of MS grade water. The diluted milk was mixed with 4 mL freshly prepared chloroform: methanol (CHCl 3 :MeOH) solution (2:1 ratio, v/v) in a 5 mL glass tube. The mixture was vortexed (Popular Traders, Ambala, India) for 10 min, followed by centrifugation at 2500× g for 10 min at 4 • C (Eppendorf, Hamburg, Germany) to facilitate phase separation. The lower organic phase was carefully transferred to a new glass tube while the upper aqueous phase was re-extracted for lipids by an additional 2 mL of the (CHCl 3 : MeOH) solution. The mixture was further vortexed and centrifuged at 2500× g for 10 min at 4 • C. The organic phase of chloroform/methanol was combined after recovery and dried in a liquid nitrogen stream. Schematic shows the workflow of the study (Supplementary Figure S1).

Tandem Mass Spectrometry Data Analysis for Lipid Identification
Automated annotation of lipid molecular ions using MS and MS/MS data was performed using SimLipid software (PREMIER Biosoft, Palo Alto, CA, USA). The annotated lipid ions were manually reviewed, verified, and graded based on the confidence level of identification achieved from the MS/MS data. We created a list of precursor m/z values of the MS/MS spectra that were not annotated with any lipid ion species and subjected it to "Bulk" Structure Searches-Search COMP_DB with a List of Precursor Ions [52] (for possible annotation of triacylglycerol (TG) ion species. The MS/MS spectra assigned with possible TG ion species were then subjected manual interpretation for possible identification of TG species that have not been catalogued in the LM database.

Relative Quantification of Lipids and Statistical Analyses
Relative quantification of the identified lipid molecular ion species was performed using peak areas of extracted ion chromatograms. The extracted ion chromatograms were generated using SimLipid and Thermo Fisher Scientific's Xcalibur software.
Data normalization, transformation, and statistical analyses pertaining to testing of hypotheses, and generation of charts were performed using Microsoft Excel in-built functions and user-defined functions. We conducted the test of normality of data including the Shapiro-Wilk test using online free tool [53]).

Multivariate Analysis
ClustVis software [49] was used for performing unsupervised multivariate statistical analysis.

Conclusions
In summary, this is the first lipidomic analysis in a healthy cohort of HF cow and J buffalo across summer and winter seasons, resulting in the high confidence identification of 1145 unique lipids species, including 260 triacylglycerols and one diacylglycerol that have not been yet catalogued in LM database. The results of the comprehensive quantitative data analysis show breed and season associated changes not previously reported in bovine milk. Such changes reflect the dynamic nature of the complex milieu of milk and provide a foundation to improve our understanding to evaluate in future studies the effect of additional factors on lipid composition. The results of the present will help to facilitate the elucidation of the biological function of these lipids and their potential usage in infant formula and milk for the betterment of human health.

Supplementary Materials:
The following are available online at http://www.mdpi.com/2218-1989/10/12/507/s1, Figure S1: Schematic of the workflow of the study. Figure S2: Typical SimLipid software graphical user interface (GUI) showing annotations of an MS/MS spectrum with a lipid ion that has been assigned Mass-ID using MS excel in-built functions. Figure S3: Typical SimLipid software GUI showing annotation of an MS/MS spectrum with a lipid ion that has been assigned Group-ID using MS excel in-built functions. Figure S4: Typical SimLipid software GUI showing annotation of an MS/MS spectrum with a lipid ion that has been assigned PA-ID using MS excel in-built functions. Figure S5: Typical SimLipid software GUI showing annotation of an MS/MS spectrum with a lipid ion that has been assigned FA-ID using MS excel in-built functions. Figure S6: Manual annotation of possible fatty acid chains of the [TG 23:0+NH4]+ ion using MS/MS data. Figure S7: XIC of [PS 36:1+H]+ ion and isotopic cluster with a monoisotopic peak at m/z 791.5597. Figure S8 (1): An example MS/MS spectrum of TG 2:0_4:0_18:0 that contains 2:0 as one of the fatty acids. Figure S8 (2): An example MS/MS spectrum of TG 3:0_16:0_18:1 that contains 3:0 as one of the fatty acids. Figure S8    Acknowledgments: Mass spectrometric data were obtained at the Bioanalytical Mass Spectrometry Facility within the Mark Wainwright Analytical Centre of The University of New South Wales. This work was undertaken using infrastructure provided by UNSW, and subsidized access to this facility is gratefully acknowledged. We also appreciate the help of Yashowardhan Verma with the analysis. Thermo Fisher Scientific's Xcalibur '.raw' files were directly read by SimLipid v.6.05 software for extraction of raw data from the LC-MS experimental runs and contained 41,576 spectra (17,264 MS1 spectra, and 24,312 MS/MS spectra). The MS/MS spectra were subjected to SimLipid MS/MS database searching using 5/10 ppm tolerances for the precursor/product ions. The database search was performed in non-targeted mode i.e., no filter based on category, class, and subclass of the lipid species was applied. However, for identification of lipid molecular ion species belonging to glycerolipids, glycerophospholipids, and sphingolipids categories using MS/MS data, the program was set to report only those lipid ions that have at least one diagnostic ion other than the parent ion species e.g., a phosphatidylcholine (PC) lipid molecular ion species will be included in the final result if the diagnostic ion of the head group at m/z 184.073 is observed in the MS/MS spectrum with its peak relative intensity higher than 50%. To avoid isotope overlapping of PC and Sphingomyelin(SM) lipid species, MS/MS spectra with dominant peak at m/z 184.073 with their precursor m/z corresponding to the M+0 peak of the isotopic cluster in the MS1 spectra are only considered for the MS/MS analysis. Similarly, PE lipids were considered only if the MS/MS spectra feature peak at m/z corresponding to neutral loss of 141.0191 from the parent ion. SimLipid database searches annotated 6710 MS/MS spectra with possible lipids ions belonging to subclasses of triacylglycerol (TG), diacylglycerol (DG), phosphatidylethanolamine (PE), cholesterol derivatives (Cho Der), sterols esters, phosphatidylcholine (PC), and Sphingomyelin (SM) from out of the total 24,312 spectra between experimental runs. It must be noted that these annotated lipid ions contain redundant lipid species and need to be verified based on observed characteristic ions in the MS/MS spectra for quality assurance as described in the following sections.

Appendix A.2. Verification and Grading of the Lipid Ion Species Reported by SimLipid Database Search
The SimLipid database contains lipid structures and associated information migrated from LIPID MAPS (LM) database (https://www.lipidmaps.org/). The program follows the same naming convention e.g., TG(12:0/12:0/18:4(6Z,9Z,12Z,15Z))[iso3] bearing the LM ID LMGL03012642, that contains detailed information on the Sn1/Sn2/Sn3 positions, double bond positions in the FA chains and cis-trans isomerism. However, peaks observed in MS/MS spectra may not provide enough data to distinguish these isomers leading to difficulty in reporting the identified lipid ions with confidence. Consequently, as a solution, we have implemented Microsoft (MS) Excel in-built formulas within the SimLipid report-an MS Excel format file containing a list of identified lipid molecular ions along with their corresponding MS/MS spectra list annotated with matched diagnostic ions-to classify the identified lipid molecular ion into four different grades based on: i.
Mass-ID (only mass-resolved lipid molecular ion species); ii.
Group-ID (lipid IDs at sum composition level with only class-specific head group identification); iii.
Detailed explanation on the process of grading lipid IDs based on characteristic ions observed in the MS/MS spectra is described as follows: i.
Mass-ID (only mass-resolved lipid molecular ion species): This is the case when a lipid molecular ion species is assigned to an MS/MS spectrum because its precursor m/z value is within a specified mass tolerance e.g., 5 ppm and MS/MS spectrum features ions corresponding to the parent ion mass e.g., [M+NH 4 ] + , [M+H] + . Since no peak corresponding to any diagnostic ion corresponding to class-specific head groups or FA chains was observed, we use a naming convention that includes the subclass name with the formula of the lipid molecular ion. The MS/MS spectrum (See Supplementary Figure S2) Figure S2). Of note, this assignment was purely on the exact mass database search. There was no diagnostic ion corresponding to the Cholesterol e.g., peak at m/z 369.352. Nevertheless, we still use the Cho Der-because LMST01010173 is the lipid molecular ion species that has the minimum mass deviation (3.5326 ppm) from the observed precursor m/z value in the SimLipid database. ii.
Group ). The diagnostic ions corresponding to the FA chains 14:0 and 18:0 were observed in the spectrum. However, there was no ion corresponding to 13:0. Hence, we assign a graded ID TG 13:0@_14:0_18:0; the suffix @ indicating that the diagnostic ion of this FA chain was not observed in the spectrum. iv.
If the MS/MS spectrum features only the peaks corresponding to diagnostic ions of one of the three FA chains of the TG lipid molecular species e.g., only the peak at m/z 481.4257, then we use the Group-ID, TG 45:0, as the graded ID instead of using the notation TG 13:0@_14:0@_18:0 (See Supplementary Figure S4).  Figure S5).
The SimLipid MS/MS database search result containing 6710 MS/MS spectra annotated with lipid molecular ions were graded using the MS Excel in-built formulas. It must be noted that as an effort to reduce false positives of lipid IDs, we considered only the DG and TG lipid ions that were identified at the FA-ID level. Similarly, only phospholipids that were identified at the Group-ID level were considered. However, for Cho Der lipid species, we considered lipid species annotated at Mass-ID level too because many sterols lipid species generate only parent ion (and neutral loss of H 2 O from it) as fragment ions in their MS/MS spectra (e.g., https://www.lipidmaps.org/data/standards/fetch_gif_ mult.php?&MASS=385&LM_ID=LMST01010066&TRACK_ID=156). The initial MS/MS database search of the raw data using SimLipid software that contains lipid species migrated from the LM database does not annotate many TG lipid species that have been commonly reported in the published literature or some of the unsaturated forms [2,23,34,46,47,55]. Furthermore, the database search results also missed out on the majority of the TG lipid species containing the majority of the reported fatty acyls that have been successfully detected and measured using different separation techniques in the published literature [28,56,57]. Consequently, a stepwise procedure was adopted to manually identify additional TG species.
Step wise procedure for characterization of novel DG and TG lipid species that are not yet cataloged in the LIPIDMAPS database is described as follows.
Step 1. Assignment of TG Composition on Each MS/MS spectrum. A mass list containing all the observed precursor m/z values of all the MS/MS spectra from the raw data files was subjected to "Bulk" Structure Searches-Search COMP_DB with a List of Precursor Ions (http: //www.lipidmaps.org/resources/tools/bulk_structure_searches.php?database=COMP_DB) using the following parameters-Ion adduct: [M+NH4]+; Specify Mass Tolerance (+/− m/z): 0.005 m/z, Specify Chains: All chains, Optionally restrict search by lipid category/class: Tri(acyl|alkyl)glycerols (TG). All the trialkylglycerols are removed from the results.
Step  Figure S6). (iii) Finally, TG 23:0 could have the possible lipid species TG 4:0_4:0_15:0 (we use the underscore, "_", between the fatty acid chains to indicate non-specificity of the Sn1/Sn2/Sn3 positions of these reported fatty acid chains). (iv) Once we establish the probable lipid species for an MS/MS spectrum, we create structures using LM and create a glycerolipid structure from LM Abbreviation". (v) (http://www.lipidmaps.org/tools/structuredrawing/StrDraw.pl?Mode=SetupGLStrDraw) wherein we enter the lipid species abbreviation and selecting the option of "Allow arbitrary chain abbreviation" as YES. (vi) The generated SDF files are imported into the SimLipid server database for automated assignment of lipid species on the MS/MS spectra across all the raw data files. (vii) The process (i)-(v) is repeated for all the MS/MS spectra with annotated TG compositions. The p-value of the mean-based F = 0.6413 with degrees of freedoms 3, and 776 is 0.5886. This accepts the null hypothesis that the all the variances of all the four groups are equal.