The study was approved by the Institutional Review Board at Oregon State University, and informed consent was obtained from all mothers. Human foremilk and hindmilk samples were collected from four mothers (2–4 months of lactation). All mothers gave birth to term infants. Mothers with any signs of clinical mastitis and current antibiotic usage were excluded from the study. Before pumping, mothers thoroughly washed their hands and cleansed their nipples with a moist single-use paper towel. Breasts were pumped using a hospital-grade Medela Symphony electrical pump (Salem, OR, USA). All parts that contacted the mother or her milk were single use. Mothers were instructed to pump the first 10 mL of milk (foremilk) into an 80-mL breastmilk collection container (6109S-100; Medela) premarked at the 10-mL line. The container was prepared with protease inhibitor solution and held on ice [1 tablet of a protease inhibitor cocktail (cOmplete™ mini protease inhibitor cocktail tablets, 11836153001; Roche, Basel, Switzerland) added to 1 mL of MilliQ water, according to the manufacturer’s instructions]. The exact components of the antiprotease mixture are not disclosed; however, it inhibits a wide array of serine, cysteine, and metalloproteases as well as calpains. When the level of milk reached 10 mL, the mother immediately stopped pumping and passed the milk bottle to the researchers. The milk was immediately mixed well by gentle rocking, and divided into five 2-mL aliquots in 2-mL Eppendorf tubes and placed on dry ice within 1 min after expression. The mother kept pumping the milk into a second container until the flow rate dropped and the mother felt her breast was near empty [in all cases, later than 3 min after flow initiation, the standard definition for hindmilk (19)]. The mother then pumped 10 mL of hindmilk into a third container, which was handled as described for foremilk. After removing aliquots, the samples were immediately transported on dry ice to storage in a −80°C freezer within 30 min of expression. Milks were stored at −80°C for a maximum of three weeks prior to peptidomic analyses.

Aliquots of each milk sample (four foremilk and four hindmilk samples) were thawed on ice for approximately 30 min before removing the milk fat by centrifugation (3,000 × g for 10 min at 4°C). The infranates (400 µL) were collected by pipette. Milk proteins were precipitated from the skimmed milk by adding 400 µL of 24% trichloroacetic acid (240 g/L) followed by agitation of a vortex mixer for 10 s. Samples were centrifuged at 4,000 × g for 10 min at 4°C, and 600 µL of the supernatant containing the peptides were collected.

To remove trichloroacetic acid, salts, oligosaccharides and lactose from the peptide solution, the peptides were extracted via C18 reverse-phase preparative chromatography in 96-well plates (Glygen, Columbia, MD, USA). In general, unless otherwise noted, the volume of solution used in the different steps was 200 µL and the centrifuge parameters were 2,000 × g for 1 min at 4°C. Each step was repeated three times. To activate the C18 reverse-phase columns, 99% acetonitrile (ACN), 0.1% trifluoroacetic acid (TFA) solution was added. Columns were equilibrated with 1% ACN, 0.1% TFA before each sample was added to a separate well and centrifuged. This step was repeated two more times to load the total sample volume (600 µL) onto the column. The columns were washed with 1% ACN, 0.1% TFA to remove trichloroacetic acid, salts, lactose and oligosaccharides. The peptides were eluted with 80% ACN, 0.1% TFA. The peptide solutions were lyophilized using a freeze dry system (Labconco FreeZone 4.5 L, Kansas City, MO, USA). After freeze drying the samples, they were rehydrated in 40 µL of 0.1% FA in water.

The bicinchoninic acid assay for protein (ThermoFisher Scientific, Waltham, MA, USA)—conducted as described by the manufacturer with bovine serum albumin as reference protein—was used for determination of the protein/peptide concentration (24).

