Exploration of the lactation function of protein phosphorylation sites in goat mammary tissues by phosphoproteome analysis

Protein phosphorylation plays an important role in lactation. Differentially modified phosphorylation sites and phosphorylated proteins between peak lactation (PL, 90 days postpartum) and late lactation (LL, 280 days postpartum) were investigated using an integrated approach, namely, liquid chromatography with tandem mass spectrometry (LC-MS/MS) and tandem mass tag (TMT) labeling, to determine the molecular changes in the mammary tissues during the different stages of goat lactation. A total of 1,938 (1,111 upregulated, 827 downregulated) differentially modified phosphorylation sites of 1,172 proteins were identified (P values < 0.05 and fold change of phosphorylation ratios > 1.5). Multiple phosphorylation sites of FASN, ACACA, mTOR, PRKAA, IRS1, RPS6KB, EIF4EBP1, JUN, and TSC2 were different in PL compared with LL. In addition, the Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis showed that the calcium signaling pathway, oxytocin signaling pathway and MAPK signaling pathway were enriched. The western blot results showed that the phosphorylation levels of ACACA (Ser80), EIF4EBP1 (Thr46) and IRS1 (Ser312) increased and JUN (Ser63) decreased in PL compared with LL. These results were consistent with the phosphoproteome results. In this study, we identified for the first time the differentially modified phosphorylation sites in goat mammary tissues between PL and LL. These results indicate that the multiple differentially modified phosphorylation sites of FASN, ACACA, mTOR, PRKAA, IRS1, RPS6KB, EIF4EBP1, TSC2, and JUN and proteins involved in the calcium signaling pathway, oxytocin signaling pathway, and MAPK signaling pathway are worthy of further exploration.


Background
Phosphorylation is one of the important posttranslational modifications of proteins; it is related to many activities of life, such as signal transduction, gene expression, cell cycle and cell apoptosis [1]. A phosphoric acid group with a strong negative charge is added to an amino acid of the protein, thereby changing its structure, activity and ability to interact with other molecules, affecting many biological processes. In eukaryotes, phosphorylation occurs primarily on serine, threonine and tyrosine [2].
Goat milk is becoming very popular with the improvement of people's living conditions and their understanding of the unique nutritional and health value of goat milk. Therefore, improving the milk production and quality of dairy goats is the focus of our work. Dynamic changes in protein phosphorylation occur in goat mammary epithelial cells during the PL period and LL period [3,4]. Previous studies have shown that protein phosphorylation regulates lactation in animals [5,6]. Although phosphoproteome analysis was conducted on the oocytes of goats [7] and the heart and liver of mice [8,9], phosphorylated proteomic analysis has not been performed in the PL and LL stages of dairy goat mammary tissues. The regulatory mechanism of phosphorylated proteins during lactation is not well known.
Therefore, in this experiment, differentially modified phosphorylation sites of phosphorylated proteins between PL and LL in goat mammary tissues were identified by tandem mass tag (TMT)/isobaric tags for relative and absolute quantitation (iTRAQ) labeling, highperformance liquid chromatography (HPLC) fractionation, affinity enrichment and liquid chromatography with tandem mass spectrometry (LC-MS/MS) analysis. A total of 1,938 (1,111 upregulated, 827 downregulated) differentially modified phosphorylation sites of 1,172 proteins were identified and analyzed with the Gene Ontology (GO) classification system and Kyoto Encyclopedia of Genes (KEGG) for differentially phosphorylated proteins.
The phosphorylation levels of ACACA (Ser80), EIF4EBP1 (Thr46), and IRS1 (Ser312) were increased, and that of JUN (Ser63) was decreased in PL compared with LL, and the data showed that the phosphoproteome results were reliable. The present study provides essential information to enhance our knowledge of the expression dynamics of phosphorylated proteins and phosphorylated sites during PL and LL, and it provides reference data for further studying the roles of protein phosphorylation during lactation in goat mammary tissues.

