Advertisement

Annual rhythms of milk and milk fat and protein production in dairy cattle in the United States

Open ArchivePublished:December 04, 2018DOI:https://doi.org/10.3168/jds.2018-15040

      ABSTRACT

      An annual pattern of milk composition has been well recognized in dairy cattle, with the highest milk fat and protein concentration observed during the winter and lowest occurring in the summer; however, rhythms of milk yield and composition have not been well quantified. Cosinor rhythmometry is commonly used to model repeating daily and annual rhythms and allows determination of the amplitude (peak to mean), acrophase (time at peak), and period (time between peaks) of the rhythm. The objective of this study was to use cosinor rhythmometry to characterize the annual rhythms of milk yield and milk fat and protein concentration and yield using both national milk market and cow-level data. First, 10 yr of monthly average milk butterfat and protein concentration for each Federal Milk Marketing Order were obtained from the US Department of Agriculture Agricultural Marketing Service database. Fat and protein concentration fit a cosine function with a 12-mo period in all milk markets. We noted an interaction between milk marketing order and milk fat and protein concentration. The acrophase (time at peak) of the fat concentration rhythm ranged from December 4 to January 19 in all regions, whereas the rhythm of protein concentration peaked between December 27 and January 6. The amplitude (peak to mean) of the annual rhythm ranged from 0.07 to 0.14 percentage points for milk fat and from 0.08 to 0.12 percentage points for milk protein. The amplitude of the milk fat rhythm generally was lower in southern markets and higher in northern markets. Second, the annual rhythm of milk yield and milk fat and protein yield and concentration were analyzed in monthly test day data from 1,684 cows from 11 tiestall herds in Pennsylvania. Fat and protein concentration fit an annual rhythm in all herds, whereas milk and milk fat and protein yield only fit rhythms in 8 of the 11 herds. On average, milk yield peaked in April, fat and protein yield peaked in February, fat concentration peaked in January, and protein concentration peaked in December. Amplitudes of milk, fat, and protein yield averaged 0.82 kg, 55.3 g, and 30.4 g, respectively. Milk fat and protein concentration had average amplitudes of 0.12 and 0.07, respectively, similar to the results of the milk market data. Generally, milk yield and milk components fit annual rhythm regardless of parity or diacylglycerol O-acyltransferase 1 (DGAT1) K232A polymorphism, with only cows of the low-frequency AA genotype (5.2% of total cows) failing to fit rhythm of milk yield. In conclusion, the yearly rhythms of milk yield and fat and protein concentration and yield consistently occur regardless of region, herd, parity, or DGAT1 genotype and supports generation by a conserved endogenous annual rhythm.

      Key words

      INTRODUCTION

      Milk production and milk component concentrations follow a seasonal pattern across the year that is well recognized by dairy producers and nutritionists. Milk yield and milk fat and protein concentrations are typically greatest in the winter and reach a nadir in the summer. Seasonal variation has long been quantified and accounted for by animal breeders, but has not been well integrated into dairy management (
      • Wood P.D.P.
      A note on the repeatability of parameters of the lactation curve in cattle.
      ). The causes of these yearly changes are not fully understood and are often attributed to environmental factors such as heat stress or changes in forage quality. However, these yearly patterns may be the result of an endogenous annual rhythm controlling milk synthesis.
      In nature, annual rhythms occur as a mechanism for organisms to predict seasonal environmental changes before they occur to allow proactive adaptations. Many animal species exhibit yearly cycles of reproductive activity, hibernation, migration, hair growth, and feeding behavior, allowing them to better prepare for changes in weather conditions and food supply (
      • Lincoln G.A.
      • Clarke I.
      • Hut R.
      • Hazlerigg D.
      Characterizing a mammalian circannual pacemaker.
      ). A classic example of annual rhythmicity in production animals is seasonal breeding observed in ewes, which exhibit estrus only from early autumn to late January, restricting lambing to the spring (
      • Shelton M.
      • Morrow J.T.
      Effect of season on reproduction of Rambouillet ewes.
      ). These annual rhythms of estrus are the result of interactions between photoperiod and an endogenous physiological timekeeping mechanism (
      • Malpaux B.
      • Robinson J.E.
      • Wayne N.L.
      • Karsch F.J.
      Regulation of the onset of the breeding season of the ewe: importance of long days and of an endogenous reproductive rhythm.
      ). Feed intake of sheep may also be regulated via annual rhythms, as the secretion of leptin, ghrelin, and orexin is affected by photoperiod length (
      • Kirsz K.
      • Szczesna M.
      • Molik E.
      • Misztal T.
      • Wojtowicz A.K.
      • Zieba D.A.
      Seasonal changes in the interactions among leptin, ghrelin, and orexin in sheep.
      ). In dairy cattle, the secretion of prolactin (
      • Chew B.P.
      • Malven P.V.
      • Erb R.E.
      • Zamet C.N.
      • D'Amico M.F.
      • Colenbrander V.F.
      Variables associated with peripartum traits in dairy cows. IV. Seasonal relationships among temperature, photoperiod, and blood plasma prolactin.
      ) and serotonin (
      • Philo R.
      • Reiter R.J.
      A circannual rhythm in bovine pineal serotonin.
      ) follow melatonin-controlled annual rhythms, with prolactin levels peaking in summer and serotonin peaking in winter. Furthermore,
      • Piccione G.
      • Messina V.
      • Scianó S.
      • Assenza A.
      • Orefice T.
      • Vazzana I.
      • Zumbo A.
      Annual changes of some metabolical parameters in dairy cows in the Mediterranean area.
      determined that circulating concentrations of BHB, bilirubin, creatinine, and triglycerides followed annual rhythms in Italian Brown dairy cattle.
      Biological rhythms, including annual rhythms, can be analyzed using cosinor rhythmometry. This technique fits a linear form of the cosine function with a 12-mo period to multiple years of data to determine if it fits an annual rhythm. The degree to which the data follow an annual rhythm was assessed using a zero-amplitude test, which determined if the linear form of the cosine function modeled that data significantly better than the linear effect of time (
      ). Cosinor rhythmometry also determined the amplitude, difference from mean to peak, and acrophase, or time at peak, of a biological rhythm (
      • Bourdon L.
      • Buguet A.
      • Cucherat M.
      • Radomski M.W.
      Use of a spreadsheet program for circadian analysis of biological/physiological data.
      ). Quantifying the annual rhythms of milk and milk component production can improve management decisions and may provide insight into physiological mechanism governing annual rhythms. The objective of our study was to quantify the annual rhythms of milk and milk component production in dairy cattle in the United States and examine cow-specific factors affecting these rhythms.

      MATERIALS AND METHODS

      USDA Milk Market Data

      Milk composition data from January 2000 through October 2015 were obtained from the Agricultural Marketing Service agency of the USDA (Washington, DC). The data set contained producer bulk tank milk butterfat and true protein concentration determined primarily by mid-infrared reflectance spectroscopy, collected by dairy processers, and compiled within each marketing order by the USDA. Orders included: Northeast (order 1), Appalachian (5), Florida (6), Southeast (7), Upper Midwest (30), Central (32), Mideast (33), Pacific Northwest (124), Southwest (126), Arizona-Las Vegas (131), and Western (135). A map of the geographical locations of each order is provided in Figure 1. Butterfat data were available for all 11 orders; however, milk protein percentages were not available for the Appalachian, Florida, Southeast, or Arizona-Las Vegas regions. Only 4 yr of data were available from the Western region because the order was terminated in 2004 (
      • Tosi G.M.
      Milk in the western marketing area; Termination of the order.
      ).
      Figure thumbnail gr1
      Figure 1Map of United States Federal Milk Marketing Orders. The Western (135) marketing order was terminated in 2004. Image adapted from map provided courtesy of the USDA Agricultural Marketing Service (Washington, DC).
      All statistical analysis was performed using the MIXED procedure of SAS 9.4 (SAS Institute, Cary, NC). Monthly butterfat and protein concentration were fit to the linear form of a cosine function with a 12-mo period according to
      • Bourdon L.
      • Buguet A.
      • Cucherat M.
      • Radomski M.W.
      Use of a spreadsheet program for circadian analysis of biological/physiological data.
      using random regression, as described by
      • Niu M.
      • Ying Y.
      • Bartell P.A.
      • Harvatine K.J.
      The effects of feeding time on milk production, total-tract digestibility, and daily rhythms of feeding behavior and plasma metabolites and hormones in dairy cows.
      . The model tested the random effect of year and the fixed effects of marketing order, the linear form of a cosine function with a period of 12 mo, and the interaction of order and the linear form of cosine functions. Denominator degrees of freedom were estimated using the Kenward-Roger method and the AR(1) covariance structure was used. The fit of the 12-mo rhythms was determined using a zero amplitude test, which employs an F-test to compare the full model containing the linear form of the cosine function to a reduced linear model testing the linear effect of month (
      ). The acrophase and amplitude of each marketing order were calculated, with subsequent significance tests, using a spreadsheet program developed in Microsoft Excel (Microsoft Corp., Redmond, WA;
      • Bourdon L.
      • Buguet A.
      • Cucherat M.
      • Radomski M.W.
      Use of a spreadsheet program for circadian analysis of biological/physiological data.
      ).

      Cow-Level Data

      Data from a previously published experiment by
      • Vallimont J.E.
      • Dechow C.D.
      • Daubert J.M.
      • Dekleva M.W.
      • Blum J.W.
      • Barlieb C.M.
      • Liu W.
      • Varga G.A.
      • Heinrichs A.J.
      • Baumrucker C.R.
      Genetic parameters of feed intake, production, body weight, body condition score, and selected type traits of Holstein cows in commercial tie-stall barns.
      were used to examine the effect of seasonal rhythms on milk production at the individual cow level. Briefly, test-day milk yield and butterfat and protein concentration were collected from 1684 individual cows in 11 Pennsylvania tiestall dairy herds during the years 2002 to 2011. Of these cows, 731 had been tested for the diacylglycerol O-acyltransferase 1 (DGAT1) K232A polymorphism using a panel of 121 candidate gene markers (
      • Dekleva M.W.
      • Dechow C.D.
      • Daubert J.M.
      • Liu W.S.
      • Varga G.A.
      • Bauck S.
      • Woodward B.W.
      Short communication: Interactions of milk, fat, and protein yield genotypes with herd feeding characteristics.
      ; Igenity, Neogen Corporation, Lincoln, NE). Butterfat and protein yield were calculated as the product of milk yield and butterfat and protein concentration. Test days when cows were less than 40 DIM were removed from the analysis, as high milk fat was expected during early lactation.
      Milk, fat, and protein yield and milk fat and protein concentration were fit to the linear form of a cosine function with a 12-mo period using the MIXED procedure of SAS 9.4. Fixed effects included herd, DGAT1 genotype, lactation group (1, 2, and 3+ lactations), and DIM, along with the random effects of cow and test year. The levels of each main effect (herd, DGAT1 genotype, and lactation group), along with the average of all herds, were fit to a 12-mo cosine function and rhythm fit, amplitude, and acrophase were determined using the methodology described above for the milk market data.

      RESULTS

      USDA Milk Market Data

      A cosine function with a 12-mo rhythm was used to test for an annual rhythm in milk fat and protein concentration using USDA Federal Milk Market data. We noted an interaction of milk marketing order and the cosine function on fat concentration, so the fit, amplitude, and acrophase of the 12-mo cosine function were determined for each milk market order (Table 1). Fat concentration in all marketing orders fit the 12-mo rhythm (P < 0.0001). Mean fat percent was lowest in the Arizona-Las Vegas region (3.56%) and greatest in the Pacific Northwest (3.73%; Table 1). The amplitude of the rhythm of fat concentration ranged from 0.07 to 0.14 percentage units between orders. The Southwest, Central (Figure 2C), Southeast, Mideast, Appalachian, and Western regions had 12-mo rhythms with the greatest amplitudes (0.13 to 0.14 percentage units) and did not differ between each other (P > 0.10). The Northeast, Upper Midwest (Figure 2B), and Pacific Northwest marketing orders all exhibited an amplitude of 0.11 percentage units, which was lower than the Southwest, Central, Southeast, Mideast, Appalachian, and Western regions (P < 0.05). The amplitudes of fat concentration rhythms were lowest in Arizona-Las Vegas (0.09 percentage units) and Florida (0.07 percentage units), with Florida having a lower amplitude than Arizona-Las Vegas (P < 0.05; Figure 2A and D).
      Table 1The amplitude and acrophase of a cosine function with a 12-mo period fit to the regional monthly averages of milk fat and protein concentration
      ItemFMMO
      Federal Milk Marketing Order.
      RegionNMeanAmplitude,
      Distance from mean to peak of the 12-mo rhythm.
      %
      Acrophase,
      Date at peak of the 12-mo fitted rhythm as month and day of month.
      date
      P-value
      P-value for the zero-amplitude test corresponding to the fit of a 12-mo cosine rhythm.
      Fat, %1Northeast1893.71
      Values that do not share a superscript differ (P < 0.05).
      0.11
      Values that do not share a superscript differ (P < 0.05).
      Dec 31
      Values that do not share a superscript differ (P < 0.05).
      <0.001
      5Appalachian1903.66
      Values that do not share a superscript differ (P < 0.05).
      0.13
      Values that do not share a superscript differ (P < 0.05).
      Jan 17
      Values that do not share a superscript differ (P < 0.05).
      <0.001
      6Florida1893.61
      Values that do not share a superscript differ (P < 0.05).
      0.07
      Values that do not share a superscript differ (P < 0.05).
      Dec 4
      Values that do not share a superscript differ (P < 0.05).
      <0.001
      7Southeast1893.66
      Values that do not share a superscript differ (P < 0.05).
      0.14
      Values that do not share a superscript differ (P < 0.05).
      Jan 3
      Values that do not share a superscript differ (P < 0.05).
      <0.001
      30Upper Midwest1893.72
      Values that do not share a superscript differ (P < 0.05).
      0.11
      Values that do not share a superscript differ (P < 0.05).
      Dec 31
      Values that do not share a superscript differ (P < 0.05).
      <0.001
      32Central1893.67
      Values that do not share a superscript differ (P < 0.05).
      0.14
      Values that do not share a superscript differ (P < 0.05).
      Jan 19
      Values that do not share a superscript differ (P < 0.05).
      <0.001
      33Mideast1893.69
      Values that do not share a superscript differ (P < 0.05).
      0.13
      Values that do not share a superscript differ (P < 0.05).
      Dec 31
      Values that do not share a superscript differ (P < 0.05).
      <0.001
      124Pacific Northwest1893.73
      Values that do not share a superscript differ (P < 0.05).
      0.11
      Values that do not share a superscript differ (P < 0.05).
      Jan 12
      Values that do not share a superscript differ (P < 0.05).
      <0.001
      126Southwest1883.64
      Values that do not share a superscript differ (P < 0.05).
      0.14
      Values that do not share a superscript differ (P < 0.05).
      Dec 31
      Values that do not share a superscript differ (P < 0.05).
      <0.001
      131Arizona-Las Vegas1903.56
      Values that do not share a superscript differ (P < 0.05).
      0.09
      Values that do not share a superscript differ (P < 0.05).
      Dec 29
      Values that do not share a superscript differ (P < 0.05).
      <0.001
      135Western513.62
      Values that do not share a superscript differ (P < 0.05).
      0.13
      Values that do not share a superscript differ (P < 0.05).
      Jan 18
      Values that do not share a superscript differ (P < 0.05).
      <0.001
      True protein, %1Northeast1903.04
      Values that do not share a superscript differ (P < 0.05).
      0.08
      Values that do not share a superscript differ (P < 0.05).
      Dec 31<0.001
      30Upper Midwest1903.04
      Values that do not share a superscript differ (P < 0.05).
      0.09
      Values that do not share a superscript differ (P < 0.05).
      Dec 30<0.001
      32Central1893.07
      Values that do not share a superscript differ (P < 0.05).
      0.10
      Values that do not share a superscript differ (P < 0.05).
      Jan 6<0.001
      33Mideast1903.06
      Values that do not share a superscript differ (P < 0.05).
      0.09
      Values that do not share a superscript differ (P < 0.05).
      Dec 30<0.001
      124Pacific Northwest1893.10
      Values that do not share a superscript differ (P < 0.05).
      0.08
      Values that do not share a superscript differ (P < 0.05).
      Dec 27<0.001
      126Southwest1893.09
      Values that do not share a superscript differ (P < 0.05).
      0.10
      Values that do not share a superscript differ (P < 0.05).
      Dec 30<0.001
      135Western513.09
      Values that do not share a superscript differ (P < 0.05).
      0.08
      Values that do not share a superscript differ (P < 0.05).
      Jan 2<0.001
      a–f Values that do not share a superscript differ (P < 0.05).
      1 Federal Milk Marketing Order.
      2 Distance from mean to peak of the 12-mo rhythm.
      3 Date at peak of the 12-mo fitted rhythm as month and day of month.
      4 P-value for the zero-amplitude test corresponding to the fit of a 12-mo cosine rhythm.
      Figure thumbnail gr2
      Figure 2Annual rhythms of fat concentration in selected Federal Milk Market Orders. Orders include (A) Florida (order 6), (B) Upper Midwest (order 30), (C) Central (order 32), and (D) Arizona–Las Vegas (order 131). Data presented as LSM with error bars representing SEM of each month and the line representing a fitted cosine function with a 12-mo period. Characteristics of fitted rhythm are displayed, including mean, amplitude (Amp; peak minus mean), acrophase (Acro; time at peak), and P-value for the zero-amplitude test corresponding to the fit of the 12-mo cosine rhythm. Annual rhythms of fat concentration for the other milk market orders are displayed in (https://doi.org/10.3168/jds.2018-15040).
      The acrophase, or time at peak, of butterfat concentration ranged from December 4 to January 18 across milk marketing orders. In 6 of 11 marketing orders, including Arizona- Las Vegas, Northeast, Upper Midwest, Mideast, Southwest, and Southeast, the 12-mo rhythm of fat concentration peaked within 4 d of January 1 (P > 0.10). No difference was observed between the acrophases of the Western, Central, and Appalachian regions, and all peaked between January 17 and 19 (P > 0.10). Florida's annual rhythm peaked on December 4, markedly earlier than the other regions (P < 0.05).
      An interaction between milk marketing order and the cosine function was observed for protein concentration, so the fit, amplitude, and acrophase of the 12-mo rhythm were compared between marketing orders. Protein concentration fit the 12-mo function in all 7 regions with protein data available (P < 0.001; Table 1). Mean protein concentration was greatest in the Pacific Northwest, Southwest, and Western regions, followed by Central and Mideast (P < 0.05). The Northeast and Upper Midwest averaged 3.04% protein, the lowest of all regions (P < 0.05). The amplitudes of annual rhythms in all regions were between 0.08 and 0.10 percentage units. The greatest amplitudes (0.10 percentage units) were observed in the Central and Southwest (Figure 3B) regions. The Upper Midwest and Mideast had intermediate amplitudes (0.09 percentage units), whereas the Northeast (Figure 3A), Pacific Northwest, and Western regions displayed amplitudes of 0.08 percentage units (P < 0.05). The acrophases of annual protein concentration rhythms occurred between December 27 and January 6 in all regions (P > 0.10). These results are consistent with those observed by
      • Bailey K.W.
      • Jones C.M.
      • Heinrichs A.J.
      Economic returns to Holstein and Jersey herds under multiple component pricing.
      , who reported that 3-yr average milk production peaks were greatest in December and January in the Mideast milk market.
      Figure thumbnail gr3
      Figure 3Annual rhythms of protein concentration in selected Federal Milk Market Orders. Orders include (A) Northeast (order 1), (B) Mideast (order 33), and (C) Southwest (order 126). Data presented as LSM with error bars representing SEM of each month and the solid line representing a fitted cosine function with a 12-mo period. Characteristics of fitted rhythm are displayed including mean, amplitude (Amp; peak minus mean), acrophase (Acro; time at peak), and P-value for the zero-amplitude test corresponding to the fit of the 12-mo cosine rhythm. Annual rhythms of protein concentration for the other milk market orders are displayed in (https://doi.org/10.3168/jds.2018-15040).

      Cow-Level Data

      Data from 11 tiestall herds in Pennsylvania were analyzed to characterize seasonal patterns of yield of milk and milk fat and protein, determine how cow-level factors, including parity or DGAT1 genotype, influence annual production rhythms, and examine variation in annual rhythms between individual herds. Milk yield of the 11 herds fit an annual rhythm with an amplitude of 1.26 kg and an acrophase of March 13. A yearly rhythm of milk yield was observed in 10 of the 11 herds (P < 0.05; Supplemental Figure S3A and B, https://doi.org/10.3168/jds.2018-15040). Eight of the 10 herds that exhibited a yearly rhythm of milk yield peaked between March 13 and May 5, whereas the remaining 2 herds peaked in mid-summer. These results agree with previous data suggesting that the average amount of milk shipped from farms is greatest in March, April, and May (
      • Bailey K.W.
      • Jones C.M.
      • Heinrichs A.J.
      Economic returns to Holstein and Jersey herds under multiple component pricing.
      ). The amplitudes of the milk yield rhythms ranged from 0.84 to 1.64 kg/d, suggesting that annual rhythms may be responsible for as much as 3.3 kg/d difference in milk production across the year.
      Overall, the average fat yield of all herds fit a 12-mo rhythm with a mean of 1,313 g/d, an amplitude of 55.3 g/d, and an acrophase occurring on February 23 (P < 0.0001; Supplemental Figure S3C and D, https://doi.org/10.3168/jds.2018-15040). All but 1 of the 11 herds exhibited a yearly rhythm, whereas 1 exhibited a rhythm with a drastically lower amplitude than the other 9 herds (16.5 g/d; P < 0.0001). The 9 remaining herds had amplitudes ranging from 97.2 to 47.3 g/d. The nonrhythmic and low-amplitude herds peaked on July 5 and 6, respectively, in contrast to the other 9 herds, which peaked between January 22 and April 3. Average protein yield also fit a yearly rhythm, with a mean of 1,081 g, an amplitude of 30.4 g, and an acrophase occurring on February 19 (Supplemental Figure S3E and F, https://doi.org/10.3168/jds.2018-15040). A yearly rhythm was observed in 10 of the 11 herds. The amplitudes of annual rhythms of protein yield were highly variable between herds and ranged from 14.3 to 47.0 g g/d. However, we noted low variability in the acrophases, with 8 herds peaking between January 31 and March 1 and 1 herd peaking on April 24.
      Fat and protein concentration fit 12-mo rhythms in all herds, as expected, and exhibited much more consistent yearly patterns than milk, fat, and protein yield (P < 0.0001; Supplemental Figure S4A and B, https://doi.org/10.3168/jds.2018-15040). The average rhythm of fat concentration for the 11 herds peaked on January 25, with an amplitude of 0.12 percentage units. Of the 11 herds, 9 peaked within 30 d of January 25, with the remaining 2 peaking on March 27 and April 7. The amplitude of the yearly rhythm averaged 0.12 percentage units, equivalent to a 0.24-unit change in fat percent throughout the year, and was highly variable between herds (0.07 to 0.19 units). Milk protein concentration similarly fit an annual rhythm in all herds (P < 0.001), with an average amplitude of 0.07 units and the average acrophase occurring on December 30 (Supplemental Figure S4C and D, https://doi.org/10.3168/jds.2018-15040). All herds peaked between December 8 and January 23, with amplitudes ranging from 0.04 to 0.11 percentage units, suggesting that annual rhythms of protein percent are highly consistent between herds.
      The effect of parity on fit of the cosine function and the mean, amplitude, and acrophase of annual production rhythms was evaluated (Figure 4). Groups included cows in first, second, or third and greater lactations. Milk yield increased with increasing parity (P < 0.05) and fit an annual rhythm in all 3 parity groups (P < 0.0001). The amplitude of the 12-mo milk yield rhythm was greater for cows in their third and greater lactations than first and second lactation (P < 0.05), but was not different between first- and second-lactation cows (P > 0.10). Furthermore, acrophase did not differ between first and second lactation, but occurred over a month earlier in cows at third and greater lactation (P < 0.05).
      Figure thumbnail gr4
      Figure 4The effect of parity on annual rhythms of milk yield and composition. Annual rhythm by parity shown for (A) milk yield, (B) fat yield, (C) fat concentration, (D) protein yield, and (E) protein concentration. Data presented as LSM with error bars representing SEM of each month and the line representing a fitted cosine function with a period of 12 mo. Mean, amplitude (Amp; peak minus mean), acrophase (Acro; time at peak), and P-value for the zero-amplitude test corresponding to the fit of a 12-mo cosine rhythm are displayed in tables below each corresponding graph. Mean, amplitudes, and acrophases that do not share a superscript (a–c) differ (P < 0.05).
      A yearly rhythm of fat yield was also observed in all lactation groups (P < 0.0001). Mean fat yield increased with increasing parity, whereas cows in third and greater lactation had 80 and 55% greater amplitudes than first- and second-lactation cows, respectively (P < 0.05). The acrophase of fat yield did not differ between cows in second or third and greater lactations (P > 0.10), but occurred 23 and 35 d later in first-lactation cows compared with those in second and third lactation, respectively (P < 0.05). Protein yield fit a 12-mo rhythm in all lactation groups (P < 0.0001), and mean protein yield increased from 936 g in first lactation to 1,130 g in second lactation and 1,260 g in third and greater lactations (P < 0.05). Amplitude was greatest in third and greater lactation (73.1 g; P < 0.05), but did not differ between cows in first and second lactation (P > 0.10). The acrophase of the yearly rhythm was similar between second and third or greater lactations, occurring on March 11 for second-lactation cows and on March 1 for cows in their third and greater lactation (P > 0.10). Alternatively, the acrophase of the annual rhythm of protein yield occurred later in first-lactation cows (April 10; P < 0.05).
      Milk fat and protein concentration fit a 12-mo rhythm in all 3 parity groups (P < 0.0001). Mean milk fat concentration decreased with increasing parity (P < 0.05), but no difference in amplitude was observed between lactation groups (P > 0.10). Similarly, parity exerted no effect on the acrophase of fat concentration rhythms, with all peaking between January 14 and 25 (P > 0.10). Milk protein concentration fit a yearly rhythm in all 3 lactation groups (P < 0.05). Mean protein concentration did not differ between first and second lactation (P > 0.10), but was 8 and 6% lower in cows at the third or greater lactation, respectively (P < 0.05). No difference occurred in the acrophase or amplitude of the rhythms of protein concentration between lactation groups, with amplitudes ranging from 0.08 to 0.09 percentage units and all acrophases occurring between December 12 and 26 (P > 0.10).
      The DGAT1 genotype influenced average milk fat concentration, with the cows of the AA genotype averaging 4.19% milk fat, KA genotype averaging 3.94% fat, and KK genotype averaging 3.58% fat (Figure 5; P < 0.05). Fat yield did not differ between the AA (1,296 g) and KA genotypes (1,256 g; P > 0.10), but was lower the KK genotype (1,178 g; P < 0.05). These results are consistent with previous reports showing that the K232A polymorphism is positively associated with milk fat percent and yield (
      • Schennink A.
      • Stoop W.M.
      • Visker M.H.P.W.
      • Heck J.M.L.
      • Bovenhuis H.
      • Van Der Poel J.J.
      • Van Valenberg H.J.F.
      • Van Arendonk J.A.M.
      DGAT1 underlies large genetic variation in milk-fat composition of dairy cows.
      ). An annual rhythm of fat concentration was observed regardless of DGAT1 genotype (P < 0.0001), and DGAT1 genotype did not affect the amplitude or acrophase of the rhythm (P > 0.10). Similarly, fat yield fit a rhythm in all 3 genotypes (P < 0.0001), with no difference in mean fat yield. Amplitude was similar between all genotypes and ranged from 54.6 to 70.0 g (P > 0.10). Likewise, acrophase was unaffected by DGAT1 genotype and all rhythms peaked between February 15 and March 1. Average milk yield did not differ between the AA and KA genotypes (P > 0.10), but was 2.4 and 1.2 kg greater in cows with the KK genotype, respectively (P < 0.05). Milk yield fit a 12-mo rhythm in cows with the KA and KK genotype, but not AA. Between the KA and KK groups, no difference was observed for either amplitude or acrophase.
      Figure thumbnail gr5
      Figure 5The effect of DGAT1 K232A polymorphisms (AA, KA, KK) on annual rhythms of milk and milk fat concentration and yield. Effect of DGAT1 polymorphism on annual rhythms of (A) milk yield, (B) fat concentration, and (C) fat yield are shown. Data presented as LSM with error bars representing SEM of each month and the line representing a fitted cosine function with a 12-mo period. Mean, amplitude (Amp; peak minus mean), acrophase (Acro; time at peak), and P-value for the zero-amplitude test corresponding to the fit of a 12-mo cosine rhythm are displayed in tables below each corresponding graph. Means, amplitudes, and acrophases that do not share a superscript (a–c) differ (P < 0.05).

      DISCUSSION

      Regional milk market and cow-level data both confirmed the presence of 12-mo rhythms of fat and protein concentration, and cow-level data demonstrated a rhythm in milk and milk fat and protein yield, which are biologically and economically important. The consistency of these rhythms between years, regions, and herds suggests that the annual patterns of production evident on dairy farms may be the result of an endogenous biological rhythm rather than the outcome of environmental factors. Heat stress, the most notable environmental factor associated with seasonal production declines, typically occurs during July and August in the United States due to high temperature-humidity index (
      • Bohmanova J.
      • Misztal I.
      • Cole J.B.
      Temperature-humidity indices as indicators of milk production losses due to heat stress.
      ). If the summer decrease in production were strictly due to the acute effects of heat stress, yearly production would be expected to remain relatively stable, with a sharp decline in the summer followed by recovery in the fall when the stress is alleviated. However, the results of our study suggest a consistent and predictable cosine decrease from the peak of production to the nadir. For example, fat and protein concentration, which peak in January in the majority of regions, begin to decrease in February and March, much earlier than would be expected if the effect is simply a consequence of heat stress. Furthermore, milk yield was shown to reach a minimum between September and November rather than mid-summer, when heat stress is expected to be the greatest. Lastly, the rhythm is very consistent between years and the degree of heat stress is expected to differ between years, especially in more northern regions (Supplemental Figure S5, https://doi.org/10.3168/jds.2018-15040).
      In both the milk market and cow-level data sets, the annual rhythms of fat and protein concentration consistently peaked in mid-winter and reached a minimum in July. These results are consistent with previous reports showing that maximal fat and protein concentration occur during December and January (
      • Bailey K.W.
      • Jones C.M.
      • Heinrichs A.J.
      Economic returns to Holstein and Jersey herds under multiple component pricing.
      ). Alternatively, research studies with grazing cattle found that the occurrence of milk fat depression is greatest in May (
      • Carty C.I.
      • Fahey A.G.
      • Sheehy M.R.
      • Taylor S.
      • Lean I.J.
      • McAloon C.G.
      • O'Grady L.
      • Mulligan F.J.
      The prevalence, temporal and spatial trends in bulk tank equivalent milk fat depression in Irish milk recorded herds.
      ). However, changes in the fermentability of grass throughout the year likely affect the risk for milk fat depression to a greater degree than the annual rhythm.
      • Dunshea F.R.
      • Walker G.P.
      • Ostrowska E.
      • Doyle P.T.
      Seasonal variation in the concentrations of conjugated linoleic and trans fatty acids in milk fat from commercial dairy farms is associated with pasture and grazing management and supplementary feeding practices.
      observed that milk concentrations of trans-10 C18:1 isomers, associated with biohydrogenation-induced milk fat depression, were greatest in August and September, which is consistent with the trough of fat yield production observed in the current experiment.
      The yearly change in photoperiod with lengthening and shortening days is a strong cue for entraining annual rhythms, and greater amplitude rhythms are often associated with greater variation in the light-dark cycle across the year (
      • Hut R.A.
      • Paolucci S.
      • Dor R.
      • Kyriacou C.P.
      • Daan S.
      Latitudinal clines: An evolutionary view on biological rhythms.
      ). In the Federal Milk Market Order data, the lowest amplitude annual rhythm occurred in Florida, followed by Arizona-Las Vegas. Florida has the lowest latitude (Tallahassee, FL: 30.4° N) of all the regions investigated, and, consequently, had the smallest change in photoperiod across the year. Similarly, Arizona-Las Vegas has a lower latitude (Las Vegas, NV: 36.2° N) than the majority of the other marketing orders. However, the amplitude of the other 9 orders did not follow a pattern based on photoperiod length. For example, the Southeast and Southwest regions had the greatest amplitude rhythms, along with the Appalachian, Mideast and Western regions. The northernmost regions, including the Northeast, Upper Midwest, and Pacific Northwest, had intermediate amplitudes. One potential cause for the inconsistency between regions is the differences in housing systems among various US regions. In Florida and Arizona, a large proportion of cows are housed in shaded dry lots with access to natural lighting, whereas in northern states cows are more likely to be housed in enclosed tiestall or freestall barns under the influence of artificial lighting. However, this does not account for the high amplitudes observed in the Southeast and Southwest, where cows are also commonly housed in open dry lots.
      Data from dairy herds in the southern hemisphere would provide additional insight into the role of change in yearly photoperiod on annual rhythms of production. Because the southern hemisphere experiences reciprocal seasonality to the northern hemisphere, their annual rhythms would be inverted if photoperiod was the primary driver. Likewise, milk synthesis should be completely arrhythmic at the equator because there is no change in day length across the year. Data from New Zealand support this conjecture, as milk fat and protein concentration are greatest during their winter (June–August), whereas milk yield peaks in their spring and early summer (September–December;
      • Auldist M.J.
      • Walsh B.J.
      • Thomson N.A.
      Seasonal and lactational influences on bovine milk composition in New Zealand.
      ). Caution should be used when interpreting results of data from most countries in the southern hemisphere, however, because a larger percentage of herds within those countries use grazing systems and differences in nutrient composition of pasture across the year may confound the effects of the annual rhythm.
      Herd-level data suggested that, unlike fat and protein concentration, which peak between December and February, fat and protein yield reach an acrophase between February and April and milk yield reaches an acrophase between March and May in most herds. Between herds, variability was present in the annual rhythms of milk, fat, and protein yield. Two herds in particular stood out for having low amplitudes of fat and protein yield and acrophases of fat, protein, and milk yield that occurred about 4 to 6 mo later than the other herds. One herd also notably failed to demonstrate a yearly rhythm of milk yield. These herds possessed no identifiable feeding or management characteristics to differentiate them from the other 8 (
      • Vallimont J.E.
      • Dechow C.D.
      • Daubert J.M.
      • Dekleva M.W.
      • Blum J.W.
      • Barlieb C.M.
      • Liu W.
      • Varga G.A.
      • Heinrichs A.J.
      • Baumrucker C.R.
      Genetic parameters of feed intake, production, body weight, body condition score, and selected type traits of Holstein cows in commercial tie-stall barns.
      ). The idiosyncratic rhythms observed in these herds are perplexing, but may be related to undocumented changes in herd management that masked the endogenous annual rhythms of the cows. Potential factors that could override annual rhythms of production within an individual herd include differences in forage quality across the year or seasonal variation in reproductive performance, leading to uneven average DIM across the year. Herds employing grazing systems are especially susceptible to seasonal changes in forage quality, with predictable changes in nutrient composition of pasture due to differences in temperature and photoperiod throughout the year.
      We examined the genotype at the DGAT1 locus to determine its influence on the yearly rhythms of fat concentration, fat yield, and milk yield. The K232A polymorphism in the DGAT1 gene is responsible for up to 50% of the genetic variation in milk fat production (
      • Schennink A.
      • Stoop W.M.
      • Visker M.H.P.W.
      • Heck J.M.L.
      • Bovenhuis H.
      • Van Der Poel J.J.
      • Van Valenberg H.J.F.
      • Van Arendonk J.A.M.
      DGAT1 underlies large genetic variation in milk-fat composition of dairy cows.
      ). The AA genotype is associated with high milk fat concentration, the KK genotype is associated with low milk fat, and the KA heterozygote results in an intermediate phenotype (
      • Winter A.
      • Krämer W.
      • Werner F.A.O.
      • Kollers S.
      • Kata S.
      • Durstewitz G.
      • Buitkamp J.
      • Womack J.E.
      • Thaller G.
      • Fries R.
      Association of a lysine-232/alanine polymorphism in a bovine gene encoding acyl-CoA:diacylglycerol acyltransferase (DGAT1) with variation at a quantitative trait locus for milk fat content.
      ). The current experiment confirmed these previous results, with AA cows producing milk with 0.25% greater milk fat concentration than KA cows, which was 0.36% greater than KK. Annual rhythms of nearly all production variables persisted regardless of DGAT1 genotype, with a lone exception being milk yield in cows with the AA genotype. The reason for the lack of an annual rhythm in this genotype is unclear, but may be a consequence of its lower prevalence. Only 1,196 test-day observations from 40 unique cows carried the AA genotype, which was equivalent to 5% of the total cows in the study. This low number of cows, combined with management variability between herds, may have contributed to the inability to detect an annual rhythm. However, despite the low number of observations, annual rhythms were detected for fat concentration and yield within the AA group. The DGAT1 genotype of the animals did not affect the shape of the annual rhythms of fat percent or fat yield, with similar amplitudes and phases in all groups. It is important to note that the interaction of DGAT1 genotype and the annual rhythm was only explored in 1 region in the current experiment and should be extended to other regions of the United States and world in future work. However, current results are similar to those observed by
      • Duchemin S.
      • Bovenhuis H.
      • Stoop W.M.
      • Bouwman A.C.
      • van Arendonk J.A.M.
      • Visker M.H.P.W.
      Genetic correlation between composition of bovine milk fat in winter and summer, and DGAT1 and SCD1 by season interactions.
      , who observed that milk fat concentration was not affected by the interaction of DGAT1 genotype and season, although they did observe an interaction between DGAT1 and summer versus winter milk for UFA concentration.
      Lactation number had little effect on the shape of annual production rhythms. Increased amplitudes of milk, fat, and protein yield occurred in cows in their third or greater lactation. These results suggest that, as cows age, greater oscillation in the yearly rhythms of production occurs. Furthermore, greater lactation number was associated with slightly earlier peaks in yearly rhythms of milk protein and fat yield. The cause of changes is unclear, but may be related to overall increase in production observed as lactation number increases.
      Previous experiments have supported that dairy cattle are influenced by annual rhythms. Although cows are not generally considered to be seasonal breeders,
      • Hansen P.J.
      Seasonal modulation of puberty and the postpartum anestrus in cattle: A review.
      observed modest effects of season and photoperiod length on reproduction in cattle. Furthermore,
      • Kendall P.E.
      • Webster J.R.
      Season and physiological status affects the circadian body temperature rhythm of dairy cows.
      observed greater daily fluctuations in body temperature in the summer compared with winter. Several circulating hormones and metabolites also follow annual rhythms across the year.
      • Petitclerc D.
      • Peters R.R.
      • Chapin L.T.
      • Oxender W.D.
      • Refsal K.R.
      • Braun R.K.
      • Tucker H.A.
      Effect of blinding and pinealectomy on photoperiod and seasonal variations in secretion of prolactin in cattle.
      observed that circulating prolactin concentrations are over 6 times greater in summer than winter, and that this effect occurs even when melatonin signaling is blocked through blinding or pinealectomy. Plasma concentrations of bilirubin, creatinine, triglycerides, and BHB follow 12-mo rhythms, with BHB peaking on April 1, bilirubin peaking on July 14, creatinine peaking on June 12, and triglycerides peaking June 16 (
      • Piccione G.
      • Messina V.
      • Scianó S.
      • Assenza A.
      • Orefice T.
      • Vazzana I.
      • Zumbo A.
      Annual changes of some metabolical parameters in dairy cows in the Mediterranean area.
      ). The annual rhythms observed in these metabolites may be regulated by the same mechanism as the annual rhythms of production. The annual patterns of production may also be a consequence of seasonal changes in feed intake. Cattle typically increase feed intake in the winter, and
      • Ueda K.
      • Mitani T.
      • Kondo S.
      Relationship of rumen fill and fermentation to diurnal and seasonal variation of herbage intake in dairy cows grazed on perennial ryegrass pasture.
      reported DM and fiber in the rumen of dairy cows is greater in autumn than spring, with a reciprocal relationship for ruminal VFA concentration. Endocrine factors regulating intake have been reported to follow an annual rhythm in other animals. Circulating leptin concentrations vary across the year in sheep and Siberian hamsters (
      • Henry B.A.
      • Blache D.
      • Dunshea F.R.
      • Clarke I.J.
      Altered “set-point” of the hypothalamus determines effects of cortisol on food intake, adiposity, and metabolic substrates in sheep.
      ), and leptin, orexin, and ghrelin appear to be modulated by photoperiod in sheep (
      • Kirsz K.
      • Szczesna M.
      • Molik E.
      • Misztal T.
      • Wojtowicz A.K.
      • Zieba D.A.
      Seasonal changes in the interactions among leptin, ghrelin, and orexin in sheep.
      ).
      A seasonal rhythm in milk yield and composition may be adaptive to improve nutrition of the calf. From a biological perspective, a clear advantage exists to providing nursing offspring with high-energy milk during the winter months. In other organisms, seasonal rhythms function to maximize reproductive success by scheduling parturition to allow neonates to be born during periods of high food availability and favorable climatic conditions (
      • Lincoln G.A.
      • Anderson H.
      • Loudon A.
      Clock genes in calendar cells as the basis of annual timekeeping in mammals—A unifying hypothesis.
      ). As an important component of mammalian reproduction, it is reasonable to expect lactation to similarly follow an annual rhythm and provide more milk and higher fat and protein to neonates in the winter when energetic demands are greater. Seasonal reproduction in sheep and hamsters appears to be coordinated by interactions with photoperiod and circadian clock genes located within the hypothalamus. Melatonin is released from the pineal gland during the dark phase of the photoperiod and binds with high affinity to the pars tuberalis (PT) of the hypothalamus (
      • Johnston J.D.
      • Tournier B.B.
      • Andersson H.
      • Masson-Pévet M.
      • Lincoln G.A.
      • Hazlerigg D.G.
      Multiple effects of melatonin on rhythmic clock gene expression in the mammalian pars tuberalis.
      ). The duration of melatonin release affects the phase of the core clock genes Per and Cry within the PT, and which acts downstream to influence the pulsatile frequency of GnRH, thus affecting reproductive cyclicity (
      • Misztal T.
      • Romanowicz K.
      • Barcikowski B.
      Melatonin-a modulator of the GnRH/LH axis in sheep.
      ). Although this mechanism is well characterized for controlling seasonal reproduction, a dearth of research exists examining the presence of mechanisms governing annual rhythms of lactation.
      The amplitude, acrophase, and period length of annual rhythms can be manipulated through photoperiod alterations.
      • Gwinner E.
      Circannuale Rhythmen bei Tieren und ihre photoperiodische Synchronisation.
      demonstrated that yearly cycles of testes growth and molting in migrating common starlings (Sturnus vulgaris) can be shorted by accelerating changes in photoperiod. Furthermore, accelerated photoperiods can speed up the reproductive cycles of rainbow trout (Oncorhynchus mykiss;
      • Bon E.
      • Breton B.
      • Govoroun M.S.
      • Le Menn F.
      Effects of accelerated photoperiod regimes on the reproductive cycle of the female rainbow trout: II Seasonal variations of plasma gonadotropins (GTH I and GTH II) levels correlated with ovarian follicle growth and egg size.
      ). Photoperiod manipulation has been extensively studied for its effects on milk production (
      • Dahl G.E.
      • Buchanan B.A.
      • Tucker H.A.
      Photoperiodic effects on dairy cattle: A review.
      ). Cows consistently increase milk yield when placed in a photoperiod with over 12 h of light, with maximal response occurring in a 16 h of light, 8 h of dark photoperiod. Interestingly, the effects of photoperiod on milk fat and protein concentration are inconsistent, with the majority of experiments demonstrating no effect. Prolactin is a common signal between lactation response and photoperiod and may play a role in is a candidate hormone for potential endocrine control of annual milk synthesis rhythms. Plasma prolactin concentrations are elevated under long-day lighting (
      • Foster R.G.
      • Kreitzman L.J.
      Seasonal Reproduction in Mammals and Birds.
      ). The effect of prolactin on seasonal rhythms of lactation has not been examined, but may be a potential area for future research.
      There appears to be a paradox between the effects of long-day lighting and annual rhythms of milk yield. In natural conditions milk increases across the winter when the duration of the light cycle is less than 12 h, whereas photoperiod research suggests that milk yield should be greatest when more than 12 h of light is present. The effect may be explained by induction of photorefractoriness to the annual rhythm by constant long-day photoperiods. In other mammalian species, long-term exposure to a constant photoperiod causes spontaneous reversion, or photorefractoriness, of a seasonal physiological response to the state expected in the opposite photoperiod (
      • Lincoln G.A.
      • Johnston J.D.
      • Andersson H.
      • Wagner G.
      • Hazlerigg D.G.
      Photorefractoriness in mammals: Dissociating a seasonal timer from the circadian-based photoperiod response.
      ). Suffolk sheep exposed to a fixed, short, 8 of light, 16 h of dark photoperiod exhibited summer-type physiology, characterized by a decrease in serum LH concentrations and anestrous (
      • Robinson J.E.
      • Karsch F.J.
      Refractoriness to inductive day lengths terminates the breeding season of the Suffolk ewe.
      ). In long-day breeding Syrian hamsters (Mesocricetus auratus), gonadal degeneration spontaneously occurs after long-term administration of a fixed, 14 of light, 10 h of dark photoperiod. The mechanism for this phenomenon is thought to be through dissociation of the endogenous annual timer in the PT with melatonin-based signals from the light-dark cycle (
      • Lincoln G.A.
      • Johnston J.D.
      • Andersson H.
      • Wagner G.
      • Hazlerigg D.G.
      Photorefractoriness in mammals: Dissociating a seasonal timer from the circadian-based photoperiod response.
      ). The long-term treatment needed to elicit the effects of long-day photoperiod on milk synthesis provides additional support for this mechanism. In other species, a fixed photoperiod must be applied for 4 to 12 wk before induction of photorefractoriness, which is consistent with data suggesting that the effect of long-day lighting on milk yield does not manifest until 4 wk after administration (
      • Dahl G.E.
      • Buchanan B.A.
      • Tucker H.A.
      Photoperiodic effects on dairy cattle: A review.
      ). Although this mechanism seems to be a promising explanation for the incongruence between the annual rhythm and photoperiod research, it has not been well examined in cows and further research should be performed to investigate this effect.

      CONCLUSIONS

      Milk yield and milk fat and protein concentration and yield fit an annual rhythm at both the regional and cow levels. These rhythms appear to exist independent of environmental effects, such as heat stress, and may be the result of endogenous circannual rhythms. Rhythms of fat and protein concentration peak in the middle of winter and reach a nadir in the summer, whereas yields of milk, fat, and protein peak in the spring. Region of the United States appears to have mild effects on annual rhythms, with lower latitude regions having lower-amplitude rhythms. Moderate variability in annual rhythms exists between herds, but DGAT1 genotype and lactation group exert little effects. The annual rhythms of production should be considered by dairy producers, consultants, and researchers when benchmarking and making management decisions. The effects of diet changes or implementation of new technologies on herd performance should be interpreted within the context of the annual rhythm. For example, 3.6% milk fat may indicate suboptimal milk fat in January but normal milk fat in July. In addition, feeding a dietary supplement in July may appear to improve milk fat percent in the following months, but the increase may be merely a consequence of the annual rhythm of production.

      ACKNOWLEDGMENTS

      The authors gratefully acknowledge Jared Daubert (Penn State University, University Park) for assistance with data collection, the USDA Agricultural Marketing Service (Washington, DC) for the use of federal milk marketing order production data, the eleven Pennsylvania dairy herds that participated in the study, and Neogen (Lansing, MI) for genotyping assistance. Research supported in part by Agriculture and Food Research Initiative Competitive Grant no. 2015-67015-23358 and 2016-68008-25025 from the USDA National Institute of Food and Agriculture (Washington, DC) , NIH Training Grant no. GM108563 (Bethesda, MD), and Penn State University including USDA National Institute of Food and Agriculture Federal Appropriations under project number PEN04539 and accession number 1000803

      REFERENCES

        • 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
        • Bailey K.W.
        • Jones C.M.
        • Heinrichs A.J.
        Economic returns to Holstein and Jersey herds under multiple component pricing.
        J. Dairy Sci. 2005; 88 (15905457): 2269-2280
        • Bohmanova J.
        • Misztal I.
        • Cole J.B.
        Temperature-humidity indices as indicators of milk production losses due to heat stress.
        J. Dairy Sci. 2007; 90 (17369235): 1947-1956
        • Bon E.
        • Breton B.
        • Govoroun M.S.
        • Le Menn F.
        Effects of accelerated photoperiod regimes on the reproductive cycle of the female rainbow trout: II Seasonal variations of plasma gonadotropins (GTH I and GTH II) levels correlated with ovarian follicle growth and egg size.
        Fish Physiol. Biochem. 1999; 20: 143-154
        • Bourdon L.
        • Buguet A.
        • Cucherat M.
        • Radomski M.W.
        Use of a spreadsheet program for circadian analysis of biological/physiological data.
        Aviat. Space Environ. Med. 1995; 66 (7487815): 787-791
        • Carty C.I.
        • Fahey A.G.
        • Sheehy M.R.
        • Taylor S.
        • Lean I.J.
        • McAloon C.G.
        • O'Grady L.
        • Mulligan F.J.
        The prevalence, temporal and spatial trends in bulk tank equivalent milk fat depression in Irish milk recorded herds.
        Ir. Vet. J. 2017; 70: 14
        • Chew B.P.
        • Malven P.V.
        • Erb R.E.
        • Zamet C.N.
        • D'Amico M.F.
        • Colenbrander V.F.
        Variables associated with peripartum traits in dairy cows. IV. Seasonal relationships among temperature, photoperiod, and blood plasma prolactin.
        J. Dairy Sci. 1979; 62 (512138): 1394-1398
        • Dahl G.E.
        • Buchanan B.A.
        • Tucker H.A.
        Photoperiodic effects on dairy cattle: A review.
        J. Dairy Sci. 2000; 83 (10791806): 885-893
        • Dekleva M.W.
        • Dechow C.D.
        • Daubert J.M.
        • Liu W.S.
        • Varga G.A.
        • Bauck S.
        • Woodward B.W.
        Short communication: Interactions of milk, fat, and protein yield genotypes with herd feeding characteristics.
        J. Dairy Sci. 2012; 95 (22365236): 1559-1564
        • Duchemin S.
        • Bovenhuis H.
        • Stoop W.M.
        • Bouwman A.C.
        • van Arendonk J.A.M.
        • Visker M.H.P.W.
        Genetic correlation between composition of bovine milk fat in winter and summer, and DGAT1 and SCD1 by season interactions.
        J. Dairy Sci. 2013; 96 (23127906): 592-604
        • Dunshea F.R.
        • Walker G.P.
        • Ostrowska E.
        • Doyle P.T.
        Seasonal variation in the concentrations of conjugated linoleic and trans fatty acids in milk fat from commercial dairy farms is associated with pasture and grazing management and supplementary feeding practices.
        Aust. J. Exp. Agric. 2008; 48: 1062-1075
        • Foster R.G.
        • Kreitzman L.J.
        Seasonal Reproduction in Mammals and Birds.
        1st ed. Yale University Press, New Haven, CT2009
        • Gwinner E.
        Circannuale Rhythmen bei Tieren und ihre photoperiodische Synchronisation.
        Naturwissenschaften. 1981; 68 (7322206): 542-551
        • Hansen P.J.
        Seasonal modulation of puberty and the postpartum anestrus in cattle: A review.
        Livest. Prod. Sci. 1985; 12: 309-327
        • Henry B.A.
        • Blache D.
        • Dunshea F.R.
        • Clarke I.J.
        Altered “set-point” of the hypothalamus determines effects of cortisol on food intake, adiposity, and metabolic substrates in sheep.
        Domest. Anim. Endocrinol. 2010; 38 (19733031): 46-56
        • Hut R.A.
        • Paolucci S.
        • Dor R.
        • Kyriacou C.P.
        • Daan S.
        Latitudinal clines: An evolutionary view on biological rhythms.
        Proc. Biol. Sci. 2013; 280 (23825204): 20130433
        • Johnston J.D.
        • Tournier B.B.
        • Andersson H.
        • Masson-Pévet M.
        • Lincoln G.A.
        • Hazlerigg D.G.
        Multiple effects of melatonin on rhythmic clock gene expression in the mammalian pars tuberalis.
        Endocrinology. 2006; 147 (16269454): 959-965
        • Kendall P.E.
        • Webster J.R.
        Season and physiological status affects the circadian body temperature rhythm of dairy cows.
        Livest. Sci. 2009; 125: 155-160
        • Kirsz K.
        • Szczesna M.
        • Molik E.
        • Misztal T.
        • Wojtowicz A.K.
        • Zieba D.A.
        Seasonal changes in the interactions among leptin, ghrelin, and orexin in sheep.
        J. Anim. Sci. 2012; 90 (22785160): 2524-2531
      1. Koukkari W.L. Sothern R.B. Chronobiometry: Analyzing for Rhythms. Springer, Netherlands, Dordrecht, the Netherlands2006
        • Lincoln G.A.
        • Clarke I.
        • Hut R.
        • Hazlerigg D.
        Characterizing a mammalian circannual pacemaker.
        Science. 2006; 314 (17185605): 1941-1944
        • Lincoln G.A.
        • Anderson H.
        • Loudon A.
        Clock genes in calendar cells as the basis of annual timekeeping in mammals—A unifying hypothesis.
        J. Endocrinol. 2003; 179 (14529560): 1-13
        • Lincoln G.A.
        • Johnston J.D.
        • Andersson H.
        • Wagner G.
        • Hazlerigg D.G.
        Photorefractoriness in mammals: Dissociating a seasonal timer from the circadian-based photoperiod response.
        Endocrinology. 2005; 146 (15919753): 3782-3790
        • Malpaux B.
        • Robinson J.E.
        • Wayne N.L.
        • Karsch F.J.
        Regulation of the onset of the breeding season of the ewe: importance of long days and of an endogenous reproductive rhythm.
        J. Endocrinol. 1989; 122 (2769153): 269-278
        • Misztal T.
        • Romanowicz K.
        • Barcikowski B.
        Melatonin-a modulator of the GnRH/LH axis in sheep.
        Reprod. Biol. 2002; 2 (14666149): 267-275
        • Niu M.
        • Ying Y.
        • Bartell P.A.
        • Harvatine K.J.
        The effects of feeding time on milk production, total-tract digestibility, and daily rhythms of feeding behavior and plasma metabolites and hormones in dairy cows.
        J. Dairy Sci. 2014; 97 (25306274): 7764-7776
        • Petitclerc D.
        • Peters R.R.
        • Chapin L.T.
        • Oxender W.D.
        • Refsal K.R.
        • Braun R.K.
        • Tucker H.A.
        Effect of blinding and pinealectomy on photoperiod and seasonal variations in secretion of prolactin in cattle.
        Proc. Soc. Exp. Biol. Med. 1983; 174 (6634713): 205-211
        • Philo R.
        • Reiter R.J.
        A circannual rhythm in bovine pineal serotonin.
        Experientia. 1980; 36 (7418825): 664-665
        • Piccione G.
        • Messina V.
        • Scianó S.
        • Assenza A.
        • Orefice T.
        • Vazzana I.
        • Zumbo A.
        Annual changes of some metabolical parameters in dairy cows in the Mediterranean area.
        Vet. Arh. 2012; 82: 229-238
        • Robinson J.E.
        • Karsch F.J.
        Refractoriness to inductive day lengths terminates the breeding season of the Suffolk ewe.
        Biol. Reprod. 1984; 31: 656-663
        • Schennink A.
        • Stoop W.M.
        • Visker M.H.P.W.
        • Heck J.M.L.
        • Bovenhuis H.
        • Van Der Poel J.J.
        • Van Valenberg H.J.F.
        • Van Arendonk J.A.M.
        DGAT1 underlies large genetic variation in milk-fat composition of dairy cows.
        Anim. Genet. 2007; 38 (17894561): 467-473
        • Shelton M.
        • Morrow J.T.
        Effect of season on reproduction of Rambouillet ewes.
        J. Anim. Sci. 1965; 24: 795-799
        • Tosi G.M.
        Milk in the western marketing area; Termination of the order.
        Federal Register, Washington, DC2004
        • Ueda K.
        • Mitani T.
        • Kondo S.
        Relationship of rumen fill and fermentation to diurnal and seasonal variation of herbage intake in dairy cows grazed on perennial ryegrass pasture.
        Anim. Sci. J. 2016; 87 (26608355): 1148-1156
        • Vallimont J.E.
        • Dechow C.D.
        • Daubert J.M.
        • Dekleva M.W.
        • Blum J.W.
        • Barlieb C.M.
        • Liu W.
        • Varga G.A.
        • Heinrichs A.J.
        • Baumrucker C.R.
        Genetic parameters of feed intake, production, body weight, body condition score, and selected type traits of Holstein cows in commercial tie-stall barns.
        J. Dairy Sci. 2010; 93 (20855024): 4892-4901
        • Winter A.
        • Krämer W.
        • Werner F.A.O.
        • Kollers S.
        • Kata S.
        • Durstewitz G.
        • Buitkamp J.
        • Womack J.E.
        • Thaller G.
        • Fries R.
        Association of a lysine-232/alanine polymorphism in a bovine gene encoding acyl-CoA:diacylglycerol acyltransferase (DGAT1) with variation at a quantitative trait locus for milk fat content.
        Proc. Natl. Acad. Sci. USA. 2002; 99 (12077321): 9300-9305
        • Wood P.D.P.
        A note on the repeatability of parameters of the lactation curve in cattle.
        Anim. Prod. 1970; 12: 535-538