Liquid chromatography separation was performed on a Waters Nano Acquity UHPLC (Waters Corporation, Milford, MA, USA) with a nanospray source. Peptides were loaded (1 µL) onto a 180 µm × 20 mm, 5-µm bead 2 G nanoAcquity UPLC trap column (reverse-phase) for enrichment and online desalting and then onto a 100 µm × 100 mm, 1.7-µm bead Acquity UPLC Peptide BEH C18 column (Waters) connected to a Orbitrap Fusion Lumos (Thermo Scientific). To reduce the effects of analysis variation, the samples were loaded alternating between the two milk types. Peptides were eluted using a gradient of 0.1% FA in water (A) and 99.9% ACN, 0.1% FA (B) with a flow rate of 500 nL/min. The 120-min gradient consisted of 3−10% solvent B over 3 min, 10−30% solvent B over 99 min, 30−90% solvent B over 3 min, 90% solvent B for 4 min, 90−3% solvent B over 1 min then finally held at 3% solvent B for 10 min. Each sample analysis was followed by a 30-min column wash.

Spectra were collected in positive ionization mode with an electrospray voltage of 2,400 V. The mass spectrometer was set to scan masses between 400 and 1,500 m/z at a resolution of 120 K. The automatic gain control target was set to 4.0 × 10⁵, with a maximum injection time of 50 ms. The fragmentation mode was set to collision-induced dissociation and the collision energy was 35%. The mass spectrometer cycle time was 3 s, with data-dependent analysis and automated precursor peak selection. Precursors were excluded (10 ppm mass error) after one fragmentation and held for 1 min. Precursor ions were selected for fragmentation based on the following criteria: most intense peaks, ion intensity threshold 5.0 × 10³ and charge state 2–7. Fragments were detected with ion trap automatic scan range.

Spectra were analyzed by database searching in Thermo Proteome Discoverer (v2.1.0.81). The data were searched against an in-house human milk protein sequence database. Phosphorylation of serine and threonine and oxidation of methionine was allowed as potential modifications. Peptides identified with high confidence (P < 0.01) were accepted. Each peptide sequence with multiple modifications was grouped into a single entry for counts. Counts measured the number of unique peptide sequences identified in a sample. Abundance measured the area under the curve of the eluted peak (ion intensity). As the abundance measurement is based on relative ion intensity, it is only an approximation of the amount of peptides in the samples.

Peptides were mapped to their parent protein sequence of β-casein, osteopontin, and α_s1-casein using an in-house tool (PepEx) (15), which can be accessed at http://mbpdb.nws.oregonstate.edu/pepex/. PepEx gives a visualization of the position of peptides identified by totaling the abundance or count of each amino acid from the peptidomic data.

The identified peptides in human milk samples were investigated for literature-identified bioactive peptides using our recently created Milk Bioactive Peptide Database (MBPDB) (5), which is a comprehensive database of all milk protein-derived bioactive peptides reported in literature and can be accessed at http://mbpdb.nws.oregonstate.edu/. The search performed was conducted as a sequence search, which searches for bioactive peptides matching the input peptide sequence. To predict the bioactive potential of each identified peptide, we used two online prediction tools—CAMPR3 (antimicrobial prediction) (25) and Peptide Ranker (general bioactivity prediction) (26). The CAMPR3 prediction was conducted by inputting the list of identified peptide sequences and searching for antimicrobial potential with four algorithms (SVM, random forest, artificial neural network, and discriminant analysis). For general bioactivity (Peptide Ranker), peptide sequences were input with no further options available and results were scores from 0 to 1 for each peptide with 1 being most likely to be bioactive.

To investigate which proteases were involved in cleaving human milk proteins into peptides, three approaches were used. In Proteome Discoverer, peptides were identified with annotated sequences. With this information, Excel was used to determine the count and abundance of which amino acids were present in the P1 and P1′ positions of the N- and C-termini of each peptide identified. The cleavage specificity of known milk proteases were mapped to these amino acids for visual interpretation in the figures [based on the most abundant amino acids listed at P1 and P1′ on Merops (27)] and included elastase (MEROPS ID: S01.131), kallikrein 6 (KLK 6) (MEROPS ID: S01.236), kallikrein 11 (KLK 11) (MEROPS ID: S01.257), plasmin (MEROPS ID: S01.233), and thrombin (MEROPS ID: S01.217). To determine which peptides were possibly released by exopeptidase activity, the number of peptides that differed from a related sequence by one amino acid at either the N- or C-termini was counted.

