2.1. Preparation of Milk Samples

Samples of donkey milk were collected from Xinjiang Yukunlun Natural Food Engineering Co. Ltd. in Yopurga County of Xinjiang Province (XJ) in Western China, and from Shandong Dong’e ejiao Company Limited in Dong’e County of Shandong Province (SD) in Eastern China, respectively. The breed of XJ donkey was Jiangyue donkey, a hybrid breed by crossing local donkeys with the introduced Guanzhong donkey from Shanxi province, while SD donkey was Dezhou donkey, a breed produced in the counties along the Bohai Sea in the northern Shandong plain. Samples were collected one batch per day from a collected mixture of donkey milk from 20 donkeys, and within five consecutive days, a total 10 batches of donkey milk were collected from the two plants. All samples were collected from bulk milk containers into sterile tubes, transported to the laboratory immediately after sampling, and then stored at −80°C until further analysis.

2.2. Protein Extraction, and Digestion TMT Labeling

Skimmed milk was prepared and then centrifuged at 10,000g for 30 min to obtain a supernatant of whey proteins, which was consequently determined by the bicinchoninic acid method and sodium dodecyl sulphate-polyacrylamide gel electrophoresis (SDS-PAGE) electrophoresis [20,21]. Protein lysate (8M urea, 1% SDS, containing protease inhibitor) was used for sample lysis [22].

Protein of 200 μg was supplemented to the volume of 180 μl with lysate. The samples were mixed with Tris(2-carboxyethyl)phosphine (final concentration was 10 mM) and were then incubated at 37°C for 60 min, followed by a reaction with appropriate iodoacetamide (final concentration was 40 mM) at room temperature for 40 min in the dark. Each tube of sample was added by precooled acetone, and incubated at −20°C for 4 h before centrifugation at 10,000g for 20 min. Trypsin digestion was performed at a substrate ratio of 1:50 (w/w), followed by an overnight incubation of enzymatic hydrolysis at 37°C [22,23].

The digested peptides were labeled with a TMT Reagent Kit (Thermo Fisher Scientific, Wilmington, USA, No. 90111) according to the manufacturer’s protocol. Samples after tagging were prepared for further analysis.

2.3. One-dimensional Separation of UPLC with High pH

Sample were fractionated using high pH reverse phase separation. The TMT-tagged peptides were re-suspended with ultra-performance liquid chromatography (UPLC) loading buffer, and separated by an ACQUITY UPLC BEH C18 Column (1.7 µm, 2.1 × 150 mm², Waters Corporation, Milford, USA) in high PH liquid phase. Separation was performed at 200 µL/min with the gradient increased for 76 min (Phase A: 2% acetonitrile, pH 10; Phase B: 80% acetonitrile, pH 10), and the ultraviolet detection wavelength was 214 nm.

2.4. LC-MS/MS Analysis

Experiments were performed on a Q-Exactive mass spectrometer (Thermo Fisher Scientific) that was coupled with Easy-nLC 1200 (Waters Corporation), and Thermo Xcalibur 4.0 (Thermo Fisher Scientific) was used as the data acquisition software. A 4-μL of each fraction was injected for nano liquid chromatography with tandem mass spectrometry (LC-MS/MS) analysis. The peptide mixture (2 μg) was loaded onto a C18-reversed phase column (75 μm × 25 cm, Thermo Fisher Scientific) in buffer A (2% acetonitrile and 0.1% Formic acid) and separated with a linear gradient of buffer B (80% acetonitrile and 0.1% Formic acid) at a flow rate of 300 nL/min for 66 min. Q-Exactive mass spectrometer was operated in the data-dependent mode to switch automatically between MS and MS/MS acquisition. Survey full-scan MS spectra (350–1300 m/z) were acquired with a mass resolution of 70 K, followed by 20 sequential high energy collisional dissociation MS/MS scans with a resolution of 35 K. In all cases, one microscan was recorded using dynamic exclusion of 18 s.

2.5. Sequence Database Searching for Protein Annotation

The software version for library search is Proteome Discoverer Software 2.1. Acquired MS/MS spectra were analyzed for protein identification and relative quantification using Proteome Discoverer (Thermo Fisher Scientific, version 2.1) against the protein database from NCBInr (version 29/08/2018; donkey), with a precursor mass tolerance of 20 ppm and fragment mass tolerance of 0.02 Da. The highest score for a given peptide mass (best match to that predicted in the database) was used to identify parent proteins. Tryptic digestion was specified with up to two missed cleavages. Peptide spectral matches were validated based on q-values at a 1% false discovery rate.

2.6. Statistics

