Environmental selection during the last ice age on the mother-to-infant transmission of vitamin D and fatty acids through breast milk

The FADS Gene Cluster

Allelic variation in the FADS gene complex corresponds strongly to geographic ancestry (14, 15), and is the most pronounced allele frequency difference between the Greenlandic Inuit and other human populations (16). Because modern Western cultures place much attention on the cardiovascular disease risks associated with diets rich in omega-6 versus omega-3 long-chain polyunsaturated fatty acids (LC-PUFAs) (17, 18), the high frequencies of these alleles in the Inuit have been primarily interpreted as an adaptation to their traditionally meat-rich (and omega-3–rich) diet (16, 19). However, and surprisingly, 95% of native Central and South Americans also show evidence of selection on the same FADS polymorphisms, yet these cultural groups do not traditionally consume diets nearly as rich in omega-3 fatty acids as do the Inuit, suggesting that their common ancestral population is characterized by one which did (13).

From an evolutionary perspective, the possibility that the selective force favoring these FADS alleles was a healthier ratio of omega-3 to omega-6 fatty acids in adults would indicate that the effect must have been very strong, as selective pressure is relatively low on genetic variation that influences later and post-reproductive years of life (20). In contrast, even very small positive effects on fertility and childhood survival have far greater selective power (20). Considering this, the selective benefits of these FADS polymorphisms are more likely connected to their significant role in modulating the relative levels of omega-3 and omega-6 fatty acids during gestation and postnatal growth of infants.

To elaborate, the fatty acids arachidonic acid, docosahexaenoic acid (DHA), and eicosapentaenoic acid (EPA) are essential to cognitive and visual development (21). Maternal polymorphisms in the FADS gene cluster have strong effects on the levels of EPA and DHA in breast milk [in humans (22⇓⇓⇓⇓–27) and in other mammals (27, 28)]. This suggests that lipid levels in a mother’s diet are optimized by specific polymorphisms in her FADS gene cluster during milk synthesis (Fig. 2A). These polymorphisms are therefore a potentially strong target of selection because fetuses and infants have a very limited ability to synthesize LC-PUFAs on their own and are dependent on their mother’s genotype to modulate their relative proportions (22⇓⇓⇓⇓–27). We propose that the phenotypic effects of FADS polymorphisms on the LC-PUFA content of breast milk was the primary target of selection in past human populations that consumed diets with compromised proportions of omega-3 and omega-6 fatty acids, such as is seen in the traditional diets of people living in the Arctic.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 2.

Overview of breast anatomy, development, milk production, and histology. (A) Genetic influences on milk production summarized from the main text. The rodent, human, and cow figures indicate the systems from which the results are known. (B) Three stages of mammary development alongside the genetic mutations known to increase ductal branching. Hormones induce the pubertal and gestation/lactation stages (denoted in pink). In vitamin D-deficient adipose tissue, ductal branching increases during these two stages of development. The embryonic stage (in green) is not hormonally induced; branching density is influenced by the ectodysplasin pathway (108) and, specifically, increased by EDAR V370A. (C) Cross-section of a ductal lumen showing the bilayer of luminal and myoepithelial cells. (D) Schematic of a desmosome, one of the main adhesive structures between mammary luminal cells. The orange lines represent desmosomal cadherins in the extracellular space. The gray lines in the purple area represent the intracellular filaments. (E) Close-up view of the desmocolin proteins that comprise the desmo-adhesome. An increase in the activity of the ectodyplasin pathway alters the relative proportions of Dsc2 and Dsc3, leading to a reduction in the adhesive strength.

EDAR.

The human EDAR V370A variant encodes a change in amino acid sequence in the highly conserved death domain of EDAR (29⇓⇓–32). EDAR functions as a protein receptor on a cell’s surface that activates the transcription factor NF-ĸB in the ectodysplasin pathway (29). Comparative studies show this pathway is functionally conserved across virtually all vertebrates, playing essential roles in the development of ectodermal structures from bird feathers to fish scales (33, 34).