To predict which enzymes were active in milk using more complex cleavage specificity by accounting for amino acids from P4 to P2′ and P4 to P4′, the two online predictors EnzymePredictor and Proteasix, respectively, were used.

EnzymePredictor¹ required input of the peptide amino acid sequences and their associated protein ID as found at uniprot.org (28). This information is used to identify the amino acids located at the P4, P3, P2, P1, P1′, and P2′ positions of both the N- and C-terminal cleavage sites of each peptide. Based on this information and the cleavage specificity patterns of the selected proteases, EnzymePredictor predicts which enzyme could yield each peptide. For each protease identified, EnzymePredictor provides the count of N-terminal-, C-terminal- and total cleavages assigned to each enzyme, the number of times an enzyme could have cleaved within all the protein sequences in the library of proteins observed and the odds ratio (comparing the number of predicted cleavages based on sequence alone with the actual peptide cleavages assigned to each protease).

Proteasix² predicted which enzymes were involved in the peptide cleavage sites using amino acid cleavage patterns from Merops (27). Proteasix required the input protein ID and start and stop amino acid number of the peptide within the protein. The enzymes searched against were plasmin (PLG), neutrophil elastase (ELANE), thrombin (F2), KLK 6, and KLK 11 based on the amino acids positioned from P4 to P4′ at the cleavage site. Predicted and observed cleavages were combined in the result output.

Analyses were carried out using the statistical program RStudio version 1.0.136. A linear mixed model with Tukey’s HSD post hoc test was used to adjust for multiple comparisons between foremilk and hindmilk. Significant difference was defined as P < 0.05. Results are presented as least square mean ± SE. For individual comparison of the peptide abundance in foremilk compared to hindmilk, a two-tailed pairwise t-test was conducted using the t-test function in Excel.

Protein concentration in breast milk was determined by the bicinchoninic acid assay and was not significantly different between the foremilk (10.2 ± 0.4 mg/mL) and hindmilk (9.5 ± 0.4 mg/mL). The coefficient of variation for total protein concentration was 8.6%. On average, 474 ± 14 peptides were identified in the foremilk and 591 ± 14 identified in the hindmilk (Figure 1A). Significantly more peptides were identified in the hindmilk than the foremilk (P = 0.0106). A total of 370.3 ± 11.7 peptides were found in both the foremilk and hindmilk of each mother and comprised 98.3 ± 1.1% of the foremilk peptide abundance and 96.7 ± 1.1% of the hindmilk peptide abundance. Of these peptides, 135 peptides deriving from 10 different proteins were identified in both fore and hindmilk of all four mothers (Tables 1–3). These accounted for 91.3 ± 3.5 and 86.8 ± 3.5% of the total abundance of peptides in foremilk and hindmilk, respectively. The abundance of all peptides in each sample were not significantly different between foremilk and hindmilk (5.9 × 10¹¹ ± 7.1 × 10¹⁰ and 7.1 × 10¹¹ ± 7.1 × 10¹⁰, respectively; Figure 1B). However, of all the peptides identified, 36 were found to be significantly more abundant in hindmilk than foremilk (P < 0.05; Table 4). No peptides were significantly more abundant in foremilk than hindmilk.