Mammary tissue collection
Guanzhong dairy goats were obtained from Longxian Goat Breeding Center, in Longxian County of Shaanxi Province, China. All procedures in our animal study were approved by the Animal Care and Use Committee of the Northwest A&F University (Yangling, China) (permit number: 17-347, data:2017-10-13). Mammary tissues were collected during the LL (280 days postpartum) and PL (90 days postpartum) periods from 18 healthy Guanzhong dairy goats (3 year olds) by surgery (9 mammary tissues in PL and 9 mammary tissues in LL). First, one side of the udder was wiped clean with alcohol, and then local anesthesia was applied to the tissue (Xilagesic 20 %, Calier, Barcelona, Spain). A 5-6 cm incision was made and the skin was pulled back using sterile forceps, exposing the mammary tissue. Samples (1 cm 2 ) were taken using a sterile scalpel blade and forceps. Pressure was applied with sterile gauze to stop any bleeding. The incision was closed with 6 to 8 surgical staples (#89,063,337, Appose ULC Skin Stapler, 35 wide; Henry Schein Inc., Melville, NY). After that, the goat was kept in a clean place during the healing period, and the skin around the wound was sterilized with iodine tincture and 75 % alcoholregularly. All collected tissues were washed with RNase-free PBS, quickly frozen in liquid nitrogen and stored at -80°C.

Protein extraction
The mammary tissues epithelial tissues of dairy goats in PL and LL were grinded into cell powder in the presence of liquid nitrogen. The cell powder was transferred to a 5 mL centrifuge tube. Subsequently, the cell powder was added to four volumes of lysis buffer (8 M urea, 1 % Protease Inhibitor Cocktail, 2 nM EDTA), followed by thrice of ultrasonic treatment on ice with a high-intensity ultrasonic processor (Scientz, Ningbo, China). The remaining debris was removed by centrifugation at 12,000 g at 4°C for 10 min. Finally, the supernatant was collected, and the protein concentration was quantified using the BCA kit (Scientz, Ningbo, China). Proteins were extracted from 18 mammary tissues of dairy goats that included nine mammary tissues in PL and nine mammary tissues in LL. Afterwards, proteins extracted from nine mammary tissues in PL were divided into three groups. Each group was made up of three parts of protein extracted from PL mammary tissues. The same method was used to construct the LL protein library.

Trypsin digestion
The protein solution was reduced with 5 mM dithiothreitol for 30 min at 56°C and alkylated with 11 mM iodoacetamide for 15 min at room temperature in the dark. Consequently, the protein solution was digested. The protein sample was diluted by adding 100 mM triethylammonium bicarbonate (TEAB). The concentration of urea in the solution should be less than 2 M. Finally, trypsin and protein were added at 1:50 mass ratio for the first digestion overnight and at 1:100 mass ratio for the second digestion for 4 h.

TMT labelling and HPLC fractionation
The production was desalted by a Strata X C18 SPE column (Phenomenex, Torrance, USA) after trypsin digestion and vacuum drying. Afterwards, peptides were reconstituted in 0.5 M TEAB and processed using a TMT kit (Scientz, Ningbo, China). Peptides were labelled using TMT for quantitation; labelling was performed as previously described [10]. One unit of TMT reagent was thawed and reconstituted in acetonitrile. The peptide mixtures were incubated for 2 h at room temperature, pooled, desalted and dried by vacuum centrifugation. The tryptic peptides were divided into segments using a Thermo Betasil C18 column (5 μm particles, 10 mm ID, 250 mm length) by high pH reversed-phase HPLC. In summary, peptides were separated into 60 fractions after more than 60 min of gradient treatment with 8-32 % acetonitrile (pH 9.0). Subsequently, the peptides were combined into eight fractions and dried by vacuum centrifuging. Reactions were quenched using 8 µL of 11 % lysine following the manufacturer's recommendations. Finally, 10 % trifluoroacetic acid (TFA, Sigma-Aldrich) was added into the peptide solutions. The volume ratio of the peptide solution to TFA was 1:10, and the mixture was pooled, desalted and dried by vacuum centrifugation.

Affinity enrichment
The mixtures were first incubated with immobilised metal affinity chromatography (IMAC) microsphere suspensions with vibration in loading buffer (50 % acetonitrile/6 % TFA). Phosphopeptide-rich IMAC microspheres were collected by centrifugation. Afterwards, the supernatant was removed. The IMAC microspheres were washed with 50 % acetonitrile/6 % TFA and 30 % acetonitrile/0.1 % TFA to remove nonspecifically adsorbed peptides. To elute the enriched phosphopeptides from the IMAC microspheres, an elution buffer containing 10 % NH 4 OH was added, and the enriched phosphopeptides were eluted with vibration. The supernatant containing phosphopeptides was collected for LC-MS/MS analysis after lyophilisation.

LC-MS/MS analysis
The tryptic peptides were dissolved in 0.1 % formic acid (solvent A), directly loaded onto a homemade reversedphase analytical column (15 cm length, 75 μm i.d.). The gradient consisted of solvent B (0.1 % formic acid in 98 % acetonitrile), which increased from 6 to 23 % over 26 min, 23-35 % in 8 min and climbing to 80 % in 3 min, followed by holding at 80 % for the last 3 min, all at a constant flow rate of 400 nL/min on an EASY-Nlc 1,000 ultra-performance liquid chromatography (UPLC) system.
Peptides were characterised using tandem mass spectrometry (MS/MS) in Q Exactive™ Plus (Thermo) coupled online to the UPLC after the peptides were subjected to NSI source. The electrospray voltage applied was 2.0 Kv. The m/z scan was 350 to 1,800 for full scan, and intact peptides were detected in the Orbitrap at a resolution of 70,000. Using a normalization collision energy setting of 28, the fragments were detected in the Orbitrap at a resolution of 17,500 for the selected MS/ MS of peptides. A data-dependent procedure was alternated between one MS scan followed by 20 MS/MS scans with a 15.0 s dynamic exclusion. The automatic gain control was set at 5E4. The fixed first mass was set at 100 m/z.

Database search
Maxquant search engine (V.1.5.2.8) was used to process the resulting MS/MS data. The UniProt Ovis aries (27,472 sequences) concatenated with a decoy database that was used for the tandem mass spectra search. The search parameters were as follows: 6-plex TMT, the minimum length of the peptides was set at 7 amino acid residues and the maximum modification number of the peptides was set at 5; fixed modification (carbamidomethylation of cysteine residues); variable modification (protein N-term acetylation); and oxidation of methionine residues and phosphorylation of serine, threonine and tyrosine residues, thereby allowing two missing cleavages by cleavage enzyme trypsin/P. The mass error was set at 20 ppm in the first search and 5 ppm in the primary search, and the fragment ions of mass tolerance were set at 0.02 Da. The minimum score for modified peptides was set at > 40, and a false discovery rate was set at < 1 %.

Bioinformatics analysis
Gene Ontology (GO) annotation data were obtained from UniProt GoA database. Proteins were divided into three categories by GO annotation, namely, a biological process, cellular compartment and molecular function. The GO with a corrected P value < 0.05 was considered significant. InterProScan (http://www.ebi.ac.uk/interpro/ ) was used to detect the enrichment of functional domains of differentially expressed proteins compared with all identified proteins. Protein domains with a corrected P value < 0.05 were considered significant. The KEGG database (https://www.genome.jp/kegg/tool/map_ pathway2.html) was used to identify enriched pathways [11,12]. The pathway with a corrected P value < 0.05 was considered significant. Based on the KEGG website, these pathways were classified into hierarchical categories. Fisher's exact test was used to examine the enrichment of the differentially phosphorylation protein against all identified phosphorylation proteins. Wolfpsort, a subcellular localisation predication wolf, was used to predict subcellular localisation. Soft motif-x was used to analyse the sequence model constituted with amino acids in specific positions of modify-21-mers (10 amino acids upstream and downstream of the site) in all protein sequences. In addition, all database protein sequences were used as background database parameter, as well as other parameters with default.

Enrichment-based clustering
Hierarchical clustering was based on the functional classification of different proteins. Categories that were at least enriched in one of the clusters with P < 0.05 were filtered after collating all the categories and their P values. This filtered P value matrix was transformed by the function x = − log10 (P value). Finally, these x values were z-transformed for each functional category. These z scores were then clustered by one-way hierarchical clustering (Euclidean distance and average linkage clustering) in Genesis. Cluster membership was visualised by a heat map using the "heatmap.2" function from the "gplots" R-package.

Western blot
Approximately 25 µg of protein extracted from the goat mammary gland epithelial tissues was separated by SDS-PAGE and then transferred onto polyvinylidene difluoride membranes (PVDF, Merck Millipore, MA, USA). The membranes were immersed in 10 % nonfat milk powder in Tris-buffered saline containing 0.1 % Tween 20 (pH 7.6) for 2.5 h at room temperature. The membrane was incubated with the corresponding primary antibodies overnight at 4°C (Table 1). In the next step, the membrane was incubated with the suitable HRPconjugated secondary antibodies against mouse and rabbit at 4°C for 2 h. The proteins were visualized using ECL prime western blotting detection reagent (Amersham, GE Healthcare Lifesciences, Sweden) by gel documentation system (Biospectrum 410, UVP). Proteins were quantified by the Quantity One program (Bio-Rad, California, USA).

Statistics
In this study, the reporter ion intensity of the modified peptide was used to calculate the modification site quantitation. First, the average quantitation of each sample was calculated in multiple replicates, and then the ratio of the average quantitation between the two samples was calculated. The ratio is used as the relative quantitation between the two samples. The next step is to calculate the significant p value of differential modification between the two samples. First, the relative quantitative values of each sample are taken as log2 transform (so that the data conformed to the normal distribution), and the p value was calculated by the two-sample two-tailed T-test method. A p value < 0.05 and a modification ratio > fold change were regarded as upregulation. A p value < 0.05 and protein ratio < 1/fold change was regarded as downregulation. Each experiment was repeated at least three times. All of the data were processed by SPSS 19.0 (Beijing, China). The western blot results are shown as the means ± SE (standard error), and the differences were compared by the two-sample two-tailed T-test method (**P < 0.01, *P < 0.05).

Sample repeatability test and mass spectrometry quality control
Differentially modified phosphorylation sites and phosphorylated proteins from the PL and LL of goats were analyzed by LC-MS/MS and TMT. The workflow of the present study is shown in Fig. 1a. The results of the sample repeatability test are shown in Fig. 1b. Figure 1b is a heatmap drawn by calculating Pearson's correlation coefficient between two pairs of samples. Pearson's correlation coefficient showed that the reproducibility of the sample was good. The mass error was detected at 0 as the center axis and was concentrated in a range of less than 10 ppm. The length of most peptide segments was between 8 and 20 amino acid residues, which indicated that the preparation of the sample was successful in our study ( Figure S1).

Analysis of differentially phosphorylated proteins
A quantitative study of differentially modified phosphorylation sites was performed in this experiment by an integrated approach involving LC-MS/MS and TMT labeling. All samples were prepared in triplicate. There were 6,979 phosphorylated sites distributed on 2,608 proteins. In addition, 5,901 phosphorylated sites distributed on 2,454 proteins were quantified. In the classified groups, the variation of the difference expressed was more than 1.5 times and less than 0.67 times the change criterion of significant upregulation or significant downregulated (P values < 0.05). A total of 1,938 (1,111 up-regulated and 827 down-regulated) differentially modified phosphorylation sites of 1,172 proteins were changed ( Fig. 2; Table S1). In addition, we found for the first time that there were multiple differentially modified phosphorylation sites between PL and LL, including FASN (modified phosphorylation sites: 2120 To identify the phosphorylated proteins, the phosphorylated proteins and their differentially modified phosphorylation sites were annotated in detail by GO and KEGG pathways (Table S2).
A total of 2,454 proteins with quantitative information were classified by GO annotation based on three categories, namely, biological process, cellular component and molecular function. The results showed that 'metabolic process', 'biological regulation' and 'catalytic activity' were significantly enriched. 'Metabolic process'  (Fig. 3, Table S3).

Functional enrichment of the differentially phosphorylated proteins
Proteins that contained differentially modified phosphorylation sites were analyzed by GO enrichment analysis. From the analysis of the results, up-and downregulated phosphorylated site corresponding proteins were primarily enriched to translation, transportation and energy metabolism, such as 'translation', 'organic substance transport', 'Golgi vesicle transport' and, 'transporter activity' (Fig. 4).
KEGG enrichment analysis revealed that the proteins that contained differentially phosphorylated sites were remarkably enriched in the 'MAPK signaling pathway', 'Oxytocin signaling pathway' and 'Calcium signaling pathway' which are highly correlated with lactation (Fig. 5).

Validation of changes in the phosphorylated protein levels by western blot
Phosphorylation levels in PL and LL of ACACA (Ser80), EIF4EBP1 (Thr46), IRS1 (Ser312) and JUN (Ser63) were tested by western blot. The results showed that the phosphorylation levels of ACACA, EIF4EBP1 and IRS1 increased and JUN decreased in PL compared with LL (Fig. 6).

Discussion
Protein phosphorylation regulates many cellular processes in eukaryotic cells by posttranslational modification [13]. Protein phosphorylation reveals the effects of temperature on wheat leaves and spikelets [14], and the importance of protein phosphorylation has been demonstrated for different chloroplast functions [15]. Protein phosphorylation has an important role in sperm quality [16] and an important connection with lactation [5,17]. The level of phosphorylation constantly changes throughout the development of mammary tissues and during the different stages of lactation [3]. Therefore, in this study, differentially phosphorylated proteins in PL and LL were detected by LC-MS/MS and TMT. These results will contribute to the elucidation of the mechanism underlying the regulation of lactation in goat mammary tissues by phosphorylation.
The production of milk increases after parturition and subsequently reaches its peak [18]. These phenomena are dependent on a series of factors, including somatotropin, insulin and functional protein. Nevertheless, the effect of protein phosphorylation on goat mammary tissues during lactation remains unclear. Therefore, differentially phosphorylated proteins were explored in goat mammary tissues. A total of 1,938 differentially modified phosphorylation sites of 1,172 proteins were found between PL and LL. The western blot results of the phosphorylation levels of ACACA, EIF4EBP1, IRS1 and JUN showed that the phosphoproteome results were reliable, and these results are consistent with other studies [19,20]. The mTOR pathway is a key regulatory pathway in lactation. PRKAA, IRS1, RPS6KB, EIF4EBP1, TSC2, FASN and ACACA are located in the mTOR pathway. Previous studies showed that activation of mTOR pathway can enhance the synthesis of lactose and triglycerides [21,22], and researchers has shown that APOC3, FASN and ACACA proteins are highly correlated with the decomposition and utilization of fat [23][24][25]. Moreover, the inactivation of PLIN4 could decrease fat accumulation [26]. Insulin regulates glucose uptake and metabolism by binding to an insulin receptor [27], and insulin receptors (IRS1 and IRS2) are necessary for insulinfunction. In a previous study, Cryer, A et al. [28] demonstrated that activation of ACACA could promote fatty acid synthesis in adipocytes. Barber, M.C et al. [29] showed that lipid synthetic capacity during early lactation is associated with activity of FASN.
However, the differentially modified phosphorylation sites of the proteins identified in this study (FASN  TSC2 (1719, 1327 and 918)) have never been explored. Therefore, we speculate that the changes in these phosphorylated sites may be strongly related to lactation. The energy metabolism of mammary tissues increases remarkably during lactation [32]. GO enrichment showed several differentially phosphorylated proteins enriched in biological processes, which were related to energy metabolism, material transport and protein translation. These results provide an overview of dynamic changes during lactation. A series of signaling pathways were enriched in PL compared with LL, such as the 'calcium signaling pathway' (chx04020), 'oxytocin signaling pathway' (chx04921) and 'MAPK signaling pathway' (chx04010). Among them, the three signaling pathways of 'calcium signaling pathway', 'oxytocin signaling  Figure S2), IRS1(full-length blots/gels are presented in Figure S3), p-ACACA(fulllength blots/gels are presented in Figure S4), ACACA(full-length blots/gels are presented in Figure S5), p-EIF4EBP1(full-length blots/gels are presented in Figure S6), EIF4EBP1(full-length blots/gels are presented in Figure S7), p-JUN(full-length blots/gels are presented in Figure S8), JUN(full-length blots/ gels are presented in Figure S9) in PL and LL were detected by western blot. Full-length blots/gels of β-actin are presented in Figure S10 pathway', and 'MAPK signaling pathway' are highly related to lactation [33][34][35][36][37]. Therefore, we speculate that the changes in these phosphorylated sites in calcium signaling pathway, oxytocin signaling pathway, and MAPK signaling pathway may be related to lactation. These phosphorylated sites will become the research focus of our future work. This study provided a better understanding of the function of phosphorylated proteins in goat mammary gland epithelial tissues and provided a reference for future research.

Conclusions
The present study was the first to conduct sequence phosphorylated proteomics in the mammary tissues of dairy goats, and the findings of this study markedly increased the current understanding of posttranslational phosphorylation of protein events that occur during the progression of lactation. The western blot results of the phosphorylation levels of ACACA, EIF4EBP1, IRS1 and JUN showed that the phosphoproteome results were reliable. Multiple differentially modified phosphorylation sites of FASN, ACACA, mTOR, PRKAA, IRS1, RPS6KB, EIF4EBP1, TSC2, and JUN and proteins of the calcium signalling pathway, oxytocin signaling pathway, and MAPK signaling pathway might have potential research value in regulating goat lactation. We will study the function of these proteins and their differentially modified phosphorylation sites during lactation in the future. Therefore, the present study provides valuable reference information for further study of the molecular mechanism of lactation in goats.