EDAR is one of four genes implicated in hypohidrotic ectodermal dysplasia (HED), a set of ∼150 syndromes characterized by mild to severe defects in ectodermally derived structures, such as hair, teeth, breasts, and sweat glands (33, 35⇓–37). In contrast to HED, the EDAR V370A allele has the opposite effect. Genome-wide association studies (GWASs) of multiple Asian populations show that EDAR V370A is correlated with hair shaft caliber (38, 39), earlobe and chin shape (40), and a suite of morphological variants on teeth (40⇓⇓–43). Knock-in mouse studies reveal that the V370A allele leads to a twofold increase in NF-ĸB activation (30, 31, 44). Just like humans, the EDAR V370A mouse has thicker hair shafts, an increase in the branching density of mammary gland ducts, and (inconsistently) an increase in the number of eccrine glands on the footpads (30, 31).

Several studies conclude that EDAR V370A experienced a bout of intense positive selection ∼20,000 y ago in northern China (29, 31, 45). Previous interpretations as to why this selection occurred focused on the associated increase in eccrine sweat gland density on the fingertips, concluding that the selection was for improved thermoregulatory sweating during a warm spell during the LGM (31) or that the increase in sebaceous glands would offer protection from the cold, dry air of the LGM (30). These interpretations rely primarily on the current distribution of the allele in living populations (Fig. 3A), which shows high frequencies of EDAR V370A in North and East Asian and Native American populations but virtual absence in other populations around the world (29⇓–31, 45).

• Download figure
• Open in new tab
• Download powerpoint

Fig. 3.

Relationship between EDAR V370A and degree of incisor shoveling. (A) Current world-wide allele frequencies of EDAR V370A (in yellow) and other EDAR haplotypes (in blue) (data from ref. 29). Because these data are from living people, the Native American data include significant admixture from European colonization, a population with essentially no occurrence of EDAR V370A. These modern data likely vastly underrepresent the occurrence of EDAR V370A among indigenous people before European contact. Current Asian allele frequencies may be higher than pre-LGM frequencies due to back-migration from the Beringian population after the end of the LGM, a migration supported by linguistic data (46). (B) Histogram showing EDAR V370A genotype and degree of incisor shoveling, demonstrating the imperfect but strongly additive nature of EDAR V370A’s influence on incisor shoveling (adapted from ref. 43). (C) Frequencies of incisor shoveling scores observed in archaeological populations from Africa, Europe, South Asia, East Asia, North America, and South America (Table S1). Purple and blue represent a lack of incisor shoveling, and as such, an individual who is EDAR V370A^−/−. Note the lack of shoveling scores of 0 and 1, and very low occurrence of 2, in the indigenous people of the Western Hemisphere, indicating a very high frequency of EDAR V370A before European contact.

The post-LGM diaspora into the Americas and the subsequent centuries of European colonization (and genetic admixture) have dramatically overwritten the allelic variation of ∼20,000 y ago when EDAR V370A experienced a bout of selective pressure. To best interpret the evolutionary history of EDAR V370A, we need information on past allele frequencies. This task is difficult from a genome sequence approach because of the paucity of ancient DNA available. However, the pleiotropic effects of EDAR V370A provide a unique opportunity to more precisely reconstruct its evolutionary history without DNA sequence data.

Temporal and Geographic Variation of EDAR V370A

As noted previously, EDAR V370A is correlated with a suite of morphological variants on teeth, the most notable of which is incisor shoveling (40⇓⇓–43). GWASs report a close relationship between genotype (EDAR V370A) and this dental phenotype in living Chinese, Japanese, and Koreans (40⇓⇓–43). In these populations, incisor shoveling of any degree correlates with the presence of one or two copies of the EDAR V370A allele in an imperfect but clearly additive manner (41⇓–43) (Fig. 3B).