The peptides derived on average from 24.5 ± 0.8 proteins in the foremilk and 27.0 ± 0.8 proteins in the hindmilk (P = 0.127). In total, peptides from 42 proteins were identified. The peptides derived mostly from β-casein (51.1 ± 7.4% of total abundance), osteopontin (31.1 ± 11.3%), α_s1-casein (11.1 ± 6.1%), polymeric immunoglobulin receptor (3.2 ± 1.4%), butyrophilin (2.6 ± 2.6%), mucin-1 (0.3 ± 0.3%), bile salt-activated lipase (0.1 ± 0.08%), fibrinogen α-chain (0.08 ± 0.07%), and perilipin-2 (0.1 ± 0.1%). The remaining proteins accounted for 0.3 ± 0.1% of total abundance. The highest number of peptides derived from β-casein, and more peptides were identified from β-casein in hindmilk (199 ± 9 peptides) than in foremilk (149 ± 9 peptides) (P = 0.0312), although there was no difference in abundance (Figure 2). The second highest number of milk peptides derived from osteopontin. The count of peptides from osteopontin was higher in hindmilk (121 ± 2) than foremilk (111 ± 2) (P = 0.0334), but there was no difference in abundance. α_s1-Casein-derived peptides were higher in hindmilk than foremilk by count (65 ± 1 vs. 53 ± 1, P = 0.00184) and tended to be higher in hindmilk than foremilk by abundance (P = 0.0587). For mucin-1, both the count (P = 0.0138) and abundance (P = 0.0129) of peptides were higher in hindmilk than foremilk (19 ± 1 vs. 12 ± 1 for count). Peptides deriving from polymeric immunoglobulin receptor and butyrophilin did not differ in count or abundance between foremilk and hindmilk.

The total number of peptides deriving from β-casein, osteopontin, and α_s1-casein accounted for most of the peptide abundance—94.7 ± 1.2 and 91.9 ± 1.2% in foremilk and hindmilk, respectively. The peptides identified in foremilk and hindmilk were highly similar, shown for these three proteins from which most peptides were derived (Figures 3–5).

To identify which of the known proteases in milk were involved in digestion of the foremilk and hindmilk, we used three bioinformatic approaches. Each peptide typically has two cleavage sites, unless the N-terminus or C-terminus of the peptide coincides with the N-terminus or C-terminus of the protein sequence. Foremilk and hindmilk had 879 ± 27 and 1,101 ± 27 cleavage sites, respectively.

At the P1 amino acid position, the amino acid most often present was lysine followed by arginine, alanine, and serine (Figures 6 and 7). All of these amino acids were present at more cleavage sites in hindmilk than foremilk by count. Compared to peptides with lysine, arginine, and alanine at P1, peptides with serine at P1 were much less abundant. Cleavage sites with a lysine and arginine at P1 correspond with the known site specificity for plasmin, KLK 6, KLK 11, and thrombin. Peptides with proline, valine, glutamine, asparagine, isoleucine, phenylalanine, and tryptophan at the P1 position were also significantly higher in hindmilk compared with foremilk, by count. Of these, tyrosine, asparagine, and isoleucine at the P1 position of cleavage sites correspond with the site specificity of elastase. By abundance, only serine (P = 0.0258) and leucine (P = 0.0004) at the P1 position were significantly greater in hindmilk compared to foremilk.

At the P1′ position of each cleavage site, the amino acids alanine (P = 0.0258), leucine (P = 0.0155), glutamic acid (P = 0.0473), proline (P = 0.0057), glutamine (P = 0.0115), tyrosine (P = 0–0249), aspartic acid (P = 0.0071), and methionine (P = 0.0253) were significantly higher in hindmilk than foremilk by count. Of these, only alanine at the P1′ position corresponds with one of our listed protease cleavage preferences (thrombin).

Though examining the P1 and P1′ amino acids can indicate matches with enzyme cleavage specificity, enzyme cleavage patterns are often more complex than their P1 and P1′ preferences. Two tools predict which proteases are most likely responsible for peptide cleavages based on a more complex cleavage specificity pattern. Based on data from EnzymePredictor analysis, plasmin was involved in the most cleavages by count in both foremilk and hindmilk (Table 5). There were no significant differences in the odds ratios (i.e., the likelihood of protease activity) for any of the proteases between foremilk and hindmilk. However, the C-terminal (P = 0.0007) and total cleavage counts (P = 0.0162) assigned to plasmin and the N-terminal (P = 0.0028) and total cleavage counts (P = 0.0085) assigned to elastase were higher in hindmilk compared with foremilk. Thrombin was only found to be responsible for one cleavage; however, the number of expected cleavages was low as well, resulting in a high odds ratio. KLK 6 and 11 are not available for search within Enzyme Predictor.