The t-test function in R language was used to calculate the p-value of the significant difference between samples. In this project, the screening standard of the significant differentially expressed proteins was p < 0.05 & (fold change < 0.83 or FC > 1.20). Principal Component Analysis (PCA) was carried out against proteomic compositions to observe differences between XJ and SD. Gene Ontology (GO) analysis of identified proteins was carried out by BLAST2GO software (version 2.5.0). Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis was conducted by KOBAS (version 2.1.1) according to the KEGG database (http://www.genome.jp/kegg/pathway.html).

3.1. Casein Proteins Profiles in Donkey Milk

Total 15 casein proteins were identified from the two sample groups (Table 1). Among them, three proteins including к-casein, α_s2-casein and casein kinase were overexpressed in the SD samples (3.2, 1.4, and 1.3 times of that in the XJ samples, respectively), while only one protein predicted к-casein was found to be overexpressing in the XJ samples. Proteins including α_s1/α_s2–casein and precursors, casein kinases and substrates, β-casein, and caseinolytic peptidase were expressed in both sample groups (Table 1). Our study revealed a characteristic caseins pattern in the two donkey samples consisting of: (i) three α_s1-CN components; (ii) three α_s2-CN components; (iii) one β-CN components; (iv) two к-CN components; and (v) four CN kinase and one caseinolytic peptidase (Table 1).

3.2. Overview of Whey Proteins and Functions in Donkey Milk

The overall whey proteins and functions were analyzed by GO and KEGG analysis, respectively (Figure 1). GO analysis was performed to identify the Biological Process (BP), Cellular Components (CC), and Molecular Functions (MF) (Figure 1A), while analysis of KEGG pathways was performed to identify the functions including six of Brite A class and 42 of Brite B class (Figure 1B). Among all of the identified whey proteins, total 46 terms of the three GO classes were clustered, with 10 different molecular functions, 15 kinds of cellular components, and 21 diverse biological processes. To be noticed, the summation of the numbers of proteins presented within the 46 terms was 1982, which was much higher than the actual identified number of proteins 620, indicating that some whey proteins performed several functions and thus belonged to multiple GO classes. For the KEGG analysis, 184, 141, 117, 87, 84, and 35 proteins were related to the six Brite A functions of organismal systems, human diseases, cellular processes, metabolism, environmental information processing and genetic information processing, respectively. The top 10 KEGG Brite B functions included immune system, transport and catabolism, global and overview metabolic pathway maps, signal transduction, cancers, endocrine system, infectious viral diseases, infectious bacterial diseases, digestive system, and cell growth and death. Further phylogenetic classification for orthologous proteins were analyzed by annotation based on Cluster of Orthologous Groups of proteins. Among the total 25 functions in the database, 23 functions were annotated (Supplementary Figure S1), with the top five functions including posttranslational modification, protein turnover, chaperones (O), General function prediction only (R), Signal transduction mechanisms (T), Carbohydrate transport and metabolism (G), and Amino acid transport and metabolism (E).

Figure 1

Functional analysis based on gene ontology level 2: Biological Process (BP), Cellular Components (CC), and Molecular Functions (MF) (A), and Kyoto Encyclopedia of Genes and Genomes (B) of whey proteins in donkey milk from China.

3.3. Common and Different Whey Proteins in Donkey Milk from the Two Groups

The comparison of proteome data of whey proteins between the XJ and SD samples were initially analyzed by PCA and correlation heatmap of all proteins detected from the two groups (Figure 2). The samples from different groups clustered separately, while the samples from the same group clustered together, suggesting that there were great differences in expression of whey proteins between the XJ and SD samples, and the samples within each group had good repeatability.

Figure 2

Correlations of expression of proteins between the Xinjiang Province (XJ) and Shandong Province (SD) samples of donkey milk. (A) Principle component analysis of proteomic data. (B) Correlation heatmap of whey proteins between the two sample groups.

Profiles of whey proteins in donkey milk were compared between the XJ and SD samples. A total of 620 proteins were identified and relatively quantified by proteomic analysis, including 438 of the common proteins in both sample groups (Figure 3A and 3B). There are 26 whey proteins with a statistically significant greater than twofold change in the samples of donkey milk from XJ and SD. Compared with those in the SD samples, 20 proteins such as Whey-acidic-protein (WAP) four-disulfide core domain protein and elafin were up-regulated, and six proteins including GDP-mannose 4,6-dehydratase and protein S100-G were down-regulated in the XJ samples (Figure 3A and Supplementary Table S1). The common proteins with the similar abundance (p > 0.05) in the two sample groups occupied approximately 70.6% of the total identified whey proteins, which indicates that whey proteins in the two groups of donkey milk shared a big part of fractions of whey proteins. However, 105 and 77 proteins were identified to overexpress in XJ and SD samples, respectively (p < 0.05), suggesting that each group had unique fraction of proteins profile (Figure 3B). The hierarchical clustering analysis of differentially expressed whey proteins further revealed two different clusters corresponding to the two sample groups, with total 182 proteins significantly more abundant in either of the XJ or in the SD samples (Supplementary Figure S2).

Figure 3

Common and different expression of whey proteins in donkey milk from Xinjiang Province (XJ) and Shandong Province (SD). (A) Expression of whey proteins. Shown is the log₂ transformed ratio of the protein abundance in the XJ samples relative to the abundance in the SD samples. Grey dots: background proteome. The level of significance was evaluated based on p < 0.05 and the fold change two times cut-off. Red dots: proteins expressed uniquely high in the XJ samples; blue dots: proteins expressed uniquely high in the SD samples. (B) Venn diagram of number of whey proteins expressed in the two sample groups. The fold change 1.2 times cut-off.

Functions of biological process, cellular components, and molecular functions associated with the differentially expressed proteins were identified based on the analysis of GO level 2 (Figure 4). Functional categorization showed that the numbers of differently expressed proteins in each of the GO term were all bigger than that of the SD samples, indicating that more proteins involved in the three targeted functions in the XJ samples rather than the SD samples. For example, 31 proteins of cellular process related to biological process, 35 extracellular region part proteins of the cellular components, and 33 binding proteins related to molecular functions were overexpressing in the XJ samples, respectively, while only another nine, nine and 18 proteins of the same GO functions were overexpressing in the SD samples, respectively (Figure 4).

Figure 4

Functional categorization based on Gene Ontology (GO) level 2 analysis of differentially abundant whey proteins in donkey milk from Xinjiang Province (XJ) and Shandong Province (SD).

In addition to a single protein and its function, identification of biological network is very important for understanding how proteins interact with each other and function within the cell. The interactions of differentially expressed whey proteins in donkey milk were investigated using the Search Tool for the Retrieval of Interacting Genes/Proteins database (Supplementary Figure S3). The combined score of the interaction networks ranged from 150 to 999, and the higher score reflects the higher probability that the interaction is supported by several types of evidences and biologically meaningful. Total 112 proteins were enriched in networks. Important proteins overexpressed in the XJ samples with a high fold change included tissue-type plasminogen activator (big and red node, accession no.: XP_014693260.1), which involved in protein–protein interactions with 74 proteins. Neuroserpin and apolipoprotein B-100 (big and green nodes, accession no.: XP_014701927.1 and XP_014721916.1, respectively) were highly expressed in the SD sample rather than in the XJ samples, which had interactions with 58 and 68 proteins, respectively.

4.1. Favorable Properties Indicated by a Low Casein/Whey Protein Ratio

Caseins were characterized with allergenic properties, and a low casein/whey protein ratio made it more favorable for human nutrition [2–4,6,24]. The casein/whey ratio in our study was similar between the samples from XJ and SD (p > 0.05), with a low level of 0.48 and 0.43 (Supplementary Figure S4), respectively, both much lower than the previously reported ratio of 1.3 in donkey milk collected from the same area of XJ, and also lower than the mean ratio of 0.67 in donkey milk from Italy [4,9]. The low level of casein/whey ratio is much more lower than other milk (such as cow, sheep, goat and mare), very close to the value of human milk (0.4), and thus plays an important role in reducing its allergenicity [2,3,6]. The lower ratio of casein/whey relating to the lower ratio of soluble calcium also provides the special property of soft curd formation for donkey milk, and thus makes donkey milk a good substitute of human milk for bottle-fed infants [2,3,6]. Therefore, the products studied in our study are suggested to be very good sources of human nutrition with a favorable casein/whey ratio.

4.2. Specificity of Important Caseins and New Isoforms

In most milk types, four gene products including α_s1-, α_s2-, β-, and к-caseins are synthesized in mammary gland in response to lactogenesis hormones and other stimuli, with α_s1-casein as the most important protein causing milk allergy [9]. Within the previously identified isoforms of α_s1-casein, only one reported type with 202 amino acids was observed in this study. Actually, numbers of amino acids of α_s1-casein were found to be diverse in different samples, for example, four α_s1-casein isoforms with 202, 201, 197 and 196 amino acids were detected in a dephosphorylated milk sample collected from an individual donkey at middle lactation stage, while a higher molecular mass with 210 amino acids long was observed in donkey milk in another research [9,25]. To be noticed, there is one new isoform of α_s1-casein consisting of 208 amino acids in our study, which has not been reported in donkey milk before. The existence in an individual milk sample of a full-length α_s2-casein with 221 amino acids and its four different isoforms of 216, 185, 214, and 183 amino acids have been reported in previous studies [3,8,26]. However, the two α_s2-casein isoforms identified in our study had 236 and 231 amino acids, which were firstly observed in donkey milk. Actually, α_s2-casein is a very minor component in donkey milk, thus the absence of these two α_s2-casein components in previous studies were very likely due to the very low amount of them in donkey milk [2,3,24]. The primary structure of two β-casein isoforms as the most abundant casein component, with 226 and 218 amino acids, respectively, have been characterized in donkey milk in literatures [25,27]. Unlikely in our study, one new isoform of β-casein with 233 amino acids was observed in both of the samples with the similar abundance. As the amount of к-casein in the fraction of Equidae is also very low, two components of к-casein both with 182 amino acids were identified in our study, indicating that our methods had a very high accuracy [3]. Moreover, the discovery of the two к-casein components in our study was the first time on account of that only one isoform of к-casein with 162 amino acids long has been reported in donkey milk previously [3]. Thus among the total seven isoforms of caseins, only one has been reported in previous studies, and the other six of them were firstly identified here in donkey milk.

4.3. Variations of Major Whey Proteins and Functions

As the major constituents of whey proteins, β-lactoglobulins and α-lactalbumin were found to have similar expression profile with fold change of no more than 1.2 (Table 2). Lysozyme was recognized as the third major whey protein in donkey milk, and much more abundant than that in other milk such as cow, goat, sheep and even human milk [2–4]. In our study, donkey milk from SD overexpressed lysozyme in comparison with that of XJ while expressed less lysozyme precursor, which would be turned into more lysozyme after post-translational modification. Even though most gram-positive bacteria can be initiated lysis by lysozyme, some coccus-shaped species are resistant to lysozyme due to the O-acetylation of its peptidoglycan by O-acetyltransferase [1,6,28]. Interestingly, our previous work found that gram-positive bacteria Macrococcus and Streptococcus were identified as the top predominant genera in the same samples of donkey milk (data not published), which is consistent with the observation of lysozyme in this study because of their resistance to lysozyme. Moreover, absence of Listeria, Salmonella, or Shigella in the same samples of donkey milk (data not published), is very likely related to antimicrobial effects of at the presence of lysozyme and lactoperoxidase in it (Table 2) [18,29]. Immunoglobulins (mainly IgG, IgA, IgM) and lactoferrin was reported to occupy about 11.5% and 4.5% of the total whey proteins in donkey milk, however, our proteomic analysis has no observation of these two enzymes but only five immunoglobulins-related components as shown in Table 2 [4,5]. Among them, IgGFc-binding protein has been suggested to be important in the immunological defense for cell protection and anti-inflammation in tissues based on its immunoglobulin G-fragment crystallizable (IgG-Fc) binding property and tissue distribution [30,31]. The J chain is an important polypeptide which regulates polymer formation of IgA and IgM, and its incorporation into polymeric IgA and pentameric IgM endows them with several special features including agglutination of bacteria, non-inflammatory reaction with limited complement-activating activity, and a high affinity for the polymeric immunoglobulin receptor of transporting them into milk [32]. In addition, serum albumin was reported to account for 6.2% of total whey proteins in donkey milk previously, however, only its precursor was determined in this study, the reason of which is unclear [2,4,5]. Noticeably, serum albumin is not synthesized in the mammary gland but may enter the milk by uptake with immunoglobulins, which may partially explain the absence of serum albumin and immunoglobulins, but further studies are necessary to do for more precise explanations [3].

4.4. Potential Factors for Proteins Differentiation between the Two Provinces

Differentiation on expression of whey proteins between the XJ and SD maybe ascribed to several factors including donkey breed, feed, climate, and lactation period [33,34]. The detailed information of the two sample groups have been summarized and shown in Supplementary Table S2. The most dominant protein WAP Four-Disulfide Core domain protein 2 (WFDC2) in XJ belongs to a large family of WFDC proteins with multiple functions including antiprotease, anti-bacterial, anti-viral and anti-inflammatory properties [35–37]. It has been shown that the expression of WFDC2 was related to lactation period, for example, the expression is elevated in early lactation, down regulated in mid lactation, subsequently increase slowly and markedly during involution [38]. Our samples were all collected from jennies with mature milk, however, specific lactation period of them might be various, which reasonably induced the significant difference on the expression of it between the two groups. Moreover, different expression of some proteins may be also related to the health status of donkey breast, for example, the higher expression proteins in XJ including clusterin, C-reactive protein and serum amyloid A protein-like have been proposed to up-regulate in milk from donkey bodies with mastitis [39]. Component of feed is another important factor which may affect the expression of some proteins, for example, protein S100-G, the expression of which was much higher in SD rather than in XJ, has been indicated to have multiple binding partners and functions and utilize calcium to propagate a response through a signaling pathway as a calcium-binding protein [40]. The up regulation of this protein in SD may be related to the addition of vitamin D3 and calcium ions in feed (Supplementary Table S2) [41]. It should be noticed that the differentiation of these proteins may not only result from single factor but simultaneous effects of multiple factors, thus further studies of concrete analysis on combined factors are necessary to be performed in the future.