We scored the degree of expression in maxillary incisor shoveling for 5,333 people from >54 archaeological populations from across Europe, Asia, and North and South America (Supporting Information and Tables S1 and S2). Incisor shoveling is absent or expressed only to a slight degree outside of North and East Asian populations (where EDAR V370A frequencies are close to zero) but close to ubiquity in Native American populations, a pattern long interpreted as the result of an unknown selective pressure (47, 48) (Fig. 3C).

Considering the GWAS results that indicate an additive effect of EDAR V370A’s contribution to incisor shoveling, we interpreted a score of 1 or greater as evidence of an individual likely carrying at least one copy of EDAR V370A. There are no reported instances of a score of 0 among the 3,183 individuals assessed from North and South America, and as such, no evidence of any Native Americans before European contact without at least one copy of the EDAR V370A allele (Fig. 3C and Supporting Information). A series of ANOVAs comparing the frequency of each shoveling grade (range: 0–7) across geographic regions demonstrates that a score of 1 has the greatest variation between geographic groups, further suggesting that a score of 1 or higher is reflective of the presence of at least one copy of EDAR V370A (Supporting Information and Table S3).

These phenotypic data strongly suggest that incisor shoveling (and, concomitantly, the EDAR V370A allele) reached near-fixation in the population ancestral to all indigenous people of the Western Hemisphere. Given the timing of the dispersal into the Americas ∼17,000 y ago (3, 49), these dental data support the conclusion that EDAR V370A underwent positive selection in the ancestral Native American population during the LGM, similar to what we see with the FADS gene cluster (13), rather than in a population in what is now China (our evolutionary quantitative genetic analysis further supports this; Supporting Information).

The Ancestral Population of Native Americans and the LGM

Indigenous people in the Western Hemisphere derive from people who occupied latitudes above 55°N in Asia before 40,000 y ago (50). Well-dated archaeological sites near the mouth of the Yana River indicate a year-round adaptation to the Beringian Arctic Zone (latitude of 70°N) by 32,000 calibrated radiocarbon years (cal) B.P. (51).

Between 28,000 and 18,000 y ago, during the LGM, as aridity increased and biological productivity reduced dramatically (52), plants and animals in many parts of the world retreated to refugia, leaving signals of genetic isolation and hybridization that are documented widely by molecular ecologists (53, 54). Human populations similarly abandoned arid regions in Africa, Eurasia, and Australia, as indicated by genetic bottlenecks and evidence of local settlement hiatus (e.g., refs. 55, 56).

Various lines of evidence suggest that the Native American lineages during the LGM were genetically isolated in one or more refugia within Beringia (2, 5⇓–7), and at least one archaeological site appears to confirm a human presence in northeastern Beringia (latitude of 67°N) during that time frame (57). Analyses of Native American genetic variation indicate that the ancestral population was isolated for up to 9,000 y before dispersal in North and South America (after 15,000 cal B.P.), with an effective population size of at least a few thousand individuals (3, 49, 58).

The “Beringian Standstill” model places the ancestral Native American population in the Beringian refugium during much or all of the LGM (2, 5⇓–7), where it was geographically isolated from the Siberian population by uninhabitable areas in northeastern Asia and unable to expand into the Americas because of the coalesced Laurentide and Cordilleran ice sheets (Fig. 1B). The genetic isolation is a key point, as evolution occurs more readily in the absence of gene flow (59). The FADS polymorphisms and EDAR V370A likely existed in the pre-LGM ancestral populations [although the allele is not present in Denisovans (60)]. However, there was a significant opportunity for selection to shift the frequencies of functional alleles dramatically when this population became genetically isolated in the extreme UV environment of Beringia, including Arctic Beringia (above 66°N).

Biological Consequences of Low UV

