The evolution of the nutrient composition of mammalian milks
Summary
- In mammals, nutrient allocation during lactation is a critical component of maternal care as milk intake promotes juvenile growth and survival, and hence maternal and offspring fitness.
- Milk composition varies widely across mammals and is hypothesized to have arisen via selection pressures associated with environment, diet and life history. These hypotheses have been proposed based on observations and/or cross-species comparisons that did not standardize for stage of lactation and did not consider evolutionary history of the species in analyses.
- We conducted the largest comparative analysis of milk composition to date accounting for phylogenetic relationships among species in order to understand the selective advantage of producing milk with specific nutritional profiles. We examined four milk constituents in association with species ecology while incorporating phylogeny in analyses.
- Phylogenetic signal was apparent for all milk constituents examined. After controlling for phylogeny, diet and relative lactation length explained the greatest amount of variation in milk composition. Several aspects of species' ecologies, including adaptation to arid environments, reproductive output and maternal body mass were not associated with milk composition after accounting for phylogeny.
- Our results suggest that milk composition is largely a function of evolutionary history, maternal nutrient intake and duration of milk production. Arriving at these conclusions was made possible by including the evolutionary relationships among species.
Introduction
Few behaviours are as complex and variable among species than parental care, which is an important adaptation that has allowed species to thrive in diverse habitats and environmental conditions (Rosenblatt & Snowdon 1996). Although species differ in many components of parental care including which parent gives the care, the quantity of time, energy and resources invested, the type of care given and length of parental care (Clutton-Brock 1991), one of the best studied parental care behaviours in mammals is maternal nutritive investment in young. Nutritional investment occurs both in utero and during lactation in mammals and supports structural, neurological and organ development and functionality of the foetus or neonate. Nutritional provisioning is essential for offspring growth and survival, yet relatively little is known about the evolution of this important aspect of parental care. Because maternal physiological investment in milk synthesis is typically far greater than investment during gestation (Gittleman & Thompson 1988), understanding the selective pressures shaping lactation strategies, including the nutritional composition of milk, will give valuable insight into factors contributing to the evolution of offspring provisioning.
Across species, substantial variation exists in the proximate constituents of milk (i.e. protein, fat, sugars) and milk energy density. For example, milk fat content varies from 0·2% in the black rhinoceros, Diceros bicornis, to 60% in some species of Phocid seals (Oftedal & Iverson 1995). Milk sugars are virtually nonexistent in some pinnipeds and >11% in some Diprotodont marsupials (Oftedal 2000) and protein varies from just slightly >1% in some primate species to almost 16% in the eastern cottontail, Sylvilagus floridanus (Oftedal & Iverson 1995). The adaptive significance of this interspecific variation is unclear (Oftedal et al. 1993), but the specific composition of milk produced by a species has been postulated to have been selected for on the basis of the immunological, thermoregulatory, osmoregulatory or nutritive needs of the young (Payne & Wheeler 1968; Oftedal & Jenness 1988; Peddemors, Demuelenaere & Devchand 1989; Blackburn 1993; Kunz et al. 1995; Oftedal & Iverson 1995; Tilden & Oftedal 1997). As such, variation among species in life histories, diet, and habitat types are hypothesized to be associated with interspecific differences in lactation strategies (Oftedal 1984; Oftedal & Iverson 1995; Kunz & Hood 2000; Hinde & Milligan 2011). Unfortunately, our understanding of factors associated with the diversity of milk composition observed across mammals is currently limited.
Comparative studies of milk composition have experienced four main issues of conceptual and statistical nature: (i) studies have been largely qualitative, with species compared based on similar characteristics. (ii) Studies have been quantitative but neglected to incorporate species relatedness in statistical analyses. Not incorporating phylogenetic hypotheses is problematic because related species tend to resemble each other and thus confound the interpretation of biological data (Blackburn 1993) and because it violates assumptions of data independence (Felsenstein 1985; Garland & Adolph 1994; Garland, Bennett & Rezende 2005). (iii) Although more recent interspecific comparisons of milk composition often account for phylogeny, these studies have been based on restricted phylogenetic groups, such as within the order primates (Hinde & Milligan 2011), within Pinnipeds (Schulz & Bowen 2005), or within Eutherian mammals (Langer 2008), limiting our understanding of broad scale patterns. (iv) Some of these comparative studies (Langer 2008; Hinde & Milligan 2011) include milk data that were collected at different times within the lactation period. Because there is temporal variation in milk composition, comparing milks of different species taken from different time periods confounds interspecific variation with intraspecific temporal variation (Oftedal 1984; Oftedal & Iverson 1995). Here, we use a phylogenetic approach to test hypotheses proposed to explain differences among mammals (including Eutherians, Prototherians and Metatherians) in milk composition at mid-lactation.
Hypotheses and predictions
We test the significance of several variables that have been explicitly or implicitly proposed to contribute to differences in milk composition across mammals in addition to a few of our own predictions. Milk components included in our analyses are fat, protein, dry matter and energy. Sugars were excluded as some groups of mammals produce milk with only little or undetectable concentrations of sugar.
Body mass
Within particular phylogenetic groups, there are significant negative relationships between body mass and milk protein, dry matter, fat or energy density (Merchant et al. 1989; Derrickson, Jerrard & Oftedal 1996; Hinde & Milligan 2011). Species with small body masses are expected to produce more concentrated and energetically dense milk because they have higher mass-specific metabolic demands and reduced digestive capacity due to smaller gastrointestinal tracts (Blaxter 1961). Thus, small-bodied species should be incapable of ingesting greater quantities of milk to meet metabolic and nutritional demands and should instead require more concentrated and energy dense milk (Blaxter 1961).
Arid-adapted
Water balance is a concern for species inhabiting xeric environments (Schmidt-Nielsen & Schmidt-Nielsen 1952); thus, the production of highly concentrated milk might be expected in order to reduce maternal water loss through milk. However, some arid-adapted species, such as the camel, Camelus bactrianus, and zebras, Equus burchelli and Equus zebra, produce relatively dilute milk, which may have evolved to facilitate evaporative cooling of offspring (Oftedal & Iverson 1995). Whether an arid-adapted species produce concentrated or dilute milk may depend on body mass. Larger-bodied species that require water for evaporative cooling are predicted to have more dilute milk, whereas smaller-bodied species are expected to produce more concentrated milk because they do not depend on water to dissipate heat (Schmidt-Nielsen & Schmidt-Nielsen 1952). Therefore, we predict an interaction between body size and arid-adaptation on milk composition with small-bodied xeric species producing more concentrated milk.
Maternal diet
Differences in diet among species may contribute to interspecific variation in milk composition (Jenness & Sloan 1970; Leigh 1994; Kunz et al. 1995; Derrickson, Jerrard & Oftedal 1996) by contributing to differences in the availability of raw materials to the mammary gland for milk synthesis. Carnivorous species are thought to produce milk higher in fat, protein and energy density because animal matter contains more fat and protein than plants and fruits, which would increase the availability of fatty acids and amino acids that can be utilized for milk synthesis by the mammary gland (McNab 1980; Morrison 1980; Kunz & Diaz 1995). Thus, we expect milk composition to vary in accordance with diet, with herbivorous species producing milk lower in concentrations of fat and protein, and lower energy density relative to omnivorous and carnivorous species.
Lactation length
Species with longer lactation durations are predicted to produce milk lower in gross energy and total solids. Producing relatively dilute milk over a long lactation period might serve to protect maternal body stores from being depleted prior to completion of neonatal growth and development (Hinde & Milligan 2011). On the other hand, transferring large amounts of nutrients and energy to offspring in a short period of time is thought to reduce maternal maintenance costs during lactation while allowing a greater proportion of energy to be transferred to the young (Fedak & Anderson 1982).
Developmental stage at birth
Martin (1984) proposed that mammals producing precocial young would have more dilute milk than species producing altricial young because precocial young begin consuming solid food at an earlier age and would not require as much nutrient transfer from the mother during lactation. However, given that precocial young achieve independence earlier and they require more energy for thermoregulation than altricial species (Hackländer, Arnold & Ruf 2002), we predict that precocial species produce more concentrated and energy dense milk.
Biome (aquatic vs. terrestrial)
Aquatic mammals require a minimum blubber thickness that needs to be developed to reduce heat loss and achieve thermal balance (Drescher 1979; Worthy 1985). Postnatal thickening of the blubber layer may be facilitated by rapid maternal transfer of fat through milk. Thus, aquatic animals are predicted to produce more concentrated milk that is higher in fat and energy density to allow for rapid deposition of an insulating subcutaneous fat layer for neonatal thermoregulation (Jenness & Sloan 1970; Oftedal, Boness & Bowen 1988; Oftedal 1993; Oftedal & Iverson 1995).
Reproductive output
Lactation is the most energetically expensive phase of mammalian reproduction (Gittleman & Thompson 1988) and higher reproductive output, such as the production of a larger litter mass relative to maternal mass, increases energetic demands and milk energy output during lactation (König, Riester & Markl 1988). Thus, females with higher reproductive output are expected to produce more concentrated milk with higher fat and protein concentrations and greater energy density because of the greater nutritive and energetic requirements of the litter (Oftedal 1993; Power, Oftedal & Tardif 2002).
Materials and methods
Data collation
Species inclusion
We followed the species inclusion criteria of Oftedal & Iverson (1995). First, milk must have been collected from at least three individuals. Second, samples must have been collected by manual palpation from the mammary gland, not from the neonate's stomach or by using vacuum systems. Third, mother and young must not have been separated for an extended period of time (>24 h) for species that do not normally have intersuckling intervals of this length. Fourth, adequate information had to be provided to determine that milk was collected at mid-lactation. Domestic dogs, cats and agricultural species were not included in our analyses because of possible artificial selection on milk composition. Humans were excluded due to substantial variation among cultures and populations in breastfeeding practices, habitat and diet, which can impact milk synthesis.
Delineation of lactation stages
Only data on milk collected around mid-lactation (i.e. peak lactation) were included following the criteria of Oftedal & Iverson (1995) for delineation of lactation stages. Mid-lactation includes the period where milk composition is stable relative to other time points during lactation, whereas the initial and final changes in milk composition are considered to be indicative of early and late lactation, respectively. For marsupials, we considered mid-lactation to begin at the plateau in milk composition as young begin emerging from the pouch and ending shortly thereafter. For phocids, where young are weaned before beginning to consume solid foods, we considered mid-lactation to include the entire period extending from when milk composition becomes stable until weaning.
Milk composition
We began by including all milk data in Oftedal & Iverson (1995; n = 100 species) and then conducted a search for articles published between 1995 and 2011 using Thomson Reuters Web of Knowledge (Thomson Reuters 2012). From these articles, 30 additional species met the criteria established above and were included in our study. All milk composition data are expressed on a wet mass basis. Unless indicated otherwise, milk energy density (kcal per g) was calculated from the equation E = 9·11F + 5·86P + 3·95L, where the units for fat (F), protein (P) and sugars (L) are grams per gram of whole milk as in Derrickson, Jerrard & Oftedal (1996). Milk composition data along with analytical techniques employed and associated references are presented in Table S1 (Supporting information).
Ecology
Data on species' ecologies were extracted from compendiums and published articles when compendiums were missing data. If ranges of values were provided, the midpoint was calculated. For multiple measures on a single species, averages were obtained and we made sure not to count data from the same source multiple times. Developmental stage at birth was based on neonatal independence in four trait categories: thermoregulatory, sensory, locomotory and nutritionally, following Derrickson (1992). Species were coded with a 1 if the trait was present within 2 days of birth and a 0 if the trait was absent: hair covering body (thermoregulatory), eyes open (sensory), ambulatory without assistance (locomotory) and consumes solid food (nutritionally). Codes were summed across the four traits to give a ranking of developmental stage at birth from 0 to 4, with species receiving a 0 being the most altricial species and those receiving a 4 being the most precocial species. Because of the small sample size of species in the number 4 ranking (n = 2), species ranked as 3 or 4 were combined. Assignment of diet type (carnivorous, omnivorous or herbivorous) was based on the predominant food type in the diet. Reports of food items occasionally eaten were not included. Species were considered to be adapted to arid conditions if the predominant habitat type is characterized by mean annual precipitation <500 mm (McGinnies, Goldman & Paylore 1968). Predominant habitat type was assigned using the habitat classification system of Olson et al. (2001), based on range and habitat descriptions from Nowak (1999) and the IUCN Red List (2012). Species occupying multiple habitat types equally were not considered to be specifically arid-adapted. Reproductive output was estimated by dividing total litter mass by maternal mass. Total litter mass was calculated by multiplying neonate mass by litter size. As in Ernest (2003), lactation length was defined as the period of time between birth and offspring independence from milk. Lactation length was analysed relative to the total period of reproduction so that relative lactation length = absolute lactation length/(absolute gestation length + absolute lactation length). Ecological data are presented in Table S2 (Supporting information).
Phylogeny
We constructed a composite phylogeny (Fig. 1) using the program Mesquite (Maddison & Maddison 2011) based on published phylogenetic trees. Our mammalian tree was taken from Bininda-Emonds et al. (2007) and was trimmed to remove species not included in our study. We added to the tree species for which we had milk composition data but that were not in the Bininda-Emonds tree (n = 5; Cervus canadensis nelsoni, Cervus elaphus hispanicus, Equus ferus przewalski, Papio anubis, Papio cynocephalus). Phylogenetic relationships among Cervus species and subspecies were obtained from Randi et al. (2001), Equus from Oakenfull, Lim & Ryder (2000) and Papio from Newman, Jolly & Rogers (2004). Branch lengths were taken from Bininda-Emonds et al. (2007) and for the five species added, branch lengths were assigned by choosing the arbitrary ultrametricize function in Mesquite. The final tree contained some soft polytomies (n = 6), whose branch lengths were set to equal zero in Mesquite.
Statistical analyses
All continuous variables were log10-transformed prior to analysis. Statistical analyses were performed using the matlab (Matlab 2011) regression v2.m program (Lavin et al. 2008) using different statistical models: ordinary least squares regression (i.e. conventional nonphylogenetic approach; OLS), phylogenetic generalized least squares regression (PGLS) and regression with an Ornstein–Uhlenbeck transformation (RegOU; Garland, Bennett & Rezende 2005). We also used SAS (2002) to confirm output obtained from OLS regressions in matlab regression v2.m.
The OLS regression assumes no correlation of residuals among species. The PGLS regression utilizes a process like Brownian motion (BM) character evolution to examine phylogenetic autocorrelation of residuals assuming BM evolution of some variables and simultaneously a lack thereof in other variables (Lavin et al. 2008). For PGLS, we included Pagel's lambda (λ), a test statistic based on a BM process that progressively removes phylogenetic structure by utilizing the variance-covariance matrix of branch lengths and proportionally reducing the amount of covariance among species (Pagel 1999; Lavin et al. 2008). A λ value of 1 corresponds to the original phylogeny used, whereas λ of 0 indicates lack of phylogenetic correlation with the trait. In regression v2.m, a bootstrapping method utilizing 2000 simulations was employed to estimate 95% confidence intervals for λ. Confidence intervals excluding 0 indicate that λ is significantly different from 0, and thus, the trait can be said to exhibit phylogenetic signal. RegOU is based on an Ornstein–Uhlenbeck (OU) process of evolution (Felsenstein 1988; Garland et al. 1993; Blomberg, Garland & Ives 2003). The OU transformation alters branch lengths to make the tree more or less hierarchical than the original (Blomberg, Garland & Ives 2003). Regression v2.m estimates the optimal OU transformation parameter (d) and a value of 1 indicates that the original phylogenetic tree best fits the data, whereas a d of 0 indicates that a star phylogeny better fits the data. A d parameter between 0 and 1 indicates that a tree with branch lengths between the original phylogeny and a star phylogeny best fits the data (Blomberg, Garland & Ives 2003).
For each milk constituent (e.g. fat, protein, dry matter and energy density), we developed progressively more complex models with the simplest models containing a single independent variable. We then ran models including all possible combinations of the independent variables that were significant from the simplest models that best fit the data. Species with any missing values for the life-history/ecology traits were removed from all models for that specific milk constituent so that different statistical models (e.g. OLS, PGLS and RegOU) within each constituent could be compared. However, models with different dependent variables could not be statistically compared because they include different species and have different sample sizes. For example, for some species data on dry matter content of milk were missing, so these species were removed from all statistical models where dry matter was the dependent variable, but these species were included in all models where the dependent variables were fat, protein and energy density (if those data were available).
The fit of the different statistical models to the data was determined using Akaike Information Criterion for small sample size (AICc). AICc values were calculated in the regression v2.m program by the equation: (−2 × ln ML likelihood) + (2 × number of parameters [p] × sample size [n])/(n−p−1). Smaller AICc values indicate better fit of that model to the data. Models with AICc values <2 units larger than the best model are also typically considered to have strong support (Burnham & Anderson 2002). Maximum likelihood ratio tests (LRTs) were also employed to compare the fit of RegOU to either its PGLS or OLS counterpart for the same model and to compare PGLS and OLS counterparts. The difference in the ln maximum likelihoods between models multiplied by 2 follows a chi-square distribution, with degrees of freedom being equal to the difference in parameters between the two models being compared. PGLS and RegOU counterparts, however, have the same number of parameters and thus 0 degrees of freedom. For these comparisons, a difference in twice the ln maximum likelihoods >3·84, which is the 95th percentile of a χ2 distribution with 1 degree of freedom, indicates a significant difference in the fit of the models to the data (Felsenstein 2004).
Maximum likelihood ratio tests were also used to compare more complex models to simpler models within OLS tests, within PGLS tests and within RegOU tests when models were a nested subset of another (for example, the model containing only biome could be compared to biome + diet but not to relative lactation length + diet). When PGLS or RegOU models are found to better fit the data than OLS models, based on AICc values or LRTs, it is indicative of phylogenetic signal in the residuals (Grafen 1989; Freckleton, Harvey & Pagel 2002; Blomberg, Garland & Ives 2003). For all statistical tests, α = 0·05.
Results
There were no significant relationships between maternal body mass and milk composition (including all four components of milk analysed) or between reproductive output and milk composition for all statistical models considered (i.e. OLS, PGLS and RegOU; all P > 0·05). Overall, based on AICc values (smaller is better), phylogenetic models (RegOU or PGLS) provided a better fit to the data than their nonphylogenetic counterparts (OLS; Tables S3–S6, Supporting information). Furthermore, LRTs showed that both phylogenetic models (i.e. RegOU and PGLS) were statistically significantly better than their nonphylogenetic counterparts (P < 0·05 for all models). Therefore, a few ecological variables including developmental stage at birth, and whether or not a species is arid-adapted, were not included in more complex models because they were only significant in the OLS statistical models. Estimates of the optimal OU transformation parameter (d) were between 0·7 and 0·9 (Tables S3–S6, Supporting information) indicating that for all models, the tree that provided the best fit to the data contained a hierarchical structure and more closely approximated the original phylogeny with untransformed branch lengths. All estimates of λ were >0·93 and all confidence intervals excluded 0, indicating phylogenetic signal in all milk constituents analysed.
For models of fat as the dependent variable, RegOU and PGLS models had the lowest AICc values (Table S3, Supporting information). Based on LRTs, RegOU models for diet, biome + diet, relative lactation length + diet and relative lactation length + diet + biome were significantly better than their PGLS counterparts (all P < 0·01). For all other models, PGLS was significantly better than their RegOU counterparts (all P < 0·01). Overall, the best models (i.e. the ones with the lowest AICc values) were the RegOU model of relative lactation length + diet followed by the RegOU models of diet alone and relative lactation length + diet + biome. The RegOU model of diet was not significantly different from the RegOU model of relative lactation length + diet (χ2 = 2·1, d.f. = 1, P = 0·14) nor was the RegOU model of diet significantly different from the RegOU model of relative lactation length + diet + biome (χ2 = 0·2, d.f. = 2, P = 0·90). In addition, the RegOU model of relative lactation length + diet was not significantly different from the RegOU model of relative lactation length + diet + biome (χ2 = 1·9, d.f. = 1, P = 0·39), indicating that biome and lactation length do not contribute significantly to milk fat concentration when incorporating diet. This was confirmed with a partial F-test (relative lactation length: partial F1,125 = 2·65, P = 0·11; biome: partial F1,125 = 1·13, P = 0·29).
For protein, the best models based on the lowest AICc values were the RegOU and PGLS models of diet (Table S4, Supporting information). For diet, the difference in AICc values between PGLS and RegOU was <2 units and PGLS was not significantly different from its RegOU counterpart (based on LRTs, all P > 0·05). For all models except diet, PGLS models were significantly better than their RegOU counterparts (based on LRTs, all P < 0·01). More complex models for effects of ecological variables on milk protein concentration were not developed because only diet had a significant effect on protein concentration for the phylogenetic models.
For dry matter, RegOU and PGLS models had the lowest AICc values (Table S5, Supporting information). PGLS models for diet, biome + diet, relative lactation length + diet and relative lactation length + diet + biome were not significantly different from the respective RegOU models (all P > 0·05). For all other models, PGLS models were significantly better than their RegOU counterparts (all P < 0·0001). Overall, the models with the lowest AICc values were the PGLS and RegOU models of relative lactation length + diet and relative lactation length + diet + biome. LRTs indicated that the PGLS model of relative lactation length + diet + biome was not significantly better than the PGLS model of relative lactation length + diet (χ2 = 1·10, d.f. = 1, P = 0·29) nor was the RegOU model of relative lactation length + diet + biome significantly different from the RegOU model of relative lactation length + diet (χ2 = 0·08, d.f. = 1, P = 0·78). A partial F-test confirmed that biome was not significantly associated with milk dry matter content (PGLS: partial F1,119 = 1·65, P = 0·20; RegOU: partial F1,119 = 1·34, P = 0·25) when relative lactation length and diet were included in the models.
For energy density, PGLS models had the lowest AICc values (Table S6, Supporting information). PGLS models were not significantly different from RegOU models for diet, biome + diet, relative lactation length + diet and relative lactation length + diet + biome (all P > 0·05). For all other models, PGLS models were significantly better than their RegOU counterparts (all P < 0·0001). PGLS and RegOU models of relative lactation length + diet and relative lactation length + diet + biome had the lowest AICc values. The PGLS model of relative lactation length + diet + biome was not significantly different from the PGLS model of relative lactation length + diet (χ2 = 2·92, d.f. = 1, P = 0·09) and the RegOU model of relative lactation length + diet + biome was not significantly different from the RegOU model of relative lactation length + diet (χ2 = 2·64, d.f. = 1, P = 0·10). A partial F-test confirmed that biome was not significantly associated with milk energy density (PGLS: partial F1,111 = 0·76, P = 0·39; RegOU: partial F1,111 = 0·80, P = 0·37) when relative lactation length and diet were included in the models.
Taking the best model for each milk constituent, diet had a significant impact on milk fat, protein and dry matter concentrations, and energy density (Table 1). Specifically, carnivores had significantly higher milk energy density and dry matter, fat and protein concentrations than omnivores and herbivores (Table 1, Fig. 2). Relative lactation length was negatively correlated with both milk dry matter concentration and energy density (Table 1).
| Variable | Coefficient | SE | Partial F | d.f. | P |
|---|---|---|---|---|---|
| Fat | |||||
| y-intercept | 0·94 | 0·15 | 39·40 | ||
| Carnivore vs. Herbivore | 0·43 | 0·09 | 23·37 | 1, 127 | <0·0001 |
| Carnivore vs. Omnivore | 0·32 | 0·06 | 27·20 | 1, 127 | <0·0001 |
| Herbivore vs. Omnivore | −0·10 | 0·07 | 1·88 | 1, 127 | 0·17 |
| Diet | 16·33 | 2, 127 | <0·0001 | ||
| Protein | |||||
| y-intercept | 0·84 | 0·12 | 52·91 | 1, 127 | |
| Carnivore vs. Herbivore | 0·19 | 0·06 | 10·86 | 1, 127 | 0·001 |
| Carnivore vs. Omnivore | 0·14 | 0·04 | 12·56 | 1, 127 | <0·0001 |
| Herbivore vs. Omnivore | −0·05 | 0·05 | 1·04 | 1, 127 | 0·31 |
| Diet | 7·57 | 2, 127 | 0·0008 | ||
| Dry Matter | |||||
| y-intercept | 1·35 | 0·07 | 35·13 | ||
| Carnivore vs. Herbivore | 0·15 | 0·03 | 20·49 | 1, 120 | <0·0001 |
| Carnivore vs. Omnivore | 0·21 | 0·04 | 22·82 | 1, 120 | <0·0001 |
| Herbivore vs. Omnivore | −0·07 | 0·04 | 3·34 | 1, 120 | 0·07 |
| Diet | 13·86 | 2, 120 | <0·0001 | ||
| Relative lactation length | −0·20 | 0·07 | 7·98 | 1, 120 | 0·006 |
| Energy | |||||
| y-intercept | 0·12 | 0·08 | 2·03 | ||
| Carnivore vs. Herbivore | 0·25 | 0·05 | 28·85 | 1, 112 | <0·0001 |
| Carnivore vs. Omnivore | 0·31 | 0·06 | 28·37 | 1, 112 | <0·0001 |
| Herbivore vs. Omnivore | −0·07 | 0·05 | 2·08 | 1, 112 | 0·15 |
| Diet | 17·98 | 2, 112 | <0·0001 | ||
| Relative lactation length | −0·21 | 0·10 | 4·45 | 1, 112 | 0·04 |
- Full models for milk constituents (fat, protein, dry matter concentration and energy density) are all regression with an Ornstein–Uhlenbeck transformation models which were the best fit to the data based on AICc values and ln maximum likelihood ratio tests (see 4).
Discussion
Numerous hypotheses have been proposed to explain why milk composition varies among mammalian species, but until now, it has been unclear to what extent species relatedness and ecology contribute to the observed interspecific variation in milk composition. Herein, we compiled the most comprehensive data set to date on the gross composition of mammalian milks, including monotremes, marsupials and placental mammals. We show that for all of the models in our study, statistical models incorporating phylogenetic relatedness provided a better fit to the data than conventional nonphylogenetic models (Tables S3–S6, Supporting information).
We also compared models using different branch length transformations: λ (PGLS with Pagel's λ) and the Ornstein–Uhlenbeck transformation (RegOU). λ estimates were >0, indicating that some of the similarity in milk composition among closely related species is due to common ancestry, and thus, milk composition can be said to exhibit phylogenetic signal. In our study, estimates of d suggested that a tree with a hierarchical structure between the one used and a star phylogeny provided the best fit to the data. LRTs showed that different models (RegOU and PGLS) performed differently depending on the particular model (Tables S3–S6, Supporting information).
Developmental stage at birth and whether a species is adapted to arid conditions were both significantly related to milk composition in nonphylogenetic (OLS) models. In these OLS models, arid-adapted species were found to produce more dilute milk and species producing young with a higher degree of altriciality produced milk higher in protein concentration. However, these variables were no longer significant when examined in the context of phylogeny (Tables S3–S6, Supporting information), suggesting that phylogenetic relationships may explain more of the variation in milk composition than stage of neonatal development or environmental conditions. The different conclusions reached when employing nonphylogenetic and phylogenetic comparative methods underscore the importance of conducting statistical analyses that account for phylogenetic relationships among species.
Body mass was not correlated with any milk components in any of the statistical models (OLS, PGLS, or RegOU models; Tables S3–S6, Supporting information). In a comparative study of the milk dry matter content of 62 species of Eutherian mammals, Langer (2008) also found no significant relationship between female body mass and per cent milk dry matter both prior to and after accounting for phylogeny. Similarly, no differences in milk composition were found when comparing species within Pteropus (Hood et al. 2001), Eulemur (Tilden & Oftedal 1997) and Equus (Oftedal & Jenness 1988) genera, despite substantial interspecific variation in body mass. In contrast, other studies have found negative allometric relationships between body mass and milk dry matter content in rodents (Derrickson, Jerrard & Oftedal 1996) and macropod marsupials (Merchant et al. 1989). However, allometric relationships on a broader range of taxa should be interpreted with caution when species relatedness is not accounted for as phylogeny may be more influential in driving these trends (Harvey & Pagel 1991). After accounting for phylogeny, Hinde & Milligan (2011) still found a negative relationship between body mass and milk energy density in primates. Although body mass might be an important predictor of milk composition among species within certain mammalian groups, it does not appear to be related to differences in gross milk composition across all mammals.
Species with greater reproductive output, estimated as mean litter mass as a proportion of maternal mass, were predicted to produce more concentrated and energy-rich milk due to the greater nutritive and energetic demands of the litter (Oftedal 1993; Power, Oftedal & Tardif 2002). However, reproductive output was not correlated with any milk components in any statistical model tested. The total amount of nutrients transferred to young, however, is not only a function of the concentration of milk constituents but also the quantity of milk produced (Landete-Castillejos et al. 2003; Hinde, Power & Oftedal 2009b). It is possible that species with higher reproductive output transfer a greater total amount of nutrients and energy to offspring rather than producing milk with greater nutrient concentrations.
Whether a species occupies an aquatic or terrestrial habitat was significantly correlated with milk dry matter and fat concentration, and energy density (Tables S3, S5, S6, Supporting information). However, when including diet and relative lactation length in models with biome, biome became non-significant and LRTs indicated that models containing biome did not fit the data significantly better than models excluding biome. Thus, relative lactation length and diet were the most influential variables contributing to milk composition after accounting for phylogeny (Table 1 and Tables S3–S6, Supporting information). This contrasts Langer (2008), where no relationship was found between diet and milk dry matter. We found that carnivorous species produce milk higher in fat, protein and dry matter concentration, and energy density than both herbivorous and omnivorous species (Fig. 2), a pattern that is also observed among frugivorous and insectivorous bats (Kunz & Stern 1995; Messer & ParryJones 1997). Carnivorous species typically ingest higher quantities of fat and protein given that animal matter contains a higher proportion of these nutrients, which is thought to contribute to variation in milk composition among species consuming different diets (McNab 1980; Morrison 1980; Kunz & Diaz 1995).
Relative lactation length was negatively correlated with milk dry matter concentration and energy density (Table 1). This lactation strategy likely exists to reduce the cost of milk synthesis and to prevent irreversible damage to maternal somatic tissues. Species that have long lactation lengths would unlikely be able to sustain the demands of lactation if they produced highly concentrated and energy dense milk, either due to ceilings on food intake rates or due to substantial self-maintenance costs associated with mobilizing greater quantities of body stores (Hinde & Milligan 2011) and therefore produce dilute milk. On the other hand, species with truncated lactation periods are thought to reduce the costs of lactation by transferring a large quantity of highly concentrated and energy dense milk. In this manner, a greater proportion of maternal energy stores can be turned into gains in neonatal mass rather than to sustaining maternal somatic maintenance (Fedak & Anderson 1982). Among the phocids, a taxonomic group for which detailed information on maternal energy transfer and pup growth during lactation are available, hooded seal (Cystophora cristata) mothers lose only 33% of stored fat to increase pup mass by 100% over a 4-day lactation period (Bowen, Boness & Oftedal 1987). For phocids with longer lactation lengths, such as the grey seal (Halichoerus grypus), females lose 84% of their stored fat over an average 18-day lactation period and pups have a 236% increase in mass (Fedak & Anderson 1982) and in the northern elephant seal (Mirounga angustirostris), females transfer c. 58% of their fat reserves to increase pup mass by almost 300% during the 26-day lactation period (Costa et al. 1986).
When interpreting relationships between lactation length and milk composition, caution should be exercised because age at weaning can be difficult to identify and has been defined inconsistently between studies. Age of weaning has been defined as the first day young are observed to consume solid foods, the time when milk transfer declines precipitously, or the time when young cease suckling (Martin 1984). We chose to only incorporate lactation length estimates that were defined as the time period between birth and when the young no longer consume milk (Ernest 2003). The cessation of milk consumption has been identified by direct observation of suckling behaviours, by examination of the amount of milk in the stomach, and by direct milk expression (Izard 1987; Schulz & Bowen 2004).
Although we excluded studies that did not adhere to the guidelines for milk sample collection, determination of lactation stage and minimum sample size (see 3) in an attempt to standardize data for comparative purposes, it is important to recognize other potential sources of error among published studies from which data were extracted. First, data were collated from studies employing different analytical techniques to determine milk composition. Solvent extraction methods such as the Roese-Gottlieb method for fat, Kjeldahl or CHN techniques for protein and the phenol-sulphuric acid method for sugars are considered the most reliable techniques for quantifying concentration of milk constituents (Oftedal & Iverson 1995). Second, species were not raised under identical conditions prior to collecting and analysing milk, which means that evolutionary differences among species will be to some degree confounded by immediate environmental effects on phenotype (Garland & Adolph 1991, 1994; Garland, Bennett & Rezende 2005). However, it is typically unfeasible to maintain animals under identical conditions, particularly when studying wild populations, necessitating the assumption that differences among species reflect evolutionary and genetically based differences (Lavin et al. 2008). Third, milk samples for some species were collected entirely from captive individuals while other studies were conducted on wild populations, which may influence milk synthesis (Munks et al. 1991; Rose & Flowers 2005). Captive species typically have ad libitum access to food, can achieve larger body sizes and can rely more heavily on food intake over mobilizing body stores to support milk synthesis, which may alter milk composition. However, some studies have found no differences in milk composition between captive and wild populations of the same species (Messer & ParryJones 1997; Power et al. 2008). Finally, substantial variation in milk composition can occur across individuals within species due to factors such as reproductive experience, maternal condition, maternal age and/or offspring sex (Georges et al. 2001; Landete-Castillejos et al. 2005; Hinde 2009a; Hinde, Power & Oftedal 2009b; Skibiel & Hood 2013). Thus, if milk was collected from a nonrandom sample of individuals, interspecific variation in milk composition may be, to some degree, confounded by intraspecific variation.
In conclusion, this study exemplifies the necessity to incorporate phylogenetic relationships among species in comparative studies. Not only did statistical models incorporating phylogenetic relationships among species provide a better fit to the data than conventional nonphylogenetic models, but also several variables that were correlated with milk composition when a star phylogeny was assumed were not significant when accounting for species relatedness. Most importantly, our results indicate that the selective pressures acting on milk composition include the shared evolutionary history among species, the diet consumed and the relative length of the lactation period. This study provides valuable insight into the factors favouring the evolution of one of the key components of mammalian parental care, nutritional provisioning during lactation.
Acknowledgements
We would like to thank Ted Garland, Jr. for access to and assistance with the regressionv2.m programme, Kevin Middleton for his help with Mesquite and Chris Hamilton for assistance with graphical presentation of the phylogenetic tree. Colleagues in the Hill laboratory at Auburn University, F. Stephen Dobson, Mary Mendonca and Mary Beth Voltura provided invaluable feedback on the manuscript. Comments received from three anonymous reviewers greatly improved the content of this manuscript.
