Advertisement

Nonlinear modeling to describe the pattern of 15 milk protein and nonprotein compounds over lactation in dairy cows

Open ArchivePublished:August 04, 2021DOI:https://doi.org/10.3168/jds.2020-20086

      ABSTRACT

      The protein profile of milk includes several caseins, whey proteins, and nonprotein nitrogen compounds, which influence milk's value for human nutrition and its cheesemaking properties for the dairy industry. To fill in the gap in current knowledge of the patterns of these individual nitrogenous compounds throughout lactation, we tested the ability of a parametric nonlinear lactation model to describe the pattern of each N compound expressed qualitatively (as % of total milk N), quantitatively (in g/L milk), and as daily yield (in g/d). The lactation model was tested on a data set of detailed milk nitrogenous compound profiles (15 fractions—12 protein traits and 3 nonproteins—for each expression mode: 45 traits) obtained from 1,342 cows reared in 41 multibreed herds. Our model was a modified version of Wilmink's model, often used for describing milk yield during lactation because of its reliability and ease of parameter interpretation from a biological point of view. We allowed the sign of the persistency coefficient (parameter c) that explained the variation in the long-term milk component (parameter a) to be positive or negative. We also allowed the short-term milk component (parameter b) to be positive or negative, and we estimated a specific speed of adaptation parameter (parameter k) for each trait rather than assumed a value a priori, as in the original model (k = 0.05). These 4 parameters were included in a nonlinear mixed model with cow breed and parity order as fixed effects, and herd-date as random. Combinations of the positive and negative signs of the b and c parameters allowed us to identify 4 differently shaped lactation curves, all found among the patterns exhibited by the nitrogenous fractions as follows: the “zenith” curve (with a maximum peak; for milk yield and 10 other N traits), the “nadir” curve (with a minimum point; for 20 traits, including almost all those expressed in g/L of milk), the “downward” curve (continuously decreasing; for 14 traits, including almost all those in g/d), and the “upward” curve (continuously increasing; only for κ-casein, in % N). Direct estimation of the k parameters specific to each trait showed the large variability in the adaptation speed of fresh cows and greatly increased the model's flexibility. The results indicated that nonlinear parametric mathematical models can effectively describe the different and complex patterns exhibited by individual nitrogenous fractions during lactation; therefore, they could be useful tools for interpreting milk composition variations during lactation.

      Key words

      INTRODUCTION

      Over the last almost 100 yr, many empirical models have been developed to describe the lactation curve of dairy cows. The curves described by these models generally have the following 2 phases: (1) a first ascending phase after calving that describes the rapid adaptation of the udder to milk production up to the maximum point (zenith value), the peak of lactation; and (2) a slowly descending phase after the peak up to the dry period that describes the persistency of lactation (
      • Gengler N.
      Persistency of lactation yields: A review.
      ;
      • Macciotta N.P.P.
      • Vicario D.
      • Cappio-Borlino A.
      Detection of different shapes of lactation curve for milk yield in dairy cattle by empirical mathematical models.
      ). Lactation models are of wide interest not only for describing and understanding the physiology of the mammary gland, but also for making management decisions at the farm level and genetic improvements at the population level (
      • Bouallegue M.
      • M'Hamdi N.
      Mathematical modeling of lactation curves: A review of parametric model.
      ).
      When the aim is to investigate the mean pattern of a large population of animals, parametric models are preferred as they can describe the lactation curve with a small number of parameters (
      • Macciotta N.P.P.
      • Dimauro C.
      • Rassu S.P.G.
      • Steri R.
      • Pulina G.
      The mathematical description of lactation curves in dairy cattle.
      ). Different parametric models are characterized by differences in their levels of complexity (number of parameters), in the ease with which physiological and economic functions can be identified with parameters, and in their statistical performances, such as precision and accuracy (
      • Silvestre A.M.
      • Petim-Batista F.
      • Colaço J.
      The accuracy of seven mathematical functions in modeling dairy cattle lactation curves based on test-day records from varying sample schemes.
      ). An example of a parametric model that offers a good compromise between these different evaluation criteria is that proposed by
      • Wilmink J.B.M.
      Adjustment of test-day milk, fat and protein yield for age, season and stage of lactation.
      . This model has sometimes been criticized for its poor flexibility, in particular in describing the lactation curves of individual cows. However, it has also been acclaimed for its accuracy and small standard error in large population studies, and for the ease of attributing meaning to its parameters (
      • Macciotta N.P.P.
      • Vicario D.
      • Cappio-Borlino A.
      Detection of different shapes of lactation curve for milk yield in dairy cattle by empirical mathematical models.
      ;
      • Strucken E.M.
      • de Koning D.J.
      • Rahmatalla S.A.
      • Brockmann G.A.
      Lactation curve models for estimating gene effects over a timeline.
      ).
      Similar to other parametric models, the Wilmink model has been applied not only to daily raw milk production, but also to milk corrected for fat or protein content, as well as to the daily production of milk fat and protein, which exhibit similar patterns throughout lactation. When fat and protein (and also other milk nutrients) are expressed in terms of milk composition (g/L or %), the pattern throughout lactation often appears to be the reverse of the curve for daily production, with an initial rapid decrease in early lactation to a minimum point (nadir value) and a subsequent gradual increase up to the dry period. In the case of fat content, the nadir point was observed about 3 wk after the peak of daily production, whereas it almost coincided with the lactation peak for protein content (
      • Wood P.D.P.
      Algebraic models of the lactation curves for milk, fat and protein production, with estimates of seasonal variation.
      ;
      • Stanton T.L.
      • Jones L.R.
      • Everett R.W.
      • Kachman S.D.
      Estimating milk, fat, and protein lactation curves with a test day model.
      ). The patterns throughout lactation of other milk nutrient contents (e.g., lactose, protein fractions, urea, fatty acids, minerals) are less well known and have only rarely been described using parametric models (
      • Arnould V.M.R.
      • Hammami H.
      • Soyeurt H.
      • Gengler N.
      Short communication: Genetic variation of saturated fatty acids in Holsteins in the Walloon region of Belgium.
      ;
      • Gołebiewski M.
      • Brzozowski P.
      • Gołebiewski Ł.
      Analysis of lactation curves, milk constituents, somatic cell count and urea in milk of cows by the mathematical model of wood.
      ;
      • Baul S.
      • Cziszter L.
      • Acatincai S.
      • Gavojdian D.
      • Erina S.
      • Marcu A.
      • Buzamat G.
      • George G.
      Seasonal influences on milk yield and composition dynamics during a normal lactation in dairy cows: Milk yield, fat and protein percentage.
      ).
      Studies dealing with the patterns of individual protein fractions during lactation are scarce in the literature, even though these nutrients are of great interest for their value for human nutrition (
      • Crittenden R.G.
      • Bennett L.E.
      Cow's milk allergy: A complex disorder.
      ;
      • Weimann C.
      • Meisel H.
      • Erhardt G.
      Short communication: Bovine κ-casein variants result in different angiotensin I converting enzyme (ACE) inhibitory peptides.
      ;
      • Sebastiani C.
      • Arcangeli C.
      • Ciullo M.
      • Torricelli M.
      • Cinti G.
      • Fisichella S.
      • Biagetti M.
      Frequencies evaluation of β-casein gene polymorphisms in dairy cows reared in central Italy.
      ) and their large influence on the cheesemaking properties of milk (
      • Wedholm A.
      • Larsen L.B.
      • Lindmark-Månsson H.
      • Karlsson A.H.
      • Andrén A.
      Effect of protein composition on the cheese-making properties of milk from individual dairy cows.
      ;
      • Cipolat-Gotet C.
      • Cecchinato A.
      • Malacarne M.
      • Bittante G.
      • Summer A.
      Variations in milk protein fractions affect the efficiency of the cheese-making process.
      ;
      • Amalfitano N.
      • Cipolat-Gotet C.
      • Cecchinato A.
      • Malacarne M.
      • Summer A.
      • Bittante G.
      Milk protein fractions strongly affect the patterns of coagulation, curd firming, and syneresis.
      ). Some studies include DIM classes in the model (
      • Ng-Kwai-Hang K.F.
      • Hayes J.F.
      • Moxley J.E.
      • Monardes H.G.
      Environmental influences on protein content and composition of bovine milk.
      ;
      • Ostersen S.
      • Foldager J.
      • Hermansen J.E.
      Effects of stage of lactation, milk protein genotype and body condition at calving on protein composition and renneting properties of bovine milk.
      ;
      • Maurmayr A.
      • Pegolo S.
      • Malchiodi F.
      • Bittante G.
      • Cecchinato A.
      Milk protein composition in purebred Holsteins and in first/second-generation crossbred cows from Swedish Red, Montbeliarde and Brown Swiss bulls.
      ), and a few studies use parametric functions.
      • Heck J.M.L.
      • Schennink A.
      • Van Valenberg H.J.F.
      • Bovenhuis H.
      • Visker M.H.P.W.
      • Van Arendonk J.A.M.
      • Van Hooijdonk A.C.M.
      Effects of milk protein variants on the protein composition of bovine milk.
      and
      • Rutten M.J.M.
      • Bovenhuis H.
      • Heck J.M.L.
      • van Arendonk J.A.M.
      Predicting bovine milk protein composition based on Fourier transform infrared spectra.
      included the Wilmink function in a mixed model to study the content of the major milk protein fractions corrected for the effect of stage of lactation, but they analyzed the protein fractions qualitatively as percent of total milk nitrogen, and not as their content in milk. Moreover, they did not show and discuss the trend of these fractions across the entire lactation period. Information on NPN constituents in milk are also scarce. Among the NPN compounds, only the content of milk urea has been more widely investigated across lactation (
      • Wood G.M.
      • Boettcher P.J.
      • Jamrozik J.
      • Jansen G.B.
      • Kelton D.F.
      Estimation of genetic parameters for concentrations of milk urea nitrogen.
      ;
      • Stoop W.M.
      • Bovenhuis H.
      • van Arendonk J.A.M.
      Genetic parameters for milk urea nitrogen in relation to milk production traits.
      ;
      • Rzewuska K.
      • Strabel T.
      Genetic parameters for milk urea concentration and milk traits in Polish Holstein-Friesian cows.
      ).
      The present study aimed to analyze the lactation pattern of the detailed profiles of milk nitrogen fractions (15 proteins and NPN components) expressed as proportions of total nitrogen (% N), total contents in milk (g/L), and daily production (g/d). Our aim was to test a nonlinear lactation model for describing detailed milk nitrogen compounds characterized by very different patterns during lactation in practical conditions through a large survey of cows reared in several multibreed herds under different dairy systems.

      MATERIALS AND METHODS

      Experimental Design, Animals, and Milk Sampling

      This work is part of the Cowplus project, a large survey designed to investigate the effects of the major sources of variation in milk composition and cheesemaking properties. Only multibreed herds were selected for the survey to avoid confounding the effects of breed and herd. These herds (41 in total; 2–4 breeds per herd) were located in the provinces of Trento and Bolzano (Bozen–Süd Tirol) in northeastern Italy and represent the different dairy systems and production intensity levels in the region. The dairy systems ranged from the very traditional type, with tied animals and the use of hay and concentrates, to the most modern ones, with loose housing and TMR with or without silages. Pasture in the area was practiced only by very traditional farms through transhumance to temporary summer farms in the highland pasture. However, milk recording was not performed in these areas, and thus the sampled cows were all kept in the lowland on indoor feeding. All of the cows of the smaller herds and the less numerous breeds were sampled (minimum of 7 cows in a herd), whereas the sampling percentage decreased with increasing herd and breed population sizes up to a maximum of 80 cows/herd per sampling date, which is the maximum number of fresh milk samples that our laboratory could analyze and evaluate for cheesemaking aptitude. For the present study, we sampled 1,342 cows once as follows: 422 purebred Holstein-Friesians, 598 Brown Swiss, 148 Simmentals, and 174 cows of local breeds (Rendena and Grey Alpine). The cows sampled were of different parities (1–12) and stages of lactation (5–365 DIM); therefore, the colostrum was excluded. The cows were divided in 3 classes of parity order, each ranging from 5 to 365 DIM (first parity: 442 cows with 168 ± 98 DIM; second parity: 327 cows with 159 ± 99 DIM; ≥ third parity: 573 cows with 156 ± 96 DIM; ±SD). Cows presenting clinical symptoms of mastitis or other health problems on the day of sampling were excluded.
      A 2-L milk sample was taken from each cow during the evening milking. The sample was collected from the cow's complete milking after stirring to ensure its representativeness. The daily milk yield of each cow was obtained with the sum of morning and evening milking records. Immediately after sampling, the milk samples were transported to the Milk Laboratory of the Department of Agronomy, Food, Natural Resources, Animals and the Environment (DAFNAE) of the University of Padova (Legnaro, Padua, Italy), where they were divided into different aliquots (after homogenization through stirring) for the following laboratory analyses: (1) fresh aliquots were designated for coagulation and curd firming analysis (50 mL) and for individual model cheesemaking (1.5 L), as illustrated and discussed by
      • Stocco G.
      • Cipolat-Gotet C.
      • Bobbo T.
      • Cecchinato A.
      • Bittante G.
      Breed of cow and herd productivity affect milk composition and modeling of coagulation, curd firming, and syneresis.
      ,
      • Stocco G.
      • Cipolat-Gotet C.
      • Gasparotto V.
      • Cecchinato A.
      • Bittante G.
      Breed of cow and herd productivity affect milk nutrient recovery in curd, and cheese yield, efficiency and daily production.
      ); (2) 50-mL fresh aliquots were analyzed for milk composition (MilkoScan FT6000, Foss), SCC (Fossomatic Minor, Foss), and bacterial content (
      • Bobbo T.
      • Ruegg P.L.
      • Stocco G.
      • Fiore E.
      • Gianesella M.
      • Morgante M.
      • Pasotto D.
      • Bittante G.
      • Cecchinato A.
      Associations between pathogen-specific cases of subclinical mastitis and milk yield, quality, protein composition, and cheese-making traits in dairy cows.
      ); and (3) aliquots of 1.5 mL were frozen at −80°C until reversed-phase (RP)-HPLC analysis.

      Analysis of Detailed Milk Protein and Nonprotein Profiles

      Detailed profiling of milk protein and nonprotein nitrogenous compounds was based on 3 groups of information that were reported in our previous study on the effects of breed (
      • Amalfitano N.
      • Stocco G.
      • Maurmayr A.
      • Pegolo S.
      • Cecchinato A.
      • Bittante G.
      Quantitative and qualitative detailed milk protein profiles of 6 cattle breeds: Sources of variation and contribution of protein genetic variants.
      ). Briefly, the first group was 3 N compounds (protein, CN, and urea) analyzed with the Fourier-transform infrared spectroscopy (FTIR) methods officially used for milk recording (
      ICAR (International Committee for Animal Recording)
      International agreement of recording practices—Guidelines approved by the General Assembly held in Cork, Ireland on June 2012.
      ). The second group was the analysis by RP-HPLC and the method described in detail by
      • Maurmayr A.
      • Cecchinato A.
      • Grigoletto L.
      • Bittante G.
      Detection and quantification of αS1-, αS2-, β-, κ-casein, α-lactalbumin, β-lactoglobulin and lactoferrin in bovine milk by reverse-phase high- performance liquid chromatography.
      , which provided results on the contents of the following 7 protein fractions: αS1-CN, αS2-CN, β-CN, glycosylated κ-CN (glyco-κ-CN), carbohydrate-free κ-CN (CF-κ-CN), β-LG, and α-LA. The third group of traits was 5 traits calculated from the other traits as follows: true protein was calculated as the sum of the 7 protein fractions obtained by HPLC analysis; total κ-CN was the sum of glyco- κ-CN and CF-κ-CN; whey proteins were the sum of β-LG and α-LA; milk NPN compounds were calculated as the difference between milk total protein and true protein; and minor NPN compounds as the difference between milk NPN compounds and urea.
      Briefly, the HPLC equipment consisted of an Agilent 1260 Infinity Series chromatograph (Agilent Technologies) equipped with an RP analytical column (Aeris WIDEPORE XB-C8; Phenomenex) packed with large-pore core-shell particles (3.6 μm, 200Å, 250 × 2.1 i.d.), and a Security Guard ULTRA Cartridge System (UHPLC WIDEPORE C8, 2.1 mm i.d.; Phenomenex).
      To bring the N compounds obtained by FTIR into line with the protein fractions measured by liquid chromatography, the HPLC protein fractions were multiplied by the ratio between the FTIR CN and the sum of all the HPLC CN fractions. The nitrogenous fractions were expressed in 3 different ways as follows: as a percentage of the total N content (% N), as total content per liter of milk (g/L), and in grams per day of lactation (g/d). This last value was obtained by multiplying the daily milk yield by the fraction content of milk, thus assuming that milk composition was not different in morning and evening milking. All the traits analyzed are listed in Table 1, which also provides summaries of the analytical method, descriptive statistics (means and SD), and the statistical analyses.
      Table 1Descriptive statistics and results of the nonlinear mixed model based on herd-date lactation curve modeling parameters, breed, and parity (significance) for milk yield (kg/d) and milk protein traits expressed as relative proportion (% N), milk content (g/L), and daily production (g/d)
      Trait
      Glyco-κ-CN = glycosylated κ-CN; CF-κ-CN = carbohydrate-free κ-CN; CN = Σ CN (αS1-CN + αS2-CN + β-CN+ κ-CN); whey proteins = Σ whey proteins (β-LG + α-LA); κ-CN = glyco-κ-CN + CF-κ-CN; NPN compounds = NPN compounds (Σ urea + minor NPN compounds); examples of minor NPN compounds include small peptides, ammonia, creatine, creatinine.
      Method
      FTIR = Fourier-transform infrared spectroscopy. For method details please see Materials and Methods.
      MeanSDBreed
      Four different breeds (Holstein, Brown Swiss, Simmental, and local breeds).
      Parity
      Three parity classes (first, second, ≥third).
      Herd-date
      Herd-date: proportion of total variance explained by the herd.
      RMSE
      RMSE = root mean square error.
      Protein composition, % N
       True proteinBy sum93.62.9
      P < 0.001.
      122.6
      CNFTIR78.61.1
      P < 0.01
      P < 0.001.
      260.8
      αS1-CNHPLC26.161.7
      P < 0.001.
      121.4
      αS2-CNHPLC7.81.3
      P < 0.001.
      P < 0.01
      71.2
      β-CNHPLC29.82.2
      P < 0.001.
      P < 0.001.
      171.7
      κ-CNBy sum15.02.4
      P < 0.001.
      P < 0.01
      62.0
      Glyco-κ-CNHPLC5.41.7
      P < 0.001.
      P < 0.001.
      61.5
      CF-κ-CNHPLC9.41.5
      P < 0.001.
      P < 0.001.
      21.2
      Whey proteinsBy sum15.12.7
      P < 0.001.
      P < 0.01
      102.4
      β-LGHPLC12.82.7
      P < 0.001.
      P < 0.01
      82.4
      α-LAHPLC2.10.4
      P < 0.001.
      320.3
       NPN compoundsBy difference6.82.7
      P < 0.001.
      122.4
      Milk ureaFTIR2.00.8
      P < 0.01
      820.3
      Minor NPN compoundsBy difference4.72.6
      P < 0.001.
      62.4
      Milk composition, g/L
       Total proteinFTIR35.64.1
      P < 0.001.
      P < 0.001.
      372.4
      True proteinBy sum33.33.9
      P < 0.001.
      P < 0.001.
      342.4
      CNFTIR27.93.1
      P < 0.001.
      P < 0.001.
      351.9
      αS1-CNHPLC9.31.2
      P < 0.001.
      P < 0.001.
      340.8
      αS2-CNHPLC2.80.6
      P < 0.001.
      P < 0.001.
      120.5
      β-CNHPLC10.61.2
      P < 0.01
      P < 0.001.
      200.9
      κ-CNBy sum5.41.2
      P < 0.001.
      110.8
      Glyco-κ-CNHPLC1.90.7
      P < 0.001.
      P < 0.01
      50.6
      CF-κ-CNHPLC3.40.7
      P < 0.001.
      P < 0.01
      110.5
      Whey proteinsBy sum5.41.2
      P < 0.001.
      P < 0.05
      190.9
      β-LGHPLC4.51.1
      P < 0.001.
      P < 0.05
      170.9
      α-LAHPLC0.80.1
      P < 0.001.
      P < 0.01
      280.1
      NPN compoundsBy difference2.41.0
      P < 0.001.
      130.9
      Milk ureaFTIR0.30.1
      P < 0.001.
      P < 0.01
      850.0
      Minor NPN compoundsBy difference1.71.0
      P < 0.001.
      50.9
      Milk yield, kg/dMilk recording24.98.7
      P < 0.001.
      P < 0.001.
      624.5
      Protein production, g/dBy multiplying
       Total proteinFTIR890316
      P < 0.001.
      P < 0.001.
      70157
      True proteinBy sum830295
      P < 0.001.
      P < 0.001.
      70146
      CNFTIR700248
      P < 0.001.
      P < 0.001.
      69125
      αS1-CNHPLC23387
      P < 0.001.
      P < 0.001.
      6943
      αS2-CNHPLC6826
      P < 0.001.
      P < 0.05
      5316
      β-CNHPLC26293
      P < 0.001.
      P < 0.001.
      6351
      κ-CNBy sum13250
      P < 0.001.
      P < 0.001.
      6229
      Glyco-κ-CNHPLC4720
      P < 0.001.
      P < 0.001.
      4014
      CF-κ-CNHPLC8434
      P < 0.001.
      P < 0.001.
      5820
      Whey proteinsBy sum13152
      P < 0.001.
      P < 0.001.
      6029
      β-LGHPLC11146
      P < 0.001.
      P < 0.001.
      5527
      α-LAHPLC197
      P < 0.001.
      P < 0.001.
      624
      NPN compoundsBy difference5729
      P < 0.001.
      P < 0.01
      2923
      Milk ureaFTIR63
      P < 0.001.
      P < 0.001.
      741
      Minor NPN compoundsBy difference3925
      P < 0.001.
      1722
      1 Glyco-κ-CN = glycosylated κ-CN; CF-κ-CN = carbohydrate-free κ-CN; CN = Σ CN (αS1-CN + αS2-CN + β-CN+ κ-CN); whey proteins = Σ whey proteins (β-LG + α-LA); κ-CN = glyco-κ-CN + CF-κ-CN; NPN compounds = NPN compounds (Σ urea + minor NPN compounds); examples of minor NPN compounds include small peptides, ammonia, creatine, creatinine.
      2 FTIR = Fourier-transform infrared spectroscopy. For method details please see Materials and Methods.
      3 Four different breeds (Holstein, Brown Swiss, Simmental, and local breeds).
      4 Three parity classes (first, second, ≥third).
      5 Herd-date: proportion of total variance explained by the herd.
      6 RMSE = root mean square error.
      * P < 0.05
      ** P < 0.01
      *** P < 0.001.

      Lactation Pattern According to Wilmink's Parametric Model

      We used Wilmink's model to fit the patterns of daily milk yield and those of all the nitrogenous fractions of milk with similar pattern during lactation. The model was developed by
      • Wilmink J.B.M.
      Adjustment of test-day milk, fat and protein yield for age, season and stage of lactation.
      to describe the curve of milk production during lactation as follows:
      yt = a + b × exp(−k × t) + c × t,


      where yt is the milk production in time t, expressed in DIM, and exp is the natural exponential function. This curve (depicted in Figure 1, zenith curve, green area under the solid line) is described by 4 parameters, each the result of the sum of 2 components and characterized by a specific variation pattern as follows: the first component (parameter a of the equation) is a long-term fraction (in the example in Figure 1, we assumed a value of 100), which is interpreted as the potential maximum daily yield of the cow at the beginning of lactation; this a fraction is assumed to decrease linearly during lactation with a negative slope (parameter c of the equation), which is multiplied by the time interval t from parturition (corresponding to the DIM) to the drying off of the cow (in our example we assumed c to be −0.10, and lactation length to be 300 d); c is interpreted as representing the persistency of lactation; the second component (parameter b of the equation) is a short-term negative fraction (in the example in Figure 1, we assumed a value of 20), which is interpreted as the adaptation fraction after parturition; the b fraction decreases with time t at a negative exponential rate (parameter k of the equation), which is interpreted as the speed of adaptation of milk production after parturition.
      Figure thumbnail gr1
      Figure 1Example of shapes of lactation curves. All curves have a +100 value for initial long-term fraction (a); a 0.10 daily persistency variation, positive or negative (c); a 20 value for initial short-term fraction, positive or negative (b); and a +0.05 value for adaptation speed (k). The signs for zenith curves are −b and −c, nadir curves are +b and +c, upward curves are −b and +c, and downward curves are +b and −c.
      In the example in Figure 1 (zenith curve), the variation in the first component over time (dependent on parameters a and c) is represented by the area under the dotted line. The variation in the second component (dependent on parameters b and k) is represented by the yellow area between the dotted and solid lines. The resulting variation in actual milk yield throughout lactation is obtained from the difference between these 2 areas, depicted by the evolution of the green area under the solid line.
      Using the Wilmink equation, we calculated the following 6 traits that could be useful in summarizing the lactation pattern:
      • (1)
        the initial value given by the sum of the initial potential value of the long-term fraction and the short-term fraction of adaptation at t = 0 [a + b; in Figure 1, the zenith curve, the initial value is 100 + (−20) = 80];
      • (2)
        the adaptation fraction representing the short component expressed as a percentage of the long-term fraction (b/a × 100; in Figure 1: 20/100 × 100 = 20%);
      • (3)
        the peak value, which is the maximum actual production reached during lactation (Zenith value; 93 in Figure 1);
      • (4)
        the peak DIM (i.e., the time in days from parturition to reach the peak value; 46 DIM in Figure 1);
      • (5)
        the final value, which is the level of production at the end of lactation [ac × 300 DIM; in Figure 1: 100 + (−0.10) × 300 DIM = 70];
      • (6)
        the persistency; value calculated as the percentage drop in production at the end of lactation compared with the initial potential value [(c × 300)/a × 100; in Figure 1: (−0.10 × 300)/100 × 100 = −30%].

      Using Wilmink's Model to Represent Different Lactation Patterns

      The lactation patterns of several milk protein profile traits were very different from those usually observed for daily milk yield. We tested Wilmink's model for all the traits studied by simply allowing the b and c parameters to be either negative or positive, whereas a was always positive and k was always negative.
      The combinations of different signs for the b and c parameters yielded the following 4 different curves (Figure 1):
      • (1)
        “zenith curve,” the classic curve of daily milk yield described above, characterized by a negative b and a negative c parameter;
      • (2)
        “nadir curve,” an inverse pattern with respect to daily milk yield, common for milk fat and protein concentrations, characterized by a positive b and a positive c parameter, resulting in an initial decreasing phase until a minimum value (nadir value) is reached, followed by an increasing phase up to the end of lactation;
      • (3)
        “upward curve,” characterized by a negative b and a positive c parameter, resulting in increasing values throughout lactation (more rapid at the beginning);
      • (4)
        “downward curve,” characterized by a positive b and negative c parameter, resulting in steadily decreasing values throughout the whole of lactation (more rapid at the beginning).

      Statistical Analysis

      It is worth noting that to reduce the complexity of the computations and for the linear modeling procedures in fitting daily milk yield,
      • Wilmink J.B.M.
      Adjustment of test-day milk, fat and protein yield for age, season and stage of lactation.
      and other authors (
      • Heck J.M.L.
      • Schennink A.
      • Van Valenberg H.J.F.
      • Bovenhuis H.
      • Visker M.H.P.W.
      • Van Arendonk J.A.M.
      • Van Hooijdonk A.C.M.
      Effects of milk protein variants on the protein composition of bovine milk.
      ;
      • Rutten M.J.M.
      • Bovenhuis H.
      • Heck J.M.L.
      • van Arendonk J.A.M.
      Predicting bovine milk protein composition based on Fourier transform infrared spectra.
      ) assumed a priori the k parameter of the Wilmink equation to be equal to 0.05. In our case, in adapting the shape of the Wilmink curve to the very different patterns exhibited by the milk protein and NPN fractions expressed in the 3 different ways, we made no a priori assumption regarding the value of the k parameter, and then followed a nonlinear modeling procedure to estimate all model parameters. The data from each milk trait were analyzed with the NLMIXED procedure (SAS Institute Inc.) using the general model:
      yijkl = a + b × exp(–k × t) + c × t + breedi + parityj + herd-datek + eijkl,


      where yijkl represents observations for each trait; a, b, c, and k are the parameters of the lactation curve described above; breedi is the fixed effect of the ith breed of the cow (4 classes: Holstein-Friesian, Brown Swiss, Simmental, local breeds); parityj is the fixed effect of the jth class of parity order (k = 1, 2, and ≥3); herd-datek is the random effect of the kth herd-date (n = 41), considered normally distributed; and eijkl is the residual random error, considered normally distributed. Other models were tested including the interactions of the breed and parity effects with the parameters of the lactation curve, but most of the studied traits did not reach the convergence. For this reason, the interactions were not included.
      Data with residual values outside the range of the mean ± 3 residual standard deviations calculated for each trait according to the previously illustrated model were excluded from the final data set. The incidence of outliers varied according to the trait, ranging from 3.7% of all data for NPN (in % N) to 12.6% for whey protein (in g/d). The distribution of outlier values during lactation was checked considering twelve 30-d classes. All traits presented outlier values distributed in at least 11 of the 12 lactation stages considered. Following the distribution of daily milk yield outliers, the majority of protein fraction traits presented the highest incidence of outliers values in the third class, and the lowest incidence in the tenth.

      RESULTS AND DISCUSSION

      Descriptive Statistics and ANOVA Results

      Descriptive statistics and the results of the ANOVA of the nonlinear mixed model for the detailed milk protein profile reported in qualitative (% N) and quantitative (g/L milk) terms, and in terms of daily production (g/d), as well as daily milk yield are presented in Table 1. The effects of the cow's breed and parity on the detailed milk protein and NPN profiles, and the variability due to the herd's production intensity level and to individual herds are illustrated and discussed in a previous study (
      • Amalfitano N.
      • Stocco G.
      • Maurmayr A.
      • Pegolo S.
      • Cecchinato A.
      • Bittante G.
      Quantitative and qualitative detailed milk protein profiles of 6 cattle breeds: Sources of variation and contribution of protein genetic variants.
      ). Here, these factors were considered possible nuisance sources to be taken into account in the statistical model to eliminate bias in the results of the milk protein and NPN profiles, which were obtained from a variety of breeds and dairy systems (
      • Jõudu I.
      • Henno M.
      • Kaart T.
      • Püssa T.
      • Kärt O.
      The effect of milk protein contents on the rennet coagulation properties of milk from individual dairy cows.
      ;
      • Gustavsson F.
      • Buitenhuis A.J.
      • Johansson M.
      • Bertelsen H.P.
      • Glantz M.
      • Poulsen N.A.
      • Lindmark Månsson H.
      • Stålhammar H.
      • Larsen L.B.
      • Bendixen C.
      • Paulsson M.
      • Andrén A.
      Effects of breed and casein genetic variants on protein profile in milk from Swedish Red, Danish Holstein, and Danish Jersey cows.
      ;
      • Maurmayr A.
      • Pegolo S.
      • Malchiodi F.
      • Bittante G.
      • Cecchinato A.
      Milk protein composition in purebred Holsteins and in first/second-generation crossbred cows from Swedish Red, Montbeliarde and Brown Swiss bulls.
      ).
      Briefly, the true protein content (33.3 g/L) by HPLC analysis represented 93% of the total protein content by FTIR analysis (35.6 g/L), which means a remaining 7% of NPN fractions (2.4 g/L) in the milk. Among the CN fractions, β-CN was the most abundant (29.8 % N; 10.6 g/L; 262 g/d), followed by αS1-CN (26.2 % N; 9.3 g/L; 233 g/d), κ-CN (15.0 % N; 5.4 g/L; 132 g/d), and αS2-CN (7.8 % N; 2.8 g/L; 68 g/d). Of the whey proteins, β-LG was more abundant (12.8% N; 4.5 g/L; 111 g/d) than α-LA (2.1 % N; 0.8 g/L; 19 g/d). Finally, the NPN compounds were divided among minor NPN compounds (4.7 % N; 1.7 g/L; 39 g/d) and milk urea (2.0 % N; 0.3 g/L; 6 g/d).
      The majority of published studies report only qualitative (% N) or quantitative (g/L milk) profiles of CN fractions and whey proteins of only 1 breed (
      • Heck J.M.L.
      • Schennink A.
      • Van Valenberg H.J.F.
      • Bovenhuis H.
      • Visker M.H.P.W.
      • Van Arendonk J.A.M.
      • Van Hooijdonk A.C.M.
      Effects of milk protein variants on the protein composition of bovine milk.
      ;
      • Schopen G.C.B.
      • Heck J.M.L.
      • Bovenhuis H.
      • Visker M.H.P.W.
      • Van Valenberg H.J.F.
      • Van Arendonk J.A.M.
      Genetic parameters for major milk proteins in Dutch Holstein-Friesians.
      ;
      • Bonfatti V.
      • Di Martino G.
      • Cecchinato A.
      • Vicario D.
      • Carnier P.
      Effects of β-κ-casein (CSN2–CSN3) haplotypes and β-lactoglobulin (BLG) genotypes on milk production traits and detailed protein composition of individual milk of Simmental cows.
      ) or just a few breeds (
      • Hallén E.
      • Wedholm A.
      • Andrén A.
      • Lundén A.
      Effect of β-casein, κ-casein and β-lactoglobulin genotypes on concentration of milk protein variants.
      ;
      • Jensen H.B.
      • Holland J.W.
      • Poulsen N.A.
      • Larsen L.B.
      Milk protein genetic variants and isoforms identified in bovine milk representing extremes in coagulation properties.
      ;
      • Gustavsson F.
      • Buitenhuis A.J.
      • Johansson M.
      • Bertelsen H.P.
      • Glantz M.
      • Poulsen N.A.
      • Lindmark Månsson H.
      • Stålhammar H.
      • Larsen L.B.
      • Bendixen C.
      • Paulsson M.
      • Andrén A.
      Effects of breed and casein genetic variants on protein profile in milk from Swedish Red, Danish Holstein, and Danish Jersey cows.
      ), and they are discussed in our previous study on this same data set (
      • Amalfitano N.
      • Stocco G.
      • Maurmayr A.
      • Pegolo S.
      • Cecchinato A.
      • Bittante G.
      Quantitative and qualitative detailed milk protein profiles of 6 cattle breeds: Sources of variation and contribution of protein genetic variants.
      ). The milk protein fraction profiles expressed in all 3 ways are, in any case, similar to those obtained from 1,264 Brown Swiss cows reared in 85 herds representative of different dairy farming systems (
      • Mota L.F.M.
      • Pegolo S.
      • Bisutti V.
      • Bittante G.
      • Cecchinato A.
      Genomic analysis of milk protein fractions in Brown Swiss cattle.
      ). The only appreciable difference was that the glycosylated form of κ-CN was not accounted for, and the average β-LG content was lower (78 vs. 111 g/d) in the Brown Swiss study.
      Regarding the herd-date effect, when the protein fractions were expressed qualitatively (% N), the proportion of phenotypic variance explained by the herd-date was in general low (6–12%), with the exception in the very high herd-date effect on milk urea (82%), confirming the importance of the feeding strategies to influence this trait, and in the moderate herd-date effect on CN (26%), β-CN (17%), and α-LA (32%), possibly due to the incidence of subclinical mastitis (
      • Caffin J.P.
      • Poutrel B.
      • Rainard P.
      Physiological and pathological factors influencing bovine α-Lactalbumin and β-Lactoglobulin concentrations in milk.
      ;
      • Schmitz S.
      • Pfaffl M.W.
      • Meyer H.H.D.
      • Bruckmaier R.M.
      Short-term changes of mRNA expression of various inflammatory factors and milk proteins in mammary tissue during LPS-induced mastitis.
      ;
      • Swanson K.M.
      • Stelwagen K.
      • Dobson J.
      • Henderson H.V.
      • Davis S.R.
      • Farr V.C.
      • Singh K.
      Transcriptome profiling of Streptococcus uberis-induced mastitis reveals fundamental differences between immune gene expression in the mammary gland and in a primary cell culture model.
      ). When we moved from the qualitative traits to the quantitative traits, the effect of the herd became more important because the qualitative traits were multiplied by the total protein content of milk, and thus were more affected by feeding and herd management (37%). This is even more true when observing the daily productions obtained by multiplication of the quantitative traits by the daily milk yield, which was affected by a very large herd-date effect (62%), as expected.
      Regarding the variability in the data, it is also worth noting that the coefficients of variation of all the milk nitrogenous fractions tended to increase from the qualitative to the quantitative profile and to daily production. This is because the nitrogenous fractions (expressed as % N) were multiplied by the total protein content in milk to obtain their contents in milk (g/L), and then by milk yield to obtain their daily production levels (in g/d), thus incorporating the variability in the total protein content of milk and the daily milk yield per cow, as seen with the herd-date effect.

      Modeling the Lactation Curves of Milk Nitrogenous Fractions

      Due to the scarcity of studies modeling milk nitrogenous fractions during lactation, the lactation curve modeling parameters reported here cannot be directly compared with previous literature. All 14 traits representing the qualitative nitrogenous compound profile (Table 2), 14 of the 15 traits representing the quantitative nitrogenous compound profile (Table 3), and all 16 traits representing the daily production of milk and production of the qualitative nitrogenous compound profile (Table 4) presented significant variations throughout lactation. Only the variation in the urea content of milk failed to reach significance level with the model applied. The main reason could be attributed to the fact that urea content is mainly affected by relationships between protein and energy of the diet (
      • Souza V.C.
      • White R.R.
      Variation in urea kinetics associated with ruminant species, dietary characteristics, and ruminal fermentation: A meta-analysis.
      ). In fact, the herd-date represented 74 to 85% of the total variance for this fraction (Table 1), according to the way of expression, and thus the variability due to animal's factors was modest. Moreover, the residuals of the model, plotted against class of DIM (Supplemental Figure S1 for traits expressed as % N, Supplemental Figure S2 for traits expressed as g/L, and Supplemental Figure S3 for traits expressed as g/d; 10.6084/m9.figshare.14977011) revealed a good fitting of the model along the lactation for all traits. Generally, daily milk yield shows greater residual variability in early lactation with respect to mid and late lactation in models accounting for the stage of lactation. This could be due mainly to 2 causes as follows: (1) the greater level of production in early lactation, which can halve from lactation peak to late lactation; and (2) the statistical model adopted. In the case of the use of lactation classes, residual standard deviation increases in early lactation also because of the difficulty to fit a trait that is changing very rapidly in few weeks. With 1-mo classes, it is possible to find cows after 1 or after 4 wk of lactation in the first class. Their daily milk yield is very different, but because they are within the same class of lactation stage, their differences fall entirely into the residual variance. This was not observed in the present study and it can be explained by the fact that a nonlinear regression model can better follow the physiological changes of daily production day by day, thus reducing the residual variability, especially in very early lactation stages characterized by sudden variations. This was also true for nitrogen fractions having a large short-term component (b parameter) or a fast adaptation speed (k parameter).
      Table 2Lactation curve modeling parameters and relative significance (P-values), shape, and initial, final, and peak values of milk nitrogen fractions expressed in percentage of total protein (% N)
      Trait
      Glyco-κ-CN = glycosylated κ-CN; CF-κ-CN = carbohydrate-free κ-CN; CN = Σ CN (αS1-CN + αS2-CN + β-CN + κ-CN); whey proteins = Σ whey proteins (β-LG + α-LA); κ-CN = glyco-κ-CN + CF-κ-CN; NPN compounds = NPN compounds (Σ urea + minor NPN compounds). Examples of minor NPN compounds include small peptides, ammonia, creatine, creatinine.
      Lactation curve equation parameterShape of lactation curveValue during lactationShort/long- term fraction, %Persistency, %
      Long- term fractionPersistencyShort- term fractionAdaptation speedInitial valuePeak DIM (zenithor nadir)Peak value (zenith or nadir)Final value (300 DIM)
      acbka + bac × 300ba×100c×300a×100
      Protein composition, % N
       True protein93.81
      P < 0.001.
      −0.002−11.620.22
      P < 0.05
      Zenith82.193493.7593.35−120
      CN79.16
      P < 0.001.
      −0.003
      P < 0.001.
      −4.89
      P < 0.001.
      0.13
      P < 0.001.
      Zenith74.274179.0278.30−6−1
      αS1-CN26.31
      P < 0.001.
      −0.0011.04
      P < 0.01
      0.03Downward27.3526.044−1
      αS2-CN7.92
      P < 0.001.
      −0.0010.510.02Downward8.437.506−5
      β-CN31.13
      P < 0.001.
      −0.005
      P < 0.01
      −3.16
      P < 0.001.
      0.04Zenith27.979030.5929.76−10−4
      κ-CN14.00
      P < 0.001.
      0.003
      P < 0.001.
      −0.540.18Upward13.4614.96−47
      Glyco-κ-CN3.23
      P < 0.05
      0.010
      P < 0.05
      1.780.01Nadir5.01644.756.195589
      CF-κ-CN11.07
      P < 0.001.
      −0.008−2.100.01Zenith8.97969.448.53−19−22
      Whey proteins14.39
      P < 0.001.
      0.0020.540.01Nadir14.938114.8215.1345
      β-LG11.81
      P < 0.001.
      0.0040.530.01Nadir12.343412.3112.98410
      α-LA2.38
      P < 0.001.
      −0.001
      P < 0.001.
      −0.370.15Zenith2.01252.331.94−15−18
       NPN compounds6.49
      P < 0.001.
      0.00110.590.22
      P < 0.05
      Nadir17.07336.546.891636
      Milk urea2.58
      P < 0.001.
      −0.001
      P < 0.001.
      −0.58
      P < 0.001.
      0.03Zenith1.99842.402.13−23−17
      Minor NPN compounds4.02
      P < 0.001.
      0.003
      P < 0.01
      11.930.22
      P < 0.05
      Nadir15.95324.114.7829719
      1 Glyco-κ-CN = glycosylated κ-CN; CF-κ-CN = carbohydrate-free κ-CN; CN = Σ CN (αS1-CN + αS2-CN + β-CN + κ-CN); whey proteins = Σ whey proteins (β-LG + α-LA); κ-CN = glyco-κ-CN + CF-κ-CN; NPN compounds = NPN compounds (Σ urea + minor NPN compounds). Examples of minor NPN compounds include small peptides, ammonia, creatine, creatinine.
      * P < 0.05
      ** P < 0.01
      *** P < 0.001.
      Table 3Lactation curve modeling parameters and relative significance (P-values), shape, and initial, final, and peak values of milk nitrogen fractions expressed in grams per liter of milk
      Trait
      Glyco-κ-CN = glycosylated κ-CN; CF-κ-CN = carbohydrate-free κ-CN; CN = Σ CN (αS1-CN + αS2-CN + β-CN+ κ-CN); whey proteins = Σ whey proteins (β-LG + α-LA); κ-CN = glyco-κ-CN + CF-κ-CN; NPN compounds = NPN compounds (Σ urea + minor NPN compounds). Examples of minor NPN compounds include small peptides, ammonia, creatine, creatinine.
      Lactation curve equation parameterShape of lactation curveValue during lactationShort/long- term fraction, %Persistency, %
      Long- term fractionPersistencyShort- term fractionAdaptation speedInitial valuePeak DIM (zenith or nadir)Peak value (zenith or nadir)Final value (300 DIM)
      acbka + bac × 300ba×100c×300a×100
      Milk composition, g/L
       Total protein30.60
      P < 0.001.
      0.026
      P < 0.001.
      19.85
      P < 0.001.
      0.15
      P < 0.001.
      Nadir50.453231.6138.506526
      True protein28.67
      P < 0.001.
      0.024
      P < 0.001.
      14.87
      P < 0.001.
      0.14
      P < 0.001.
      Nadir43.543229.6235.985226
      CN24.21
      P < 0.001.
      0.020
      P < 0.001.
      12.09
      P < 0.001.
      0.14
      P < 0.001.
      Nadir36.303124.9730.135024
      αS1-CN8.15
      P < 0.001.
      0.006
      P < 0.001.
      5.60
      P < 0.001.
      0.14
      P < 0.001.
      Nadir13.75348.4010.006923
      αS2-CN2.55
      P < 0.001.
      0.001
      P < 0.001.
      2.590.23Nadir5.14282.582.8210211
      β-CN9.42
      P < 0.001.
      0.007
      P < 0.001.
      3.30
      P < 0.05
      0.17
      P < 0.01
      Nadir12.72269.6411.453522
      κ-CN4.31
      P < 0.001.
      0.005
      P < 0.001.
      1.72
      P < 0.05
      0.12
      P < 0.05
      Nadir6.02314.495.764034
      Glyco-κ-CN0.820.005
      P < 0.05
      0.740.01Nadir1.56341.542.3390178
      CF-κ-CN3.04
      P < 0.001.
      0.001
      P < 0.001.
      0.340.12Nadir3.38303.083.361111
      Whey proteins4.47
      P < 0.001.
      0.004
      P < 0.001.
      1.63
      P < 0.05
      0.10
      P < 0.05
      Nadir6.10374.675.763629
      β-LG3.66
      P < 0.001.
      0.004
      P < 0.001.
      1.57
      P < 0.05
      0.10
      P < 0.05
      Nadir5.23373.864.934335
      α-LA0.74
      P < 0.001.
      0.0000.430.23Nadir1.17390.740.75581
      NPN compounds2.04
      P < 0.001.
      0.002
      P < 0.001.
      5.670.24
      P < 0.05
      Nadir7.70272.102.6127828
      Milk urea0.34−0.000−0.100.01Zenith0.251940.270.27−28−16
      Minor NPN compounds1.27
      P < 0.001.
      0.002
      P < 0.001.
      6.080.25
      P < 0.05
      Nadir7.35271.321.8148143
      1 Glyco-κ-CN = glycosylated κ-CN; CF-κ-CN = carbohydrate-free κ-CN; CN = Σ CN (αS1-CN + αS2-CN + β-CN+ κ-CN); whey proteins = Σ whey proteins (β-LG + α-LA); κ-CN = glyco-κ-CN + CF-κ-CN; NPN compounds = NPN compounds (Σ urea + minor NPN compounds). Examples of minor NPN compounds include small peptides, ammonia, creatine, creatinine.
      * P < 0.05
      ** P < 0.01
      *** P < 0.001.
      Table 4Lactation curve modeling parameters and relative significance (P-values), shape, and initial, final, and peak values of milk yield and daily production of milk nitrogen fractions expressed in grams per day per cow
      Trait
      Glyco-κ-CN = glycosylated κ-CN; CF-κ-CN = carbohydrate-free κ-CN; CN = Σ CN (αS1-CN + αS2-CN + β-CN + κ-CN); whey proteins = Σ whey proteins (β-LG + α-LA); κ-CN = glyco-κ-CN + CF-κ-CN; NPN compounds = NPN compounds (Σ urea + minor NPN compounds). Examples of minor NPN compounds include small peptides, ammonia, creatine, creatinine.
      Lactation curve equation parameterShape of lactation curveValues during lactationShort/long- term fraction, %Persistency, %
      Long- term fractionPersistencyShort- term fractionAdaptation speedInitial valuePeak DIM (zenith or nadir)Peak value (zenith or nadir)Final value (300 DIM)
      acbka+ bac × 300ba×100c×300a×100
      Milk yield, kg/d30.41
      P < 0.001.
      −0.040
      P < 0.001.
      −3.74
      P < 0.001.
      0.02
      P < 0.05
      Zenith26.683227.2218.32−12−40
      Protein production, g/d
       Total protein917
      P < 0.001.
      −0.659
      P < 0.001.
      4960.20
      P < 0.05
      Downward141371954−22
      True protein851
      P < 0.001.
      −0.612
      P < 0.001.
      5660.28Downward141766867−22
      CN723
      P < 0.001.
      −0.533
      P < 0.001.
      3150.21Downward103856344−22
      αS1-CN242
      P < 0.001.
      −0.182
      P < 0.001.
      183
      P < 0.05
      0.18
      P < 0.001.
      Downward42418776−23
      αS2-CN73
      P < 0.001.
      −0.075
      P < 0.001.
      1060.23
      P < 0.01
      Downward17951144−31
      β-CN277
      P < 0.001.
      −0.216
      P < 0.001.
      1150.27Downward39221242−23
      κ-CN127
      P < 0.001.
      −0.067
      P < 0.001.
      590.20Downward18610746−16
      Glyco-κ-CN39
      P < 0.001.
      0.014
      P < 0.05
      1250.29
      P < 0.01
      Nadir16427394332311
      CF-κ-CN103
      P < 0.001.
      −0.130
      P < 0.05
      −230.01Zenith80458262−22−38
      Whey proteins133
      P < 0.001.
      −0.092
      P < 0.001.
      680.21Downward20110551−21
      β-LG109
      P < 0.001.
      −0.063
      P < 0.001.
      980.19
      P < 0.05
      Downward2079189−17
      α-LA21
      P < 0.001.
      −0.025
      P < 0.001.
      170.26Downward381482−36
      NPN compounds58
      P < 0.001.
      −0.038
      P < 0.001.
      2320.28
      P < 0.01
      Downward29047401−20
      Milk urea9
      P < 0.001.
      −0.012
      P < 0.01
      −30.01Zenith67565−30−42
      Minor NPN compounds37
      P < 0.001.
      −0.0162570.29
      P < 0.001.
      Downward29432703−13
      1 Glyco-κ-CN = glycosylated κ-CN; CF-κ-CN = carbohydrate-free κ-CN; CN = Σ CN (αS1-CN + αS2-CN + β-CN + κ-CN); whey proteins = Σ whey proteins (β-LG + α-LA); κ-CN = glyco-κ-CN + CF-κ-CN; NPN compounds = NPN compounds (Σ urea + minor NPN compounds). Examples of minor NPN compounds include small peptides, ammonia, creatine, creatinine.
      * P < 0.05
      ** P < 0.01
      *** P < 0.001.
      If Wilmink's model is well known for describing the lactation pattern of daily milk yield (
      • Silvestre A.M.
      • Martins A.M.
      • Santos V.A.
      • Ginja M.M.
      • Colaço J.A.
      Lactation curves for milk, fat and protein in dairy cows: A full approach.
      ;
      • Macciotta N.P.P.
      • Dimauro C.
      • Rassu S.P.G.
      • Steri R.
      • Pulina G.
      The mathematical description of lactation curves in dairy cattle.
      ;
      • Bouallegue M.
      • M'Hamdi N.
      Mathematical modeling of lactation curves: A review of parametric model.
      ), modifying it to allow for both positive and negative signs for persistency c parameter) of the long-term milk fraction, and for the entity of the short-term (adaptation) fraction (b parameter) gives a good representation of the lactation pattern of all the qualitative and quantitative nitrogenous milk profiles and also their daily production levels. Combinations of the signs of the c and b parameters allowed us to identify 4 possible lactation curve shapes, as illustrated in Figure 1 and described in Materials and Methods.
      • Macciotta N.P.P.
      • Vicario D.
      • Cappio-Borlino A.
      Detection of different shapes of lactation curve for milk yield in dairy cattle by empirical mathematical models.
      designated these as “abnormal” lactation curves for the daily milk yield of individual cows, whereas here they appeared to be normal curves of the nitrogenous fractions throughout lactation.
      These 4 types of lactation curve were identified with different frequencies in the patterns of the milk nitrogenous fractions, but they seemed to be associated mostly with how the fractions were expressed. The zenith curve type, as expected, was a good fit for the daily milk yield data from the 1,342 cows sampled for this study (Figure 2). It is worth noting that in a preliminary statistical analysis carried out with the same model, where the lactation effect was not described by a parametric model but by a DIM factor with 11 monthly classes, the root mean square error of the daily milk yield was ± 4.6 kg/d versus ± 4.5 kg/d in this model.
      Figure thumbnail gr2
      Figure 2The typical Wilmink's curve (zenith type) fitted to represent the daily milk yield (kg/d) along lactation of the 1,342 cows sampled in the present study (long-term fraction = 30.41 kg/d; short-term fraction = −3.74 kg/d; persistency = −0.040 kg/d; adaptation speed = 0.02).
      Comparing the results of the present study (Figure 2) with those of
      • Wilmink J.B.M.
      Adjustment of test-day milk, fat and protein yield for age, season and stage of lactation.
      original study, we found the long-term component of daily milk production to be, on average, higher (a = 30.41 kg/d) in our multibreed herds than in
      • Wilmink J.B.M.
      Adjustment of test-day milk, fat and protein yield for age, season and stage of lactation.
      Dutch Friesian herds (a = 25.6 kg/d), reflecting improvements in genetics and management in the dairy sector in the last 3 decades or so. Moreover, the linear decrease in the long-term milk production component (persistency) was less steep in our study than in the
      • Wilmink J.B.M.
      Adjustment of test-day milk, fat and protein yield for age, season and stage of lactation.
      study (−0.04 vs. −0.06 kg/d).
      We calculated a much slower speed of adaptation (k = 0.02) of milk production after calving than the a priori assumption of k = 0.05 in
      • Wilmink J.B.M.
      Adjustment of test-day milk, fat and protein yield for age, season and stage of lactation.
      . Conversely, the short-term component of milk yield (adaptation fraction) was greater in Wilmink's study (1987) than in ours (−5.60 vs. −3.74 kg/d, respectively). The combination of the differences in these 2 parameters resulted in a less steep ascending phase of lactation in our study, but a similar interval from calving to peak (peak DIM of about 1 mo in both studies).
      It is worth noting that we found a significant interaction between the effects of lactation stage and the level of production intensity in the herd in our previous study (
      • Stocco G.
      • Cipolat-Gotet C.
      • Bobbo T.
      • Cecchinato A.
      • Bittante G.
      Breed of cow and herd productivity affect milk composition and modeling of coagulation, curd firming, and syneresis.
      ). The cows in high-production herds (mainly modern farms with loose animals, milking parlors, and total mixed diets) exhibited a longer phase of increasing daily milk yield (i.e., a slower adaptation of milk production) and a delayed peak of production. In contrast, the cows reared in low-production herds (mainly traditional farms, sometimes with tied cows, and feed based on hay and some compound feed) exhibited almost no peak of lactation (no adaptation period). Farming system, breed, and parity of cows are known to modify the curve of lactation in terms of daily milk yield (
      • Wood P.D.P.
      Breed variations in the shape of the lactation curve of cattle and their implications for efficiency.
      ;
      • García S.C.
      • Holmes C.W.
      Lactation curves of autumn- and spring-calved cows in pasture-based dairy systems.
      ;
      • Macciotta N.P.P.
      • Dimauro C.
      • Rassu S.P.G.
      • Steri R.
      • Pulina G.
      The mathematical description of lactation curves in dairy cattle.
      ), and possibly they could affect also the pattern of nitrogenous compounds during lactation. The interactions of breed and parity effects with the parameters of the lactation curve were tested, but the model failed to converge for many traits. The major problems seemed to be caused by the estimation of the short-term (adaptation) fraction (b parameter) and its speed (k parameter). These 2 parameters depend almost exclusively on differences observed in the first 1 to 2 mo of lactation. This means that only one-tenth or one-fifth of the cows included in the data set contributed to the estimation of 7 coefficients per parameter (interactions with 4 classes of breed and 3 classes of parity). For this reason, breed and parity effects were included in the statistical model merely to ensure the average results of the lactation curves of the detailed milk protein and NPN compound profiles were unbiased. However, it is evident that this model of lactation curve could be a very useful tool to analyze the effect of farming system, breed, and parity (and interactions) on long- and short-term lactation fractions, speed of adaptation after calving, and persistency during lactation, provided that the number of animals sampled and the experimental design are adequate.

      Effects of the Ways in Which the Nitrogenous Fractions During Lactation Are Expressed

      In the case of the qualitative milk nitrogenous compound profile, as the sum of the proportions of all the compounds is a constant (100%), the decreasing patterns of some substances should be compensated for by increasing patterns of others, and the presence of positive peaks for some substances should be accompanied by the presence of negative peaks for others. Therefore, it is not surprising that the shapes of the curves of the nitrogenous substances expressed as percent N are highly variable. It can be seen from Table 2 that all 4 types of curves were represented by these traits: 6 traits followed the zenith type curve, 5 traits followed the inverse pattern of the nadir type curve, one followed the upward type curve, and 2 followed the downward type curve. It should be kept in mind that these traits represent the evolution during lactation of the relative emphasis, from a physiological point of view, of the secretion in milk of different nitrogenous substances.
      The modeling of the same milk nitrogenous substances is completely different when expressed quantitatively in grams per liter of milk. In this case, the sum of all the nitrogenous fractions is not a constant, but is equal to the total protein content of milk (i.e., the content of N multiplied by 6.38). Therefore, the pattern of the nitrogenous compounds that is quantitatively expressed, or at least those nitrogenous fractions in greater quantities, is expected to be not much different from that of the total sum of all the nitrogenous compounds (i.e., the total protein content in g/L of milk). It can be seen from Table 3 that the total protein content of milk and almost all of the nitrogenous traits expressed in grams per liter are indeed characterized by a nadir type of lactation curve. Only the lesser abundant nitrogenous component in milk (urea) showed a zenith type curve. These results are consistent with those of
      • Ng-Kwai-Hang K.F.
      • Hayes J.F.
      • Moxley J.E.
      • Monardes H.G.
      Variation in milk protein concentrations associated with genetic polymorphism and environmental factors.
      , who found that all the major milk protein fractions were characterized by a nadir-like pattern during lactation. From a physiological point of view, it is known that the increase of milk yield after calving is accompanied by a decrease in fat and protein content of milk. This is due to a “dilution” effect. The retention of water in milk is strictly related to lactose secretion (in fact, its variability is very modest). After calving, lactose secretion increases at a more rapid rate than protein and fat, and this causes their decrease when expressed in grams per liter (
      • Fox P.F.
      • Uniacke-Lowe T.
      • McSweeney P.L.H.
      • O'Mahony J.A.
      Lactose.
      ). After the peak of lactation, the decrease in daily lactose secretion is the main determinant of the reduction of water secretion, and then of daily milk yield, whereas the decrease of daily fat and protein secretion is slower, causing their increase in milk content (g/L) until the end of lactation.
      In the case of daily production of the different nitrogenous compounds, their patterns throughout lactation were expected to be influenced mainly by the lactation curve of daily milk protein yield. This, in turn, was calculated by multiplying the daily milk yield (kg/d), which showed a zenith curve (Figure 2), by the milk protein content (g/L), which showed a nadir curve (Table 3). The result was a downward curve with a decreasing daily protein production (g/d) throughout the entire lactation (Table 4). As expected, almost all the milk nitrogen compounds (12 out of 15 traits) showed a downward curve throughout lactation, with the exception of CF-κ-CN and milk urea that was characterized by a zenith curve, similar to milk yield, and glyco-κ-CN by a nadir curve, similar to milk protein content (Table 4). The daily production allowed us to represent the total metabolic effort of the cows during lactation.

      Pattern of Milk Proteins Throughout Lactation

      Qualitatively, true protein showed a zenith curve characterized by a very large long-term fraction (a = 93.81% of total N) and very good persistency (the decrease c being only −0.002% N/d), so that the value at the end of lactation was almost equal to that at the beginning (93.35% N; Table 2). Its short-term fraction was also very large (b = −11.62% N), the largest among all the protein fractions, but was accompanied by the fastest adaptation speed (k = 0.22). The combined results of these equation parameters are depicted in Figure 3, which shows the very rapid increase in true protein from the initial value of 82.19% N to the peak value of 93.75% N at 34 d after parturition (Table 2). Even though the first 5 d of lactation were not sampled, and so colostrum composition did not affect the lactation modeling, it seemed clear that the early phases of the different curves captured the transition from colostrum to milk production from a physiological point of view. In particular, the short-term fraction and the adaptation speed were the most informative modeling parameters for this scope. In fact, the Wilmink's equation has also been appreciated for its ability to model the early lactation stage through the k parameter (
      • Bouallegue M.
      • M'Hamdi N.
      Mathematical modeling of lactation curves: A review of parametric model.
      ).
      Figure thumbnail gr3
      Figure 3Pattern of total protein, true protein, and total CN content during lactation expressed in percentage of total protein (% N), grams per liter of milk (g/L), and daily production (g/d).
      Expressed in terms of content in milk, true protein closely paralleled the nadir curve of total protein (Figure 3), with a very rapid decrease from the initial values of 50.45 g/L for total protein and 43.54 g/L for true protein to the minimum points of 31.61 and 29.62 g/L, respectively, about 1 mo after calving, then increasing again to the final values of 38.50 and 35.98 g/L (Table 3). Regarding the parameters of the equation, the positive short-term component represented more than 50% of the long-term component and showed a rapid adaptation speed for both traits. The positive persistency was also large, increasing the long-term component by about a quarter during lactation (Table 3).
      It is well known that, from a biological point of view, during the period from calving to the peak of lactation, the cow cannot cope with the increased nutrient demand for milk production generated by the proliferation of the secretory cells of the mammary glands. The animal finds itself in a situation of negative energy balance due to limited ingestion capacity. Therefore, it is unable to secrete nutrients at the same rate as it produces milk, leading to a dilution effect. In particular, the nadir point of milk protein content coincides with the maximum milk yield, whereas the nadir point of fat is reached 3 wk later (
      • Wood P.D.P.
      Algebraic models of the lactation curves for milk, fat and protein production, with estimates of seasonal variation.
      ;
      • Stanton T.L.
      • Jones L.R.
      • Everett R.W.
      • Kachman S.D.
      Estimating milk, fat, and protein lactation curves with a test day model.
      ;
      • Quinn N.
      • Killen L.
      • Buckley F.
      Empirical algebraic modelling of live weight of Irish dairy cows over lactation.
      ). After the peak, the energy balance changes and production decreases, and thus there is a concentrating effect, which explains the rise in the nutrient concentration after the nadir point (
      • Auldist M.J.
      • Walsh B.J.
      • Thomson N.A.
      Seasonal and lactational influences on bovine milk composition in New Zealand.
      ).
      When expressed as daily production per cow, true protein yield also paralleled total protein production (Figure 3) in a downward curve that decreased very rapidly in the first few weeks, and less rapidly toward the end of lactation. In these cases, the positive short-term components represented more than half of the long-term components, with a very rapid adaptation speed. However, in these cases negative, persistency was −22% of the long-term component of both traits (Table 4).
      Unfortunately, no other studies have applied a lactation curve model to the detailed protein profile expressed as daily production, but
      • Wilmink J.B.M.
      Adjustment of test-day milk, fat and protein yield for age, season and stage of lactation.
      also observed a downward curve for the pattern of protein yield during lactation. This was confirmed by
      • Silvestre A.M.
      • Martins A.M.
      • Santos V.A.
      • Ginja M.M.
      • Colaço J.A.
      Lactation curves for milk, fat and protein in dairy cows: A full approach.
      , who applied Wood's model (
      • Wood P.D.P.
      Algebraic models of the lactation curves for milk, fat and protein production, with estimates of seasonal variation.
      ) instead of Wilmink's. At the same time,
      • Miciński J.
      • Jastrzebski M.
      • Klupczynski J.
      Yield and composition of milk from Polish Holstein-Friesian and Jersey cows in particular months of the first lactations as dependent on milk protein polymorphism (short communication).
      found a decrease in total CN and whey protein yields from the first to the last stage of lactation in a Jersey population, whereas the patterns for the Holstein-Friesians were more similar to a zenith curve.

      Pattern of CN Throughout Lactation

      As expected, given that they account for 84% of true protein, CN showed a zenith curve similar to that of true protein when expressed as percent N (Table 2 and Figure 3), a nadir curve when expressed in grams per liter of milk (Table 3 and Figure 3), and a downward curve when expressed in grams per day (Table 4 and Figure 3). The only appreciable difference between CN and true protein was noted when they were expressed as percent N, as the short-term component of the former was about half that of the latter, but it had a lower adaptation speed; therefore, the peak was not delayed by much.
      The zenith shape of the total CN pattern with its constant decrease in the second phase of lactation is consistent with the decrease in CN number at the end of lactation observed in some studies, which is connected with the increase in plasmin activity (
      • Auldist M.J.
      • Hubble I.B.
      Effects of mastitis on raw milk and dairy products.
      ;
      • Maurmayr A.
      • Pegolo S.
      • Malchiodi F.
      • Bittante G.
      • Cecchinato A.
      Milk protein composition in purebred Holsteins and in first/second-generation crossbred cows from Swedish Red, Montbeliarde and Brown Swiss bulls.
      ). In fact, the epithelial cells in mammary glands are known to increase their permeability toward the end of lactation, leading to the release of more somatic cells in milk that contain plasminogen activators (
      • France T.C.
      • O'Mahony J.A.
      • Kelly A.L.
      The plasmin system in milk and dairy products.
      ). However, the lactation values of the present study reveal that the CN number increased by +4.0% N from d 0 to d 300 of lactation, which is due to a rapid increase of +4.7% N in the very initial stage of lactation caused by the short-term fraction and its adaptation speed, and a gradual decrease of −0.7% N regulated by the persistency of the long-term fraction after the zenith was reached at 41 d postpartum (very close to the milk yield peak at 32 d shown in Table 4). This important, rapid adaptation at the beginning of lactation cannot be captured when the statistical model for analyzing milk nitrogen composition includes large classes of DIM (1 mo or more) or when the beginning of lactation is excluded from the database (
      • Ng-Kwai-Hang K.F.
      • Hayes J.F.
      • Moxley J.E.
      • Monardes H.G.
      Variability of test-day milk production and composition and relation of somatic cell counts with yield and compositional changes of bovine milk.
      ;
      • Ostersen S.
      • Foldager J.
      • Hermansen J.E.
      Effects of stage of lactation, milk protein genotype and body condition at calving on protein composition and renneting properties of bovine milk.
      ;
      • Maurmayr A.
      • Pegolo S.
      • Malchiodi F.
      • Bittante G.
      • Cecchinato A.
      Milk protein composition in purebred Holsteins and in first/second-generation crossbred cows from Swedish Red, Montbeliarde and Brown Swiss bulls.
      ).
      In regards to individual CN, it is not unexpected that β-CN, the most abundant in terms of relative weight, showed a similar pattern to that of total CN, irrespective of how it was expressed (Figure 4). The β-CN is the most sensitive fraction to the proteolytic action of the plasmin, which degrades this CN in γ-CN and other proteose-peptones (
      • France T.C.
      • O'Mahony J.A.
      • Kelly A.L.
      The plasmin system in milk and dairy products.
      ). In contrast, αS1-CN and αS2-CN showed a very different pattern to that of β-CN when expressed as percent N (Figure 4). These were the only nitrogenous compounds showing a downward curve (positive short-term component and negative persistency, no peak). In any case, their short-term fractions were small (4% of the long-term component for αS1-CN, 6% for αS2-CN), their adaptation speed and persistency were slow (Table 2), and thus their curves were almost flat. The αS1-CN and αS2-CN pattern was similar to the nadir curve of CN and β-CN when expressed in grams per liter (Table 3 and Figure 4), and similar to the downward curve of these when expressed in grams per day (Table 4 and Figure 4).
      Figure thumbnail gr4
      Figure 4Pattern of αS1-CN, αS2-CN, and β-CN content during lactation expressed in percentage of total protein (% N), grams per liter of milk (g/L), and daily production (g/d).
      The situation with regard to κ-CN is more complicated because of the high incidence of glycosylation of this protein fraction (about one-third of the total; Table 1). The glyco-κ-CN and CF-κ-CN showed very different patterns throughout lactation, irrespective of how they were expressed. Qualitatively, in percent N, glyco-κ-CN showed a nadir curve, and CF-κ-CN showed a completely different zenith curve (Figure 5). Their sum gave total κ-CN, which followed an upward curve, the only one of all the nitrogenous compounds, however expressed, to do so. The glyco-κ-CN was characterized by a relatively small long-term component (a = 3.23% N) that showed a large positive persistency with an increase of 89% at the end of lactation (Table 2). The short-term component was large (55% of the long-term component) and positive with a slow adaptation speed, which explained the relative flatness of the lactation curve in the first 3 mo and the subsequent tendency to increase significantly until the end of lactation (Figure 5). The unglycosylated form followed an almost reverse pattern, with a relatively large long-term component characterized by negative persistency (−22% at the end of lactation), and a smaller, negative short-term component with a slow adaptation speed. As the short-term fractions of the 2 κ-CN forms had similar values and opposite signs, their sum showed a very small negative short-term component with a rapid speed of adaptation. This caused it to disappear within a few weeks, and the lactation pattern was dominated by a large long-term component increasing slowly during lactation (persistency +7%), explaining the upward curve.
      Figure thumbnail gr5
      Figure 5Pattern of total κ-CN content, and its glycosylated (glyco-κ-CN) and carbohydrate-free (CF-κ-CN) forms, during lactation expressed in percentage of total protein (% N), grams per liter of milk (g/L), and daily production (g/d).
      Regarding the individual CN fractions, we confirmed the increase in β-CN at the expense of αS1-CN and αS2-CN observed by
      • Ostersen S.
      • Foldager J.
      • Hermansen J.E.
      Effects of stage of lactation, milk protein genotype and body condition at calving on protein composition and renneting properties of bovine milk.
      and
      • Maurmayr A.
      • Pegolo S.
      • Malchiodi F.
      • Bittante G.
      • Cecchinato A.
      Milk protein composition in purebred Holsteins and in first/second-generation crossbred cows from Swedish Red, Montbeliarde and Brown Swiss bulls.
      ; the only difference was that in the present study, κ-CN tended to increase from the beginning to the end of the lactation.
      • Kroeker E.M.
      • Ng-Kwai-Hang K.F.
      • Hayes J.F.
      • Moxley J.E.
      Effects of environmental factors and milk protein polymorphism on composition of casein fraction in bovine milk.
      also observed a zenith-like pattern for the β-CN relative concentration and a consequent drop in the αS fraction, whereas κ-CN was more stable throughout the entire lactation. The zenith shape of β-CN, with its slight decrease after the peak at 90 d, is consistent with an increase in plasmin proteolytic activity during the second phase of lactation (
      • Auldist M.J.
      • Hubble I.B.
      Effects of mastitis on raw milk and dairy products.
      ). When expressed in grams per liter, all 3 κ-CN traits showed a nadir curve, but with some evident differences between them (Figure 5). The glyco-κ-CN showed a small long-term component increasing relatively rapidly so that persistency was +178% (Table 3), and an almost similar positive short-term component with a slow adaptation speed. The result was a curve that was almost flat at the beginning, and increased rapidly after the third month of lactation. In contrast, CF-κ-CN showed a large, almost constant long-term component with a small short-term component rapidly disappearing (Table 3). The result of both was a curve with a large long-term component increasing significantly throughout lactation, and a relatively large short-term component that rapidly disappeared. The negative peaks of all the 3 curves were reached about 1 mo after lactation.
      Last, when expressed in grams per day, glyco-κ-CN followed a nadir curve, and CF-κ-CN followed a zenith curve, whereas their sum, unlike κ-CN in percent N, showed a downward curve (Figure 5). In this case, glyco-κ-CN was characterized by an extremely large short-term component (323% of the long-term component) that was quickly eliminated by the high adaptation speed (k = 0.29). This explained the sharp heavy drop over 1 mo in the glycosylated fractions, whereas the small positive persistency (11%) explained the almost flat pattern of the lactation curve after the nadir point. In CF-κ-CN, the short-term component lasted longer, but was smaller (−22% of the long-term component), reaching a zenith peak almost equal to the initial value of the curve. The negative effect of persistency on the long-term component was greater than that of the glycosylated fraction (−38% at the end of lactation; Table 4). This explained why we observed a downward curve for total κ-CN daily yield, with a sharp drop in the first phase of lactation and a more gradual decrease in the second, instead of an upward curve as with percent N.

      Pattern of Whey Proteins Throughout Lactation

      The 2 major whey proteins also followed reverse patterns throughout lactation when expressed qualitatively as percent N, with a nadir curve observed for β-LG and a zenith curve for α-LA. Both the b and c equation parameters were positive for β-LG and negative for α-LA (Table 2). It is worth noting that the short-term component of α-LA was proportionally greater (4 times) than that of β-LG, but its adaptation speed was also much greater (k = 0.15 vs. k = 0.01, respectively), so that its effect on the lactation curve disappeared very rapidly (Figure 6).
      Figure thumbnail gr6
      Figure 6Pattern of whey proteins content during lactation expressed in percentage of total protein (% N), grams per liter of milk (g/L), and daily production (g/d).
      The 2 major whey proteins followed a nadir curve when expressed as grams per liter of milk (Table 3 and Figure 6), in line with the pattern of total protein in grams per liter, and both followed a downward curve when expressed as daily production (Table 4 and Figure 6), in line with the pattern of daily protein yield. However, the pattern of α-LA was flatter than that of β-LG throughout lactation, which could be explained by its close connection with the lactose content in milk, which fluctuates less than the other principal milk components during lactation, except in the presence of clinical disease (
      • Fox P.F.
      • Uniacke-Lowe T.
      • McSweeney P.L.H.
      • O'Mahony J.A.
      Lactose.
      ). As α-LA represents only about 15% of their sum (Table 1), it is not surprising that the 3 patterns of whey proteins were very similar to those of β-LG (Figure 6), irrespective of how they were expressed.
      For the whey proteins, we confirmed the nadir shape reported by
      • Ostersen S.
      • Foldager J.
      • Hermansen J.E.
      Effects of stage of lactation, milk protein genotype and body condition at calving on protein composition and renneting properties of bovine milk.
      , who found that the nadir point was reached at the same time as the CN zenith, whereas
      • Maurmayr A.
      • Pegolo S.
      • Malchiodi F.
      • Bittante G.
      • Cecchinato A.
      Milk protein composition in purebred Holsteins and in first/second-generation crossbred cows from Swedish Red, Montbeliarde and Brown Swiss bulls.
      , who included large DIM classes in their statistical model, observed an increase only after the nadir, without the initial prior decrease. In the present study, total whey proteins took twice as long as the CN to reach the peak (81 vs. 41 d). However, the peaks of the individual whey proteins were almost the same as that of total CN. The difference between total CN and whey protein in terms of the DIM peak highlighted the fact that, even with a similar short-term to long-term fraction ratio (6% for CN, 4% for whey protein), the 2 fractions reached their peaks at different times because of their very different speeds of adaptation (k = 0.22 for CN vs. k = 0.01 for whey proteins). This showed the potential importance of estimating the k parameter instead of fixing it a priori (
      • Wilmink J.B.M.
      Adjustment of test-day milk, fat and protein yield for age, season and stage of lactation.
      ).

      Pattern of NPN Compounds Throughout Lactation

      In the case of NPN compounds, we also found that milk urea followed a completely different pattern to that characterizing the minor NPN compounds. Also, given that urea represents a small proportion of all NPN compounds, the sum of all NPN compounds followed patterns that were almost parallel to those of the minor NPN compounds (Figure 7).
      Figure thumbnail gr7
      Figure 7Pattern of the nonprotein N fractions (NPN) content during lactation expressed in percentage of total protein (% N), grams per liter of milk (g/L), and daily production (g/d).
      As milk urea accounts for about one-fifth of total NPN compounds, to give a clearer idea of its pattern throughout lactation, the curves are shown without the other NPN traits in Figure 8. Note that milk urea was characterized by zenith curves, regardless of how it was expressed. With all 3 means of expression (Tables 2, 3, and 4), urea exhibited similar short-term components (−23% to −30% of the long-term components) that were characterized by a low adaptation speed (k parameter 0.01–0.03), which led to delayed lactation peaks as follows: between the second and third month of lactation for urea expressed as percent N (Table 2) and in grams per day (Table 3), and at about 6 mo after calving when expressed in grams per liter of milk (Table 4). Persistency of the long-term component was always negative; at the end of lactation, it was −16% for urea expressed as percent N, −17% when expressed in grams per liter, and −42% expressed in grams per day. However, considering that the equation parameters for the urea content of milk (g/L) were not as significant as those for CF-κ-CN, its zenith curve could not be considered different from a flat, constant value throughout lactation. It is worth noting that, in the previous study on the main sources of variation in the same data set, the phenotypic variability in milk urea was shown to be influenced mainly by factors associated with herd management and feeding rather than with the animal, such as parity and DIM (
      • Amalfitano N.
      • Stocco G.
      • Maurmayr A.
      • Pegolo S.
      • Cecchinato A.
      • Bittante G.
      Quantitative and qualitative detailed milk protein profiles of 6 cattle breeds: Sources of variation and contribution of protein genetic variants.
      ).
      Figure thumbnail gr8
      Figure 8Pattern of milk urea content during lactation expressed in percentage of total protein (% N), grams per liter of milk (g/L), and daily production (g/d).
      The patterns of minor (and also of all) NPN compounds showed a completely different pattern (Figure 7). They followed a nadir curve when expressed as percent N and grams per liter, and a downward curve in grams per day (Figure 7). Minor NPN compounds (and also all NPN compounds) were characterized by the greatest proportional short-term components (+297% to +700% of the long-term component) and the highest adaptation speed (k from 0.22–0.29). This resulted in very high values at calving, which decreased sharply during the first couple of weeks, and then tended to stabilize thereafter (Figure 7).
      The NPN measured as percent N by
      • Ostersen S.
      • Foldager J.
      • Hermansen J.E.
      Effects of stage of lactation, milk protein genotype and body condition at calving on protein composition and renneting properties of bovine milk.
      exhibited a more constant decrease from the beginning to the end of lactation than in the present work.
      • Ng-Kwai-Hang K.F.
      • Hayes J.F.
      • Moxley J.E.
      • Monardes H.G.
      Percentages of protein and nonprotein nitrogen with varying fat and somatic cells in bovine milk.
      also observed a nadir curve for NPN compounds expressed in grams per liter, with a less dramatic drop in the first month of lactation (nadir point at ~30 d).

      CONCLUSIONS

      In the present work we found that a nonlinear parametric model was able to successfully describe the pattern throughout lactation of 14 milk nitrogenous compounds expressed qualitatively (% N), quantitatively (g/L milk), and as daily production (g/d per cow). The results show that individual nitrogenous fractions follow very different, complex patterns throughout lactation. In particular, 4 different lactation curve shapes were observed. The zenith shape (the classical lactation curve of milk production), characterized by an initial increase to a zenith and a subsequent decrease, described 10 milk nitrogenous traits, mainly in terms of % N. The nadir shape, an inverse lactation curve with an initial decrease to a nadir point and a subsequent increase, described 20 traits, which included almost all those expressed in g/L of milk. The downward curve, continuously decreasing, described 14 traits, including almost all those in g/d; and the upward curve, continually increasing, described only the total κ-CN content in % N. This modeling approach helps in biological interpretation of the secretion pattern of different milk nitrogenous fractions. Moreover, it stimulates the discussion on the possible inclusion of novel traits in selection indices of dairy populations.

      ACKNOWLEDGMENTS

      The authors thank the Autonomous Province of Trento (Italy), who provided the funds for the project. The authors have not stated any conflicts of interest.

      REFERENCES

        • Amalfitano N.
        • Cipolat-Gotet C.
        • Cecchinato A.
        • Malacarne M.
        • Summer A.
        • Bittante G.
        Milk protein fractions strongly affect the patterns of coagulation, curd firming, and syneresis.
        J. Dairy Sci. 2019; 102 (30772026): 2903-2917
        • Amalfitano N.
        • Stocco G.
        • Maurmayr A.
        • Pegolo S.
        • Cecchinato A.
        • Bittante G.
        Quantitative and qualitative detailed milk protein profiles of 6 cattle breeds: Sources of variation and contribution of protein genetic variants.
        J. Dairy Sci. 2020; 103 (33069399): 11190-11208
        • Arnould V.M.R.
        • Hammami H.
        • Soyeurt H.
        • Gengler N.
        Short communication: Genetic variation of saturated fatty acids in Holsteins in the Walloon region of Belgium.
        J. Dairy Sci. 2010; 93 (20723713): 4391-4397
        • Auldist M.J.
        • Hubble I.B.
        Effects of mastitis on raw milk and dairy products.
        Aust. J. Dairy Technol. 1998; 53: 28-36
        • Auldist M.J.
        • Walsh B.J.
        • Thomson N.A.
        Seasonal and lactational influences on bovine milk composition in New Zealand.
        J. Dairy Res. 1998; 65 (9718493): 401-411
        • Baul S.
        • Cziszter L.
        • Acatincai S.
        • Gavojdian D.
        • Erina S.
        • Marcu A.
        • Buzamat G.
        • George G.
        Seasonal influences on milk yield and composition dynamics during a normal lactation in dairy cows: Milk yield, fat and protein percentage.
        Lucr. Stiint. Zooteh. Biotehnol. 2014; 47: 260-265
        • Bobbo T.
        • Ruegg P.L.
        • Stocco G.
        • Fiore E.
        • Gianesella M.
        • Morgante M.
        • Pasotto D.
        • Bittante G.
        • Cecchinato A.
        Associations between pathogen-specific cases of subclinical mastitis and milk yield, quality, protein composition, and cheese-making traits in dairy cows.
        J. Dairy Sci. 2017; 100 (28365113): 4868-4883
        • Bonfatti V.
        • Di Martino G.
        • Cecchinato A.
        • Vicario D.
        • Carnier P.
        Effects of β-κ-casein (CSN2–CSN3) haplotypes and β-lactoglobulin (BLG) genotypes on milk production traits and detailed protein composition of individual milk of Simmental cows.
        J. Dairy Sci. 2010; 93 (20655450): 3797-3808
        • Bouallegue M.
        • M'Hamdi N.
        Mathematical modeling of lactation curves: A review of parametric model.
        in: M'Hamdi N. Lactation in Farm Animals: Biology, Physiological Basis, Nutritional Requirements, and Modelization. IntechOpen Limited, 2020
        • Caffin J.P.
        • Poutrel B.
        • Rainard P.
        Physiological and pathological factors influencing bovine α-Lactalbumin and β-Lactoglobulin concentrations in milk.
        J. Dairy Sci. 1985; 68 (3842846): 1087-1094
        • Cipolat-Gotet C.
        • Cecchinato A.
        • Malacarne M.
        • Bittante G.
        • Summer A.
        Variations in milk protein fractions affect the efficiency of the cheese-making process.
        J. Dairy Sci. 2018; 101 (30122415): 8788-8804
        • Crittenden R.G.
        • Bennett L.E.
        Cow's milk allergy: A complex disorder.
        J. Am. Coll. Nutr. 2005; 24 (16373958): 582S-591S
        • Fox P.F.
        • Uniacke-Lowe T.
        • McSweeney P.L.H.
        • O'Mahony J.A.
        Lactose.
        in: Dairy Chemistry and Biochemistry. Springer International Publishing, 2015: 21-68
        • France T.C.
        • O'Mahony J.A.
        • Kelly A.L.
        The plasmin system in milk and dairy products.
        in: Kelly A.L. Larsen L.B. Agents of Change. Springer International Publishing, 2021: 11-55
        • García S.C.
        • Holmes C.W.
        Lactation curves of autumn- and spring-calved cows in pasture-based dairy systems.
        Livest. Prod. Sci. 2001; 68: 189-203
        • Gengler N.
        Persistency of lactation yields: A review.
        Interbull Bull. 1996; 12: 87-96
        • Gołebiewski M.
        • Brzozowski P.
        • Gołebiewski Ł.
        Analysis of lactation curves, milk constituents, somatic cell count and urea in milk of cows by the mathematical model of wood.
        Acta Vet. Brno. 2011; 80: 73-80
        • Gustavsson F.
        • Buitenhuis A.J.
        • Johansson M.
        • Bertelsen H.P.
        • Glantz M.
        • Poulsen N.A.
        • Lindmark Månsson H.
        • Stålhammar H.
        • Larsen L.B.
        • Bendixen C.
        • Paulsson M.
        • Andrén A.
        Effects of breed and casein genetic variants on protein profile in milk from Swedish Red, Danish Holstein, and Danish Jersey cows.
        J. Dairy Sci. 2014; 97 (24704225): 3866-3877
        • Hallén E.
        • Wedholm A.
        • Andrén A.
        • Lundén A.
        Effect of β-casein, κ-casein and β-lactoglobulin genotypes on concentration of milk protein variants.
        J. Anim. Breed. Genet. 2008; 125 (18363977): 119-129
        • Heck J.M.L.
        • Schennink A.
        • Van Valenberg H.J.F.
        • Bovenhuis H.
        • Visker M.H.P.W.
        • Van Arendonk J.A.M.
        • Van Hooijdonk A.C.M.
        Effects of milk protein variants on the protein composition of bovine milk.
        J. Dairy Sci. 2009; 92 (19233813): 1192-1202
        • ICAR (International Committee for Animal Recording)
        International agreement of recording practices—Guidelines approved by the General Assembly held in Cork, Ireland on June 2012.
        ICAR, 2012
        • Jensen H.B.
        • Holland J.W.
        • Poulsen N.A.
        • Larsen L.B.
        Milk protein genetic variants and isoforms identified in bovine milk representing extremes in coagulation properties.
        J. Dairy Sci. 2012; 95 (22612926): 2891-2903
        • Jõudu I.
        • Henno M.
        • Kaart T.
        • Püssa T.
        • Kärt O.
        The effect of milk protein contents on the rennet coagulation properties of milk from individual dairy cows.
        Int. Dairy J. 2008; 18: 964-967
        • Kroeker E.M.
        • Ng-Kwai-Hang K.F.
        • Hayes J.F.
        • Moxley J.E.
        Effects of environmental factors and milk protein polymorphism on composition of casein fraction in bovine milk.
        J. Dairy Sci. 1985; 68: 1752-1757
        • Macciotta N.P.P.
        • Dimauro C.
        • Rassu S.P.G.
        • Steri R.
        • Pulina G.
        The mathematical description of lactation curves in dairy cattle.
        Ital. J. Anim. Sci. 2011; 10: e51
        • Macciotta N.P.P.
        • Vicario D.
        • Cappio-Borlino A.
        Detection of different shapes of lactation curve for milk yield in dairy cattle by empirical mathematical models.
        J. Dairy Sci. 2005; 88 (15738251): 1178-1191
        • Maurmayr A.
        • Cecchinato A.
        • Grigoletto L.
        • Bittante G.
        Detection and quantification of αS1-, αS2-, β-, κ-casein, α-lactalbumin, β-lactoglobulin and lactoferrin in bovine milk by reverse-phase high- performance liquid chromatography.
        ACS Agric. Conspec. Sci. 2013; 78: 201-205
        • Maurmayr A.
        • Pegolo S.
        • Malchiodi F.
        • Bittante G.
        • Cecchinato A.
        Milk protein composition in purebred Holsteins and in first/second-generation crossbred cows from Swedish Red, Montbeliarde and Brown Swiss bulls.
        Animal. 2018; 12 (29307328): 2214-2220
        • Miciński J.
        • Jastrzebski M.
        • Klupczynski J.
        Yield and composition of milk from Polish Holstein-Friesian and Jersey cows in particular months of the first lactations as dependent on milk protein polymorphism (short communication).
        Arch. Tierzucht. 2008; 51: 216-223
        • Mota L.F.M.
        • Pegolo S.
        • Bisutti V.
        • Bittante G.
        • Cecchinato A.
        Genomic analysis of milk protein fractions in Brown Swiss cattle.
        Animals (Basel). 2020;
        • Ng-Kwai-Hang K.F.
        • Hayes J.F.
        • Moxley J.E.
        • Monardes H.G.
        Environmental influences on protein content and composition of bovine milk.
        J. Dairy Sci. 1982; 65 (6890960): 1993-1998
        • Ng-Kwai-Hang K.F.
        • Hayes J.F.
        • Moxley J.E.
        • Monardes H.G.
        Variability of test-day milk production and composition and relation of somatic cell counts with yield and compositional changes of bovine milk.
        J. Dairy Sci. 1984; 67 (6715630): 361-366
        • Ng-Kwai-Hang K.F.
        • Hayes J.F.
        • Moxley J.E.
        • Monardes H.G.
        Percentages of protein and nonprotein nitrogen with varying fat and somatic cells in bovine milk.
        J. Dairy Sci. 1985; 68 (3842864): 1257-1262
        • Ng-Kwai-Hang K.F.
        • Hayes J.F.
        • Moxley J.E.
        • Monardes H.G.
        Variation in milk protein concentrations associated with genetic polymorphism and environmental factors.
        J. Dairy Sci. 1987; 70 (3584600): 563-570
        • Ostersen S.
        • Foldager J.
        • Hermansen J.E.
        Effects of stage of lactation, milk protein genotype and body condition at calving on protein composition and renneting properties of bovine milk.
        J. Dairy Res. 1997; 64 (9161914): 207-219
        • Quinn N.
        • Killen L.
        • Buckley F.
        Empirical algebraic modelling of live weight of Irish dairy cows over lactation.
        Livest. Sci. 2006; 103: 141-147
        • Rutten M.J.M.
        • Bovenhuis H.
        • Heck J.M.L.
        • van Arendonk J.A.M.
        Predicting bovine milk protein composition based on Fourier transform infrared spectra.
        J. Dairy Sci. 2011; 94 (22032392): 5683-5690
        • Rzewuska K.
        • Strabel T.
        Genetic parameters for milk urea concentration and milk traits in Polish Holstein-Friesian cows.
        J. Appl. Genet. 2013; 54 (23934506): 473-482
        • Schmitz S.
        • Pfaffl M.W.
        • Meyer H.H.D.
        • Bruckmaier R.M.
        Short-term changes of mRNA expression of various inflammatory factors and milk proteins in mammary tissue during LPS-induced mastitis.
        Domest. Anim. Endocrinol. 2004; 26 (14757184): 111-126
        • Schopen G.C.B.
        • Heck J.M.L.
        • Bovenhuis H.
        • Visker M.H.P.W.
        • Van Valenberg H.J.F.
        • Van Arendonk J.A.M.
        Genetic parameters for major milk proteins in Dutch Holstein-Friesians.
        J. Dairy Sci. 2009; 92 (19233812): 1182-1191
        • Sebastiani C.
        • Arcangeli C.
        • Ciullo M.
        • Torricelli M.
        • Cinti G.
        • Fisichella S.
        • Biagetti M.
        Frequencies evaluation of β-casein gene polymorphisms in dairy cows reared in central Italy.
        Animals (Basel). 2020; 10 (32033348): 252
        • Silvestre A.M.
        • Martins A.M.
        • Santos V.A.
        • Ginja M.M.
        • Colaço J.A.
        Lactation curves for milk, fat and protein in dairy cows: A full approach.
        Livest. Sci. 2009; 122: 308-313
        • Silvestre A.M.
        • Petim-Batista F.
        • Colaço J.
        The accuracy of seven mathematical functions in modeling dairy cattle lactation curves based on test-day records from varying sample schemes.
        J. Dairy Sci. 2006; 89 (16606753): 1813-1821
        • Souza V.C.
        • White R.R.
        Variation in urea kinetics associated with ruminant species, dietary characteristics, and ruminal fermentation: A meta-analysis.
        J. Dairy Sci. 2021; 104 (33455789): 2935-2955
        • Stanton T.L.
        • Jones L.R.
        • Everett R.W.
        • Kachman S.D.
        Estimating milk, fat, and protein lactation curves with a test day model.
        J. Dairy Sci. 1992; 75 (1500567): 1691-1700
        • Stocco G.
        • Cipolat-Gotet C.
        • Bobbo T.
        • Cecchinato A.
        • Bittante G.
        Breed of cow and herd productivity affect milk composition and modeling of coagulation, curd firming, and syneresis.
        J. Dairy Sci. 2017; 100 (27837976): 129-145
        • Stocco G.
        • Cipolat-Gotet C.
        • Gasparotto V.
        • Cecchinato A.
        • Bittante G.
        Breed of cow and herd productivity affect milk nutrient recovery in curd, and cheese yield, efficiency and daily production.
        Animal. 2018; 12 (28712377): 434-444
        • Stoop W.M.
        • Bovenhuis H.
        • van Arendonk J.A.M.
        Genetic parameters for milk urea nitrogen in relation to milk production traits.
        J. Dairy Sci. 2007; 90 (17369239): 1981-1986
        • Strucken E.M.
        • de Koning D.J.
        • Rahmatalla S.A.
        • Brockmann G.A.
        Lactation curve models for estimating gene effects over a timeline.
        J. Dairy Sci. 2011; 94 (21183055): 442-449
        • Swanson K.M.
        • Stelwagen K.
        • Dobson J.
        • Henderson H.V.
        • Davis S.R.
        • Farr V.C.
        • Singh K.
        Transcriptome profiling of Streptococcus uberis-induced mastitis reveals fundamental differences between immune gene expression in the mammary gland and in a primary cell culture model.
        J. Dairy Sci. 2009; 92 (19109270): 117-129
        • Wedholm A.
        • Larsen L.B.
        • Lindmark-Månsson H.
        • Karlsson A.H.
        • Andrén A.
        Effect of protein composition on the cheese-making properties of milk from individual dairy cows.
        J. Dairy Sci. 2006; 89 (16899662): 3296-3305
        • Weimann C.
        • Meisel H.
        • Erhardt G.
        Short communication: Bovine κ-casein variants result in different angiotensin I converting enzyme (ACE) inhibitory peptides.
        J. Dairy Sci. 2009; 92 (19389946): 1885-1888
        • Wilmink J.B.M.
        Adjustment of test-day milk, fat and protein yield for age, season and stage of lactation.
        Livest. Prod. Sci. 1987; 16: 335-348
        • Wood G.M.
        • Boettcher P.J.
        • Jamrozik J.
        • Jansen G.B.
        • Kelton D.F.
        Estimation of genetic parameters for concentrations of milk urea nitrogen.
        J. Dairy Sci. 2003; 86 (12906064): 2462-2469
        • Wood P.D.P.
        Algebraic models of the lactation curves for milk, fat and protein production, with estimates of seasonal variation.
        Anim. Prod. 1976; 22: 35-40
        • Wood P.D.P.
        Breed variations in the shape of the lactation curve of cattle and their implications for efficiency.
        Anim. Prod. 1980; 31: 133-141