Most vitamin D necessary for human health and reproductive success is produced through biosynthesis initiated by skin exposure to UV-B photons (61, 62). The lower rates of cutaneous synthesis of vitamin D that result from low UV environments have a wide range of health consequences. Vitamin D is among a large family of fat-soluble micronutrients that accumulate in adipose tissue proportional to circulating serum levels, and is released slowly as serum levels reduce (63). The bioactive form of vitamin D plays a secosteroid hormonal role in permitting absorption of dietary calcium through the lining of the gut (64). Through the Vitamin D receptor (Vdr), vitamin D also regulates expression of more than 220 genes in the human genome, with significant involvement in immune function and autoimmune disorders (65, 66). It is therefore not surprising that variation in immune function and disease risk correlates with latitude and UV-B exposure (67). Vitamin D also plays an important immunomodulatory role at the maternal/fetal interface (the placenta) (68) and is associated with sex-specific variation in birth weight (69).

However important the immunological role of vitamin D, adipose tissue appears to be the main phenotypic target of Vdr (63, 65, 66, 70). Vitamin D deficiency significantly compromises the metabolic function of adipose tissue, with a wide range of deleterious health effects (63). Body fat has a remarkable degree of phenotypic plasticity that enables a wide range of essential physiological functions, from maintaining body heat, to responding to cold, to forming the pink adipocytes in the mammary glands during pregnancy in anticipation of lactation (71⇓–73) (Fig. 2A).

Nutritional studies have long shown interrelationships between vitamin D and LC-PUFAs (74⇓–76). For example, the vitamin D synthesized through skin exposure to UV-B (vitamin D3) has positive effects on lipid profiles that cannot be replicated with dietary vitamin D2 supplementation (77). These interrelationships result from the effects of vitamin D on adipose tissue and immunological function, which are modulated by environmental exposure to UV-B.

Human populations that dispersed into geographic regions above 30° latitude must have experienced selection for reduced cutaneous eumelanin pigmentation to facilitate vitamin D biosynthesis (78, 79). Selection for depigmented skin phenotypes has occurred at least twice in modern human evolution and once in Neanderthals (78, 80), with selective pressure readily detected in genomic analyses (81⇓–83).

At latitudes above 48° (e.g., 450 km north of Hokkaido, Japan; Fig. 1A), sunlight reaching the Earth’s surface contains almost no UV-B except for low levels at or near the summer solstice (84, 85). Year-round human habitation of even more extreme high latitudes leads to vitamin D deficiency that cannot be counterbalanced by depigmentation alone (84) (Fig. 1A). In these geographic regions, diets centered on vitamin D-rich foods, such as marine mammals, oily fish, reindeer, or caribou, have been adopted across a range of traditional Arctic cultures (85). Archaeological evidence shows that these dietary innovations correlate with the earliest occupation of these latitudes (55, 86, 87) and shifts in allele frequencies of the genes involved in fatty acid synthesis (16).

Dietary practices can mitigate vitamin D deficiency for older children and adults, but pregnant and nursing mothers and their breast-feeding infants are still at significant risk because they need particularly high levels of vitamin D to avoid the litany of deleterious effects of vitamin D insufficiency (88, 89).

Fat and Vitamin D Content in Milk

Human breast milk varies greatly in fat content, with the relative proportions of omega-3 and omega-6 fatty acids varying by ethnicity independent of maternal diet (90), most likely because of polymorphisms in the FADS gene cluster as discussed earlier (22⇓⇓⇓⇓–27). The LC-PUFAs in milk are drawn from the specialized pink adipocytes found only in mammary glands, which are characterized by a higher number of intracellular lipid droplets (71⇓–73) and are dependent on healthy vitamin D levels to function properly (63).

Vitamin D content in human milk also varies between individuals and across populations (90⇓⇓–93). Although the effects of dietary supplementation, sun exposure, and genetic variation affect the vitamin’s concentration within an individual, there is nevertheless a strong correlation between maternal and infant circulating vitamin D levels (94, 95), and both mother and infant levels are related to each other’s sun exposure (93, 96). Importantly, vitamin D levels in milk are about threefold lower than the circulating levels in the mother (97). Endogenously produced vitamin D3 from UV-B exposure appears to be preferentially transferred into milk compared with the vitamin D2 consumed through diet (93). It is thus reasonable to suggest that infants in extremely low UV environments would be particularly dependent on the mother’s potential capacity to transfer vitamin D from her adipose stores across the placenta and also into her breast milk.