Proteasix also predicted that plasmin cleaved the largest number of peptides in the dataset. By this analysis, there were no differences for any of the predicted protease activities between foremilk and hindmilk; however, cleavages assigned to elastase tended to be higher in hindmilk compared with foremilk (P = 0.0572) (Table 6). More cleavage sites were unassigned in hindmilk than foremilk (P = 0.0137), as only 15.5 ± 0.4% of potential cleavages were assigned by Proteasix in hindmilk compared with 17.6 ± 0.4% in foremilk. Several cleavage sites were assigned to two or more of plasmin, KLK 6 and thrombin as Proteasix was not able to distinguish among these enzymes for a number of the cleavage sites. On average, 14 cleavage sites overlapped between plasmin and thrombin, 3 overlapped between KLK 6 and plasmin, 1 overlapped between thrombin and KLK 6 and 1 cleavage site overlapped among all three enzymes.

A large portion of the identified peptides could not be assigned to a specific enzyme. However, several of these peptides matched other peptides with sequential removal of single amino acids. These could be potentially assigned to exopeptidase activity. N-terminal sequential removal of amino acids, potentially representing aminopeptidase activity, was significantly higher in hindmilk (357.0.5 ± 9.4) than foremilk (294.5 ± 9.4) (P = 0.0190) by count, but it was not higher by abundance (4.6 × 10¹¹ ± 4.9 × 10¹⁰ vs. 3.7 × 10¹¹ ± 4.9 × 10¹⁰). C-terminal sequential removal of amino acids, potentially representing carboxypeptidase activity, was similarly higher by count in hindmilk (384.0 ± 10.2 peptides) than foremilk (302.8 ± 10.2) (P = 0.0109), but not abundance (3.5 × 10¹¹ ± 3.5 × 10¹⁰ vs. 2.7 × 10¹¹ ± 3.5 × 10¹⁰).

Searching the MBPDB for identified peptide homology with known bioactive peptides revealed two peptides with 100% alignment to a known bioactive peptide (Table 7). These two peptides derived from β-casein f105-117 (SPTIPFFDPQIPK) and f185-211 (QELLLNPTHQIYPVTQPLAPVHNPISV), which were identified in 7 and 5 of the eight milk samples, respectively. β-Casein f105-117 can increase cell proliferation (29), and β-casein f185-211 is antimicrobial (30). The abundances of these two peptides were not different between the foremilk and hindmilk. With an 80% homology threshold, an additional five peptides matched to β-casein f105-117 and an additional 21 peptides matched to β-casein f185-211. The total abundance of the peptides with >80% homology to the antimicrobial β-casein f185-211 was higher in hindmilk than foremilk (P = 0.0429), but not the total count (P = 0.1602). Neither the total abundance nor the total count of the peptides with >80% homology to β-casein f105-117 differed between foremilk and hindmilk. This search also identified six peptides with >80% homology to known ACE-inhibitory peptides and one additional peptide with >80% homology to a known antimicrobial peptide.

Bioactivity prediction software (CAMPR3) predicted a higher number of antimicrobial peptides in hindmilk (81 ± 4) compared to foremilk (62 ± 4) (P = 0.0381). The abundance of these peptides tended to be higher in hindmilk than foremilk (P = 0.0927) (Table 7). However, the percentage of predicted antimicrobial peptides among the total identified peptides did not differ between foremilk (13.2 ± 0.6% by count, 1.9 ± 1.4% by abundance) and hindmilk (13.7 ± 0.6% by count, 4.3 ± 1.4% by abundance).

The count of general bioactive peptides as predicted by Peptide Ranker was also higher in hindmilk (65 ± 2) compared to foremilk (56 ± 2) (P = 0.0306). The abundance of these peptides tended to be higher in hindmilk than foremilk (P = 0.0773) (Table 7). However, the percentage of predicted bioactive peptides compared to the total did not differ between foremilk and hindmilk by count (11.7 ± 0.3 vs. 11.0 ± 0.3%, respectively). The percentage abundance of predicted general bioactive peptides compared to the total tended to be higher in hindmilk (13.4 ± 0.4%) than foremilk (11.8 ± 0.4%) (P = 0.0555).