Fat, Vitamin D, and the Mammary Gland

Given that vitamin D is intricately intertwined with adipose physiology, it is not surprising that vitamin D plays a role in virtually all aspects of mammary gland function, as evinced by the presence and activation of the Vdr in mesenchymal and epithelial breast tissue (98, 99). Vitamin D also appears to influence ductal branching and, with it, milk content.

There are three distinct phases of breast development (Fig. 2B): (i) hormone-independent embryonic development during which buds form, sprout, and start to branch into the precursor of the fatty stroma to form a ductal tree; (ii) the onset of pubertal hormonal signaling when terminal end buds form at the tips of the mammary ducts and start to invade the fat pad from which they will extract milk fats, vitamins, and other nutrients (100⇓–102); and (iii) the repeatable hormone-induced cycles of pregnancy, lactation, and weaning (101, 102).

Vitamin D-deficient adipose tissue has been modeled by an adipose-specific Vdr knockout (Vdr KO) mouse, revealing significant and, from an evolutionary perspective, fascinating sex-by-diet effects (99). Female mice that are unable to bioactivate vitamin D in adipocytes exhibit increased visceral adipose tissue overall, but they have a normal mammary fat pad mass. Although the mass of the mammary fat pad is unaffected by the loss of Vdr, there is an increase in the ductal branching density within it (99). Surprisingly, however, this only occurred when the mice are fed a high-fat diet (reminiscent of traditional human diets in Artic cultures) (99). These results suggest that ductal branching is stimulated by adjacent vitamin D-deficient adipose tissue in response to high-fat diets during the hormonally induced stages of mammary development (98, 99, 103). In sharp contrast to the results for female mice, there are no effects in males with the adipose-specific Vdr deletion in either the normal or high-fat diet (99).

The results of the Vdr KO study, bolstered by additional research (reviewed in ref. 73), reveal a sex-specific relationship between vitamin D deficiency and adipose tissue that includes an increase in mammary ductal branching during the hormone-induced stages of breast development (Fig. 2B). This physiological increase in ductal branching in response to very low levels of vitamin D indicates that other genetic mechanisms that increase ductal branching may be positively selected for in populations experiencing chronic vitamin D deficiency, such as in UV-poor environments at high latitudes.

A completely independent line of evidence suggests that variation in ductal branching influences milk content (Fig. 2A). Variation in milk vitamin D levels is heritable in cows and shows evidence of positive selection over the short time frames associated with animal domestication, demonstrating that the underlying genetic mechanisms are highly responsive to selective pressure (104, 105). Milk vitamin D content varies by cattle breed, with some producing notably higher concentrations irrespective of UV exposure (reviewed in ref. 106). Functional associations of 31 genes expressed in cows producing different milk protein and fat concentrations include transcriptomes of genes involved in mammary gland bud elongation (107), suggesting that mammary duct branching contributes to the genetic underpinnings of variation in milk composition.

EDAR V370A’s Influence on Ductal Branching

The relationship between vitamin D and ductal branching calls for a renewed consideration of EDAR V370A’s influence on mammary duct branching density as the phenotypic target of selection. As with other ectodermally derived structures (e.g., teeth), the earliest phase of mammary development results from epithelial–mesenchymal tissue interactions at embryonic day 11 in the mouse (100⇓–102). Ductal branching is established via NF-ĸB signaling, a downstream target of the ectodysplasin pathway (108). The increase in branching observed in mouse models of EDAR V370A (30, 31, 44) is likely the result of the EDA/EDAR/NF-ĸB pathway’s influence on the proteins that modulate the strength of the junctions that adhere the myoepithelial cells of the mammary ducts (109⇓–111) (Fig. 2 C–E), facilitating or hindering branching during morphogenesis (112). The more elaborate ductal branching induced by EDAR V370A during the embryonic stage of mammary development likely enhances the physiological effect of the increased branching induced by vitamin D deficiency in the later stages of breast development (during puberty and gestation/lactation).