THE INTERACTION BETWEEN GENETICS AND CLIMATE ON CRANIOFACIAL VARIATION: EXAMINING THE CAUSATIVE FORCES OF MACROMORPHOSCOPIC TRAIT EXPRESSION By Amber Plemons A DISSERTATION Submitted to Michigan State University in partial fulfillment of the requirements for the degree of Anthropology—Doctor of Philosophy 2022 ABSTRACT THE INTERACTION BETWEEN GENETICS AND CLIMATE ON CRANIOFACIAL VARIATION: EXAMINING THE CAUSATIVE FORCES OF MACROMORPHOSCOPIC TRAIT EXPRESSION By Amber Plemons Anthropologists have an extensive history using cranial form to measure group relatedness in past and present populations to answer a range of questions concerning population histories and cultural practices. However, most biological distance studies using skeletal remains do not consider extrinsic forces influencing modern human variation. Researchers have explored evolutionary and plastic responses in cranial form using measurements of the cranium and mandible, but these studies generally drew inferences through population comparisons or using inadequate statistical and biological models that so often lead to conflicting findings or confounding interpretations. To fill this gap in our current understanding of modern human variation, global craniofacial morphological, climatic, and genetic datasets are combined to measure the magnitude and directionality of several climate variables on craniofacial form, while controlling for population structure (e.g., microevolutionary forces and population history). Craniofacial morphological data from the Macromorphoscopic Databank (MaMD) and microsatellite data from Pemberton (2013) represent populations that overlap in geographic space. Finally, climate data from the National Oceanic and Atmospheric Association (NOAA) and Climate Research Unit (CRU) websites were obtained for weather stations in close proximity to populations under study. This combined dataset is used to explore the interaction between climate and genetics on craniofacial variation across 11 geographic regions using a mixed model approach known as Bayesian Sparse Factor Analysis of Genetic Covariance Matrices (BSFG). Data analysis follows the methods described by Katz and colleagues (2016) but expands their study through the exploration of selection processes using additional climate variables, including coldest month and driest month averages and annual ranges of temperature and absolute humidity. Overall, the study found significant correlation between genetic and phenotypic data indicating MMS traits can serve as genetic proxies in biodistance analyses. Several traits had higher heritability estimates (malar tubercle, zygomaticomaxillary suture course, postbregmatic depression and anterior nasal spine). Features associated with the nasal complex and facial breadth, particularly anterior nasal spine, nasal bone contour, and interorbital breadth, had strong associations to climate. These climate findings correspond to previous research on nasal form and environment where cold-dry environments select for high, narrow noses. Further evidence of selective forces in MMS traits are apparent with the reduction of these features in more variable climates where the respiratory system experiences less stress. The evolutionary mechanisms behind craniometric data have been explored extensively. Such studies use a full suite of traits that capture overall size and shape of the human cranium; however, MMS traits focus on macroscopic assessments primarily in the midfacial skeleton. MMS trait data are particularly important for expanding our understanding of natural selection whereby a large portion of cranial evolutionary research has centered around the neutral evolutionary processes. The wealth of research demonstrating the nasal complex is highly responsive to climate due to respiratory stress emphasizes the importance of exploring the proportion of genetics and environments on MMS trait manifestation. This project provides an evolutionary foundation of the neutral evolutionary and selection processes controlling systematic patterns of global craniofacial variation in the MaMD. This dissertation is dedicated to Autumn Allison. Thank you for making this journey full of wonderful memories and worth the hard ones. iv TABLE OF CONTENTS LIST OF TABLES ........................................................................................................................ vii LIST OF FIGURES ..................................................................................................................... viii INTRODUCTION .......................................................................................................................... 1 Project Design ............................................................................................................................. 6 Organization of Chapters .......................................................................................................... 10 BACKGROUND .......................................................................................................................... 12 Theoretical Frameworks for Modes of Human Variation......................................................... 13 Anthropological Genetics and Evolutionary Forces ............................................................. 13 Molecular Markers for Human Genetic Variation ............................................................ 16 Population Structure.......................................................................................................... 18 Phenotypic Plasticity............................................................................................................. 24 Human Phenotypic Variation, Climate, and Genetics .............................................................. 28 Evolution of Body Size and Shape ....................................................................................... 29 Cranial Variation and Plasticity ............................................................................................ 31 Nasal Morphology and Climate ............................................................................................ 35 MATERIALS AND METHODS .................................................................................................. 38 Materials ................................................................................................................................... 38 Macromorphoscopic Traits ................................................................................................... 38 Microsatellite Data ................................................................................................................ 41 Climate Data ......................................................................................................................... 42 Methods..................................................................................................................................... 44 Data Preparation and Summary Statistics ............................................................................. 44 Preliminary Analyses ............................................................................................................ 46 BSFG Model Setup and Genetic Output Information........................................................... 47 Data Analyses for Climate Effects ........................................................................................ 48 RESULTS ..................................................................................................................................... 50 Macromorphoscopic Trait Data ................................................................................................ 50 Summary Genetic Data ............................................................................................................. 55 Genetic Effects on MMS Variation .......................................................................................... 58 Temperature Effects on Craniofacial Morphology ................................................................... 61 Humidity Effects on Craniofacial Morphology ........................................................................ 64 Effects of Climatic Fluctuations on Craniofacial Morphology................................................. 66 Interaction of Genetics and Climate on Craniofacial Morphology ........................................... 71 Genetic Effects ...................................................................................................................... 71 Climate Effects...................................................................................................................... 73 High Heritability Traits ..................................................................................................... 74 Low Heritability Traits ..................................................................................................... 77 v DISCUSSION AND CONCLUSION .......................................................................................... 82 Contextualizing Causative Forces of MMS Trait Expression .................................................. 82 MMS Genetic Relatedness .................................................................................................... 82 Cranial Metric Mixed Models for Genetics and Climate...................................................... 84 Cranial Nonmetric Trait Heritabilities .................................................................................. 85 Limitations and Future Directions ............................................................................................ 87 Dichotomized Traits.............................................................................................................. 87 Regional Population Assignments ........................................................................................ 88 Unmatched Datasets.............................................................................................................. 88 Temporal Variation between Skeletal Collections ............................................................... 89 Broader Evolutionary Significance ....................................................................................... 90 APPENDIX ................................................................................................................................... 94 LITERATURE CITED ................................................................................................................. 98 vi LIST OF TABLES Table 1. Macromorphoscopic traits, their abbreviations, and the range of possible score values (character states). .......................................................................................................................... 40 Table 2. Summary of MMS and genetic populations with sample sizes along with weather station locations. Regional populations refer to the metapopulations assignments that encompass all three lines of data for that region. ................................................................................................. 44 Table 3. Sectioning points for dichotomizing MMS traits. .......................................................... 45 Table 4. MMS trait frequencies for dichotomized scores by Regional Population. ..................... 53 Table 5. MMS population distance matrix. Upper triangle provides Mahalanobis Mean Measure of Divergence distances while lower triangle indicates significant population distances (*) and non-significance (NS). .................................................................................................................. 54 Table 6. Genetic distance matrix for regional populations. Upper triangle provides delta-mu squared Euclidean distances from group means. Lower triangle provides pairwise estimates for the genetic matrix, A. .................................................................................................................... 56 Table 7. Summary climate data used for BSFG analyses. ............................................................ 58 Table 8. Genetic variance-covariance matrix (G) among traits calculated from narrow sense heritability estimates of factors. .................................................................................................... 59 Table 9. Phenotypic (MMS) variance-covariance matrix (P) calculated among traits from broad sense heritability estimates of factors. .......................................................................................... 60 Table 10. Heritability measures of MMS traits. ........................................................................... 61 vii LIST OF FIGURES Figure 1. Distributions of MMS imputation iterations using mice R package (Blue line is true data; Red lines are imputation iterations). .................................................................................... 45 Figure 2. Histogram of missing data distribution for MMS traits. The pattern of missing is depicted in the right plot with gray representing missing and green representing present data ... 51 Figure 3. Hierarchical clustering of regional populations for non-imputed and imputed MMS data. ............................................................................................................................................... 51 Figure 4. Correlation plot of MMS traits. ..................................................................................... 52 Figure 5. 2D multidimensional scaling plot of MMS population distances using ordinal scaling. ....................................................................................................................................................... 55 Figure 6. Procrustes plot of genetic and MMS distance matrices demonstrating group relatedness within and between the datasets. ................................................................................................... 57 Figure 7. Heatmap of traitwise covariance matrix for MMS traits with inclusion of both temperature and genetic effects. Yellow indicates high covariance and blue indicates low. ....... 60 Figure 8. Climate coefficients (black dots) with 95% confidence intervals for coldest month averages. Values right of the red line have positive directionality, and those to the left have negative directionality. .................................................................................................................. 62 Figure 9. Comparison of temperature effects on contrasts of variance for ANS (left) and INA (right). The X-axis plots temperature difference for pairwise regional populations and Y-axis in contrast of variance. ANS is the only trait demonstrating a strong correlation with temperature when population structure is removed. The remainder of the traits resembles INA with weak correlation. .................................................................................................................................... 63 Figure 10. Procrustes plot of group relatedness for MMS Mean Measure of Divergence Distances (red) and contrasts of variance (black) for temperature effects, phenotypic and residual error variance, and population structure. ............................................................................................... 64 Figure 11. Climate coefficients (black dots) with 95% confidence intervals for driest month averages. Values right of the red line have positive directionality, and those to the left have negative directionality. .................................................................................................................. 65 Figure 12. Procrustes plot of group relatedness for MMS Mean Measure of Divergence Distances (red) and contrasts of variance (black) for humidity effects, phenotypic and residual error variance, and population structure. ............................................................................................... 66 Figure 13. Plot of annual monthly temperature averages for (1986-2015) for each regional population. The x-axis is temperature in degrees Celsius from coldest to warmest. The red line indicates 18-degree breaking point for ratio calculations. ............................................................ 67 viii Figure 14. Plot of annual humidity temperature averages for (1986-2015) for each regional population. The x-axis is absolute humidity values from driest to wettest. .................................. 68 Figure 15. Climate coefficients (black dots) with 95% confidence intervals for annual temperature ratios. Values right of the red line have positive directionality, and those to the left have negative directionality. ......................................................................................................... 70 Figure 16. Climate coefficients (black dots) with 95% confidence intervals for annual ranges of absolute humidity. Values right of the red line have positive directionality, and those to the left have negative directionality. ......................................................................................................... 71 ix INTRODUCTION This dissertation examines evolutionary and ecogeographic forces impacting macromorphoscopic (MMS) trait expression. I first assess and measure the causative forces (climate and genetic) of human craniofacial variation through the lens of biological evolutionary theory and the application of statistical frameworks adopted from quantitative genetics. Next, I explore the influence of climatic variables on MMS trait expression to develop new theoretical concepts outlining the role climatic fluctuations play on phenotypic expression. This project will provide an understanding of evolutionary forces controlling systematic patterns of global craniofacial variation, particularly as expressed in the in the Macromorphoscopic Databank (MaMD). Cranial morphology in modern humans is known to be correlated with neutral evolutionary forces and, at the same time, studies have shown that climate influences facial shape in the skeleton (Hubbe et al 2009; Maddux et al. 2017) and in soft tissue morphology (Weiner 1954; Maddux et al. 2017; Zaidi et al. 2017). Thus, underlying confounding influences of plasticity and selection are inconsequential in our evolutionary theoretical works cannot be assume. This project measures the magnitude and directionality of several climate variables on MMS trait data while controlling for population structure obtained from microsatellite data using a large world-wide genetic dataset (Pemberton et al. 2013). A statistical framework developed for quantitative genetics capable of analyzing high-dimensional data while selecting the most biologically-plausible models will quantify the relationship between MMS trait data, population structure, and climate variables. The climate variables used in this project extend beyond previous anthropological studies focused only on cold temperature and build upon studies investigating the relationships between climate, genetics, and skeletal morphology. Based on previous studies and preliminary analyses, there is clear evidence for 1 population variation in the MaMD and for potential ecogeographical patterning in MMS trait expression. However, there is currently no empirical support to explain why these patterns exist. This research will attempt to isolate the distinctive evolutionary forces impacting expression in craniofacial morphology and explore new theoretical questions regarding the interaction of evolutionary and plastic responses in the human facial skeleton. The MaMD and genetic datasets provide a unique opportunity for exploring larger theoretical questions to contribute to conversations regarding evolution and plasticity in anthropology, ecology, biology, and related disciplines. Additionally, quantifying the directionality and magnitude of climate and genetic influences on macromorphoscopic trait expression provides empirical evidence to explain how and why human variation persists, arming us with the educational tools necessary for breaking down longstanding perspectives of biological race. ****** The relationships among past and living peoples are frequently quantified to measure divergence using patterns of human skeletal variation (i.e., biodistance) (Buikstra et al. 1990; Hefner et al. 2016). Biodistance studies are applied to populations at various scales of analysis to understand cultural adaptations through interpretations of postmarital residence, within group social structures, between group interactions, and small-scale migration, as well as to estimate population affinity of unidentified human remains in forensic contexts (Buikstra et al. 1990; Pietrusewsky 2017; Hefner et al. 2016). The term population affinity is an approach adopted by biological anthropologists for examining group similarities and differences between both social and biological groups. These groups can be altered in the level of refinement or nature based on the research goals and questions or case context in forensic investigations, allowing for anthropologists take multiple 2 influences of human variation into consideration when examining biological distance between and within groups or group membership for an individual (Winburn and Algee-Hewitt, 2021). This term is gaining popularity in biological anthropology replacing the outdated, inaccurate terms, race and ancestry. Race was used for decades to correlate human genotypes and phenotypes using visual characteristics and, later, craniometrics, classifying humans primarily into 3 broad races: African, Asian, and European. After anthropologists established a consensus that race is social rather than biological, many shifted to using the term ancestry. Ancestry is viewed as geographic patterns of human variation due largely to gene flow; unfortunately, when this term is used in anthropological research, it often relies on this broad, three-group classification system developed from race studies, meaning there is a new term given without a considerable shift in theory or methodology (Edgar and Pilloud 2021; Winburn and Algee- Hewitt 2021). Anthropologists are now shifting to the use of population affinity to demonstrate a shift in theoretical and methodological approaches to human variation and biodistance research. Continued discussion on how we are thinking about and applying concepts of human variation are of the utmost importance to best understand what these types of data are truly measuring and what we are aiming to achieve from these data, both theoretically and practically. Biodistance and affinity-related studies identify systematic human variation in cranial form; however, much of the research using biodistance analysis to answer questions about skeletal remains relies on comparisons to skeletal traits and/or measurements recorded in other skeletal populations without consideration of the causative forces influencing the observed variation. There is consensus that variation in cranial phenotypes significantly reflects underlying genomic variation resulting from neutral evolutionary forces (mutation, gene flow, and genetic drift) (Reyes-Centeno et al. 2017; Raff 2019; Relethford 2019; Reyes-Centeno and Rathmann 3 2020); therefore, anthropological studies of biodistance often rely on patterns of metric and nonmetric variation as a proxy for genetic relationships between and within groups of people. However, phenotypic plasticity and natural selection are confounding forces affecting human physical form at both developmental and evolutionary scales. Phenotypic plasticity refers to the ability of organisms to change their phenotypes in response to the local environment (Roberts 2011; Garland and Kelly 2018; Perry et al. 2018). Plasticity is important in human evolution and may even increase the speed of evolutionary change via selection mechanisms rather than hinder it, particularly when an individual moves into novel environments (Perry et al. 2018; Kuzawa and Bragg 2012). While general rules about the relationship between climate and body size and shape are widely accepted in anthropology (e.g., Bergmann’s Rule and Allen’s Rule), very little research has been achieved in detecting the correlation between climate and cranial shape and size, particularly regarding craniofacial form. Several foundational anthropological studies have debated over the role of climate in shaping cranial form by inferring evolutionary and plastic responses of cranial form via population comparisons of cephalic indices across generations and geographic birthplaces (Boas, 1911, 1912, 1940; Gravlee et al., 2003a, 2003b; Sparks and Jantz, 2003, 2003b). The authors of these studies also debate over the role in which genetics influences variation between these study groups and whether plasticity is significant enough to conceal those genetic influences in phenotype. However, due to the lack of genetic data in these studies, there is no way to definitively claim the degree of either genetic or environmental influence (Relethford 2004). Very few studies to date specifically test for interactions between cranial form, genetics, and climate using comprehensive models. 4 While population affinity metrics have been used extensively and deemed highly useful in anthropological research, there is little understanding on the evolutionary forces responsible for patterns of variation in human craniofacial form. The utility of these metrics is now heavily debated, due in large part to the lack of evolutionary backing (Bethard and DiGangi 2020; Stull et al. 2020; DiGangi and Bethard 2021). In fact, some forensic anthropologists argue for the complete discontinuation of population affinity metrics in forensic investigations, arguing that the lack of an evolutionary framework results in a reliance on traditional racial typologies developed during European colonialism and propagated as part of the eugenics movement (Bethard and DiGangi 2020; DiGangi and Bethard 2021). However, a greater majority consider that population affinity is critical to establish a biological profile for identifying unknown persons and is the pinnacle component to all other aspects of the biological profile (age, sex, and stature) that are inter-dependent due to ecogeographic variation (Stull et al. 2020). In other words, because variation in age, sex, and stature is population-specific, it cannot be disassociated from population affinity (Winburn and Algee-Hewitt 2021). Anthropologists are being more cognizant in their estimations of biodistance and population affinity by consciously incorporating considerations of human variation and evolutionary mechanisms behind phenotypic variation in our analyses (Ross and Pilloud 2021). To resolve current debates on the utility of population affinity metrics and satisfy callouts of theoretical standards, it is necessary to establish an evolutionary framework for understanding population variation in human cranial form. Recently established, genetic databases and opensource climate databases provide new opportunities to build complex models for exploring the relationship between the genome, climate, and cranial form. The recent development of the Macromorphoscopic Databank (MaMD)—a diverse, global dataset of macromorphoscopic (MMS) trait data—permits the 5 reliable application of craniofacial morphology as a measure of relatedness. Macromorphoscopic traits are a series of quasi-continuous cranial morphological traits used for estimating population affinity in forensic contexts (Hefner 2018). Until recently, the use of macroscopic craniofacial traits in biodistance analyses had no evolutionary basis due to the lack of adequate reference samples to investigate these relationships. Hefner (2009, 2018) developed standard data collection procedures and reference data (MaMD) for recording and analyzing craniofacial morphology. These standardized procedures allowed for the development of this large reference dataset with reliable trait data, making it possible to investigate global cranial morphological variation in a meaningful way. Project Design This project utilizes a combined dataset of microsatellite, macromorphoscopic trait, and climate information which correspond in geographic space. These data are used to approximate the causative forces behind global, regional human craniofacial variation. The directionality and magnitude effects of climate variables on craniofacial morphology will be measured, controlling for population structure of modern humans. However, by nature of the models applied, important information for genetic effects on MMS expression will also be acquired. The hypotheses and research aims will build upon each other to address the overarching research question: Do genetics and climate influence MMS trait expression and, if so, to what degree? Results from this study will help reveal if there is an evolutionary foundation for MMS traits and potentially unravel various mechanisms by which these forces are controlling craniofacial morphology. Research Goal 1: The first goal is to assess the relationship between temperature and MMS trait expression by estimating direction and magnitude of regression coefficients for coldest month averages. 6 Specific Aim: BSFG models are used to calculate the climate coefficient from average coldest month temperature values to assess the direction and magnitude of the effect on cranial morphology. This model removes the population structure from the effects to assess actual correlation more accurately between these data. The BSFG method has been applied to craniometric data focusing on the correlation of coldest month temperatures with cranial measurements (Katz et al, 2016). This phase of the project replicates previous research by Katz and colleagues aiming to quantify effects of coldest month temperatures on cranial MMS traits. Although several studies explore the relationship between cranial and soft tissue morphology and climate (Weiner 1954; Relethford 1994; Roseman, 2004; Roseman and Weaver, 2004; Katz et al. 2016; Maddux et al. 2017; Zaidi et al. 2017), I will expand this research by testing macroscopic morphological features in a high dimensional model for the first time and using larger sample sizes. The BSFG will be modeled from three lines of data: 1) MMS trait scores from the MaMD (Hefner 2018), 2) microsatellite loci from Pemberton and colleagues (2013), and 3) temperature from NOAA (National Oceanic and Atmospheric Association) online datasets (https://www.ncdc.noaa.gov/cdo-web/datasets). These datasets are described below. Combining these lines of data in a robust mixed model provides the statistical framework needed to disentangle the notoriously complex relationships between climate, genetics, and cranial form. Katz and colleagues (2016) found populations in colder climates tend toward larger crania, particularly for cranial vault breadth measures, despite accounting for population structure. Additionally, soft tissue studies have reported significant correlations between various nasal traits and climate (Zaidi et al., 2017). Therefore, similar directionality in MMS trait expression due to average coldest temperature data is expected regardless of population structure, 7 specifically for the breadth traits (i.e., Nasal Aperture Breadth and Interorbital Breadth) and traits associated with nasal morphology and the overall size of the nasal cavity (i.e., Nasal Bone Shape, Anterior Nasal spine, and Inferior Nasal Aperture). We know, for instance, that on average nasal breadth is narrower in European populations compared to African populations (Hefner, 2009). This corresponds with other studies suggesting cooler climates select for narrower nasal cavities related to cold air inhalation (Yokley, 2009; Noback et al., 2011; Holton et al., 2012; Maddux et al., 2016). The correspondence between these findings suggest predictable directionality of nasal breadth with temperature; however, the magnitude of cold temperatures varies across studies. I hypothesize other MMS traits will exhibit greater genetic effects and weaker temperature correlation with unreliable directionality (e.g., orbit shape and postbregmatic depression). These are traits that are known to have strong geographic patterning but their function in not clear (Hefner and Linde, 2018). There is, of course, potential for some traits to have both weak genetic and climate effects which would indicate a trait as having no evolutionary basis. Research Goal 2: For the second goal, I aim to determine the relationship between humidity and MMS trait expression by estimating direction and magnitude of regression coefficients for driest month averages. Specific Aim: I expand on Katz and colleagues’ (2016) and the first research goals of this study on temperature effects to include the effects of humidity on craniofacial variation. There is consensus that relative humidity does not significantly influence nasal soft tissue morphology (Maddux et al., 2016; Zaidi et al., 2017), so those data will not be used in this study. However, there has been significant influence of absolute humidity on facial morphology documented in soft tissue morphology of the nose and internal nasal structures (Yokley, 2009; Noback et al., 8 2011; Zaidi et al., 2017), but the magnitude of these effects are not quantified in the facial skeleton. Maddux and colleagues (2016) note that respiratory pressures that shape the nasal complex are highly influenced by absolute humidity regardless of temperature; therefore, the expectation is for absolute humidity predictors to exhibit significant magnitude effects on MMS trait expressions involving the nasal complex. Research Goal 3: The final goal is to explore if high degrees of climatic annual fluctuation variation demonstrate intermediate expression of MMS traits. Specific Aim: Measuring the relationships between MMS trait expression, population structure, and various climate variables adds to the current conversation on plasticity and evolutionary responses in the facial skeleton. For example, the evolution of adaptation concept has been neither fully explored nor applied to biological anthropological research. Richard Potts presented the notion of variability selection (VS), which he defines as the “process considered to link adaptive change to large degrees of environment variability” (Potts, 1999; Potts, 2013; Potts and Faith, 2015). The VS theory contends that great fluctuations in the environment result in increased survivability through the differential expression of more adaptive traits. Therefore, he agrees that natural selection is geared to survival in specific environments, but his theory diverges from others emphasizing flexibility in trait selection in fluctuating environments. This aspect of morphology will be explored using the difference between number of warm to cold month averages for temperature variation and the range of driest and wettest month averages for humidity variation. These difference values saliently capture experienced temperature and humidity per geographic region to determine if higher degrees of variation in climatic experience results in intermediate traits to facilitate survival. This intermediate expression would allow the population to be more adaptive to climatic fluctuations for increased survivability. However, this 9 pattern is expected to be expressed only in MMS traits found to have high climatic influence despite population structure. Those in stable climates should theoretically reflect directionality of trait expression observed in the previous temperature and humidity models. The interaction of evolutionary and developmental plastic responses on climate change is not well understood but is currently a keen area of ecological research (Kelly, 2019). The MaMD, in conjunction with global genetic information, provides a unique opportunity to examine theoretical questions concerning the evolution of human variation. The MaMD does not have the temporal refinement one may hope for to explore such questions, but the level of spatial coverage provided by the MaMD may suffice. Ecologists have studied the use of space-for-time substitutions. Those studies indicate space is a reliable alternative for time when studying biodiversity if there is concordance in scale between temporal and spatial variation (Blois et al., 2013). Therefore, the MaMD is likely the most appropriate dataset for examining relationships between cranial morphology and climate due in large part to the refinement and geographic coverage of ancestry and place of origin information for individuals in the database. Organization of Chapters Beyond this introduction chapter, which presents project goals and the significance of my research, chapter two provide a literature review for pertinent background needed to understand the project design and findings. This chapter includes two subsections. First, Theoretical Frameworks for Modes of Human Variation discusses in detail the ways human genetic variation is created through each evolutionary mechanism, as well as the role of plasticity on human phenotypic variation. The second section, Human Phenotypic Variation, Climate, and Genetics, presents a series of anthropological and human biological studies on the relationship between these three variables and how these studies assisted in generating questions and research aims for this dissertation. 10 Chapter 3 includes a discussion of the materials and methods. Special attention is given to sample collection for each line of data (genetic, climate, and MMS), statistical procedures, and detailed description of how the BSFG model and data are structured. The results and findings of the analytical methods are presented in Chapter 4. Results are organized by preliminary and summary data followed by results for each research question and analytical method. Chapter 5 is a discussion on the implications of the results from this dissertation research and how these findings can be expanded upon to gain a more nuanced understanding of the evolutionary and ecogeographic forces of human craniofacial morphology. Finally, conclusions of what is learned, and future research goals are summarized in Chapter 6. 11 BACKGROUND This chapter is divided into two overarching topics: the theoretical frameworks for modes of human variation and the relationship between human phenotypic variation, climate, and genetics. The first section provides a brief review of the field of anthropological genetics regarding how genetic information has been applied to anthropological questions and how theoretical applications have changed with advances in genetic technology. Following a brief disciplinary history, the different molecular markers used to detect variation in the human genome are detailed, in addition to how these markers are used to estimate population structure and reveal population histories. A final discussion is provided on phenotypic plasticity and the role plasticity plays in shaping human variation alongside genetics. Human Phenotypic Variation, Climate, and Genetics comprises three subsections: Evolution of Body Size and Shape, Cranial Variation and Plasticity, and Nasal Morphology and Climate, and illuminates how the larger evolutionary concepts have been applied in current literature. The primary focus of this dissertation is to determine the magnitude and directionality of climatic influences on craniofacial morphological expressions. Therefore, I review current literature regarding the relationship between the climate and the human body for both skeletal and soft tissue anatomy, with specific attention given to cranial form. A general lack of consensus concerning the magnitude of plasticity’s impact on human cranial form warrants a review of the literature on established relationships between climate and body shape and size. These known correlations in body form and climate demonstrate a reasonable degree of plasticity in human body, indicating a great likelihood that climate also significantly influences cranial form. Finally, studies are presented demonstrating correlations between soft tissue morphology 12 in the human face and climate which lead to predictable patterns of human variation; a pattern we also expect in the craniofacial skeleton. Theoretical Frameworks for Modes of Human Variation Human variation can be observed in many biological features and elements, but for the purposes of this review section, I focus primarily on cranial and facial phenotypic variation. Phenotypic variation is attributed to environmental forces due to developmental plasticity and to genetic variation caused by evolutionary mechanism (e.g., mutation, gene flow, genetic drift, and natural selection) (West-Eberhard 2008). Developmental plasticity refers to the ability of organisms to change their phenotype in response to the local environment which occurs during the ontogeny of a trait within individuals (West-Eberhard 2003). The relationship between plasticity and evolutionary forces is heavily researched, but that relationship is not necessarily well understood. We have only limited knowledge concerning how the environment and plasticity contribute to (or hinder) the rate of human adaptation. Likewise, the impacts of environmental sensitivity on phenotypes are heavily debated and criticized (West-Eberhard 1989, 2021; Perry et al. 2018; Fox et al. 2019). Anthropological research traditionally accepted that phenotype reflects genotype (Relethford 2019), but there is significant ecological evidence demonstrating that plasticity is a critical component of Darwinian evolution, significantly influencing phenotypic variation (West-Eberhard 2021). Anthropological Genetics and Evolutionary Forces Anthropological genetics applies analysis of the human genome and concepts of genetic variation drawn from population genetics to answer questions regarding human evolution and human variation, the diaspora of anatomically modern humans out of Africa, and the biocultural interaction between humans and disease (Crawford 2007; Crawford 2019). Population genetics 13 employs processes of natural selection, genetic drift, mutation, and gene flow to study genetic compositions of populations, while anthropological genetics applies those processes of genetic variation to understand human history, behavior, and modern variation. The linear progression of the field of anthropological genetics parallels human genetics (Crawford 2019). The earliest studies applying genetic information began using ABO blood groups to distinguish between human ‘races’ (Marks 2012; Crawford 2019). Later, during the 1950s and 1960s, much of the research by geneticists and physicians contained anthropological components focusing on heritability of diseases and the influences of inbreeding on genetic drift, particularly in remote, or non-Westernized, societies (O’Rourke 2003; Crawford 2019). These early studies were smaller in scale, but the shift from race studies to variation between regional, local, and cultural populations was underway. In 1969, the term anthropological genetics was coined by Derek F. Roberts, a British anatomist who trained several prominent anthropological geneticists and was responsible for much of the foundational work in the field (Crawford 2018, 2019). Several other figures— such as Frank Livingston, Eugene Giles, and Richard Lewontin—were instrumental in initiating anthropological genetics. These researchers combined methods and concepts developed by geneticists and physicians with cultural data about human populations to better understand human variation and its relationship to the genome (Crawford 2007). For example, Livingston (1971) examined the effects of agricultural practices on malaria prevalence rates. (Crawford 2018). Giles and colleagues (1970) studied fluctuations in gene frequencies in populations in New Guinea to examine genetic drift in small populations. Lewontin, a population geneticist and respected figure in the discipline, explored 15 protein loci to measure genetic variation within and between human populations, historically reporting that 85% of genetic variation exists 14 within populations and only 15% exists between populations (Lewontin 1972). This idea—there is more variation within populations than between populations—is foundational in anthropological theory today (Marks 2012). This new disciplinary pathway made it possible for researchers to be trained extensively in both anthropology and genetics (Crawford 2018). The discovery of the structure of DNA by Watson, Crick, Franklin, and Wilkins in 1953 and DNA sequencing by Sanger in 1977 marked eras of dramatic change in the anthropological research (Crawford 2019). Now anthropological geneticists could explore human evolution in more meaningful ways than with simple blood groups and protein markers alone (Marks 2009). Approaches like genetic mapping and quantitative genetics permitted reconstruction of sex-specific migration patterns and the use of molecular genetic markers, such as Short Tandem Repeats (STRs) and Single Nucleotide Polymorphisms (SNPs), to better understand variation in human phenotypes (O’Rourke 2003; Crawford 2007). The two main areas of research in anthropological genetics are population structure and population history. These studies dominated the field from the late 1960’s to early 1980’s and relied heavily on DNA markers to explore the influence of mate selection and population subdivision on genetic variation (population structure) and human origins and migration (population history) (Relethford 2019). In opposition to earlier works using blood groups and other racially driven genetic studies, this research attempted to identify clusters of populations based on allelic frequencies at various loci with no preliminary knowledge or assumptions about those groups (Crawford 2007; Relethford 2019). Over the last several decades, accessibility to full human genomes from around the world has made it possible examine genetic variation within and between populations in greater detail (e.g., loci-specific analyses of differential allelic expressions). These refined levels of genetic analyses, along with 15 the application of more robust computational capabilities, have allowed anthropologists to ask longstanding questions in novel ways, as well as expand the breadth of their research in new ways. We can expect the field of anthropological genetics to only continue to improve as our understanding and methods of genetic analyses continue to advance (Gokcumen 2018). Molecular Markers for Human Genetic Variation There are several molecular markers used in anthropological genetic studies depending on the goals of the researchers. Most of the markers used in studies of population structure and history are in the autosomal regions, but some studies focus on sex chromosomes (i.e., Y chromosome and mitochondrial DNA, or mtDNA) to trace parental lineages. Autosomal markers are generally preferred in anthropological genetics as they allow for examination of genetic variation in the whole population rather the sex-specific variation (Rubicz et al. 2007). These autosomal markers can exhibit variation at single nucleotides or along segments of DNA (Rubicz et al. 2007; Stoneking 2016). Human genetic research on single nucleotide variations is known as SNPS. SNPs occur from mutations or during DNA replication and exhibit a relatively slow rate of evolution. Due to that slow rate, SNPs have not changed among modern human groups and can illuminate common ancestry between individuals sharing the same marker (Rubicz et al. 2007). SNP analysis is a fast and cheap option to examine genetic variation (Stoneking 2016). SNPs are also useful when samples are degraded since they require fewer base pairs than short tandem repeats (STRs) (Dumache et al. 2016). Tandem repeats are repetitive DNA sequences at specific loci, which mutate at a much faster rate than SNPs. The two types of tandem repeats used in anthropological genetics include microsatellites and minisatellites (Stoneking 2016). Microsatellites, or STRs, are repeated sequences of 1-6 base pairs with generally less than 350 base pairs in total. The mutation rate of 16 microsatellites increases with the length of repeated sequence likely due to the increased instability of longer repeats (Li et al. 2004; Rubicz et al. 2007; Stoneking 2016). The faster mutation rate makes them ideal for recent genetic changes in populations, compared to the SNPs described above which are better for examining genetic events that occurred earlier in human history (Rubicz et al. 2007; Vieira et al. 2016). Minisatellites are like microsatellites but have an even faster rate of mutation affecting 10 to 100 base pairs with up to 1,000 total base pairs. In general, microsatellites are more commonly used in anthropological genetic studies and minisatellites are used in forensic contexts for identification (Rubicz et al. 2007, Stoneking 2016); however, more sensitive microsatellites are beginning to replace minisatellites in forensic investigations (Kowalczyk 2018). STRs of unknown, forensically significant individuals can be uploaded into a national database known as CODIS (Combined DNA Index System) which uses a series of autosomal loci for comparing unknown samples to samples from family members of missing individuals (Dumache et al. 2016). Sex chromosomes are used to trace male (Y chromosome) and female (mtDNA) lineages in human history and to classify haplogroups based on diagnostic mutations. Sex chromosomal markers provide less genetic information than their autosomal counterparts but are becoming increasingly useful in anthropological and genetic research (Stoneking 2016). By examining the portion of the Y chromosome that does not recombine with an X chromosome, the paternal line can be examined exclusively. The Y chromosome was believed to be of no use for examining variation for many years until advances in SNPs and STRs research found significant variation for Y chromosomes (Underhill et al. 2000; Stoneking 2016). The ability to distinguish Y chromosomes and trace them through paternal lines is incredibly useful in forensic genealogy, particularly among patrilocal populations (Calafell and Larmuseau 2017). Patrilocality has led to 17 higher rates of female migration, Y chromosomes display greater geographic variability than mtDNA. Y chromosomes have also been useful in population studies of group relatedness and comparison of clans, lineages, and tribes given their use in tracing paternal lineages, degree of variation, and geographic patterning (Kottak 2002; Chaix et al. 2004; Sanchez-Faddeev et al. 2013; Calafell and Larmuseau 2017). Alternatively, mtDNA, are passed down from mother to child and can be used for tracing maternal lineages and determining population origins. This marker has a high mutation rate. It also has greater survivability than nuclear DNA, making it useful for degraded samples such as skeletal remains recovered from archaeological contexts to understand population origin and migration histories of these past populations (Nesheva 2014; Dumache et al. 2016). Y chromosomes and mtDNA can be combined to capture both parental lineages for more accurate and informed population origins (Underhill and Kivislid 2007). Population Structure Population structure is the genetic structure of populations with “deviations from panmixia (random mating) in a population of infinite size, and the effect of these deviations on allele and genotype frequencies” (Relethford 2019: pp.134). In other words, population structure is the measure and distribution of allelic frequencies between groups. In instances of random mating, frequencies of alleles would be similar across groups; however, random mating is uncommon due to assortative mating practices, geographic distance, or other geographic barriers limiting group contact. This leads to differential distribution of alleles across space which can then be measured and used to study population histories across time and space. Caveats to panmixia include inbreeding, finite population size, and population subdivision. Finite population size results in genetic drift, or random chance fluctuations in allele frequencies which causes nondirectional generational changes and may encourage inbreeding which can change 18 allelic frequencies in the population (Relethford 2019). Population subdivision refers to the total population, which is required for panmixia, subdivided into continental and regional populations which are then subdivided into local populations (i.e., fragmenting the total population into smaller, finite groups). These subdivisions can occur due to inconsistent environmental regions, geographic barriers, and sociocultural factors, such as religion or language barriers resulting in nonrandom mating (Mielke and Fix 2007). This is important because as population subdivision occurs and population size reduces, the impact on modes of genetic variation increase (Mielke and Fix 2007; Relethford 2019). Genetic drift will undoubtably alter allele frequencies in every generation, but the degree of generational impact will be determinate by the size of the population (Templeton, 2019). This phenomenon is essentially the law of averages where genetic drift can have a substantial impact on allele frequencies in small populations in a short timeframe, but there will unlikely be significant changes in frequencies in large populations (Honnay, 2013; Templeton, 2019). Genetic drift directly impacts population structure by homogenizing finite population with decreased allele frequencies (reducing within population variation), resulting in increased genetic divergences between populations due to less overlap in present alleles. In the absence of migration, mutation is the only way to keep populations undergoing drift from eventually becoming genetically homogenous; thus, given fixation takes longer to occur in larger populations, the population will likely never be homozygous through drift because mutation is occurring generationally (Templeton 2019). The rate of mutation must be greater than or equivalent to the rate of genetic drift or a population becomes at risk of extinction due to homogeneity (Honnay 2013). 19 If a population dramatically reduces in size in a short period of time, many alleles will be removed from the population reducing genetic variation, which is referred to as the bottleneck effect. Genetic variation through drift can also be reduced through the founder effect. The founder effect is when a population is generated from a small group of individuals meaning only a small number of alleles are present when the population originates (Templeton 2019; Wells and Dlugosch 2019). Population geneticists can use what they know about bottleneck and founder effects on genetic variation to unravel population histories given the rate of change is greatest at the front wave. The front wave of these events shows rapid reduction in genetic diversity, a pattern known as allele surfing. Migration and selection also play key roles in genetic evolution and variation leading to changes in population structure. Migration in population genetics describes genetic variation that occurs over a long period of time due to a population introducing new alleles into another population (Milgroom 2003). As populations of organisms move across the landscape and encounter other populations, the migrating population will contribute alleles to the original population through mating practices and alter allele frequencies of the original population. While genetic drift creates genetic divergence between populations, migration decreases this divergence by (re)introducing alleles into the populations. Without migration and gene flow, mutation and selection would continue to increase genetic divergence (Milgroom, 2003; Templeton 2019). Population geneticists can learn about specific organism behavior and patterns of migration by targeting specific lines of inheritance in the genome. For example, human populations have been documented to migrate based on cultural practices of patrilocality or matrilocality. In patrilocal societies, females will demonstrate greater dispersal while matrilocality results in greater dispersal of males, and these patterns of gene flow can be detected by observing either mtDNA 20 or Y-DNA. Matrilocal societies will exhibit a smaller Fst, proportion of total genetic variance, for Y-DNA than mtDNA and vice versa in patrilocal societies (Templeton 2019). Migration can be estimated directly by tracking movement of specific alleles or multilocus genotypes and indirectly by examining the degree of genetic variation between populations (Milgroom 2003). Through selection, alleles present in the gene pool through mutation, recombination, genetic drift, and migration are edited to increase overall fitness of the population (Hamilton 2009; Templeton 2019). Natural selection is the “differential reproduction of organisms as a function of heritable traits that influence adaptation to the environment” (Carey 2003: pp.202). There are three modes by which natural selection influence genetic variation in polygenic traits: directional, stabilizing, and disruptive (Carey 2003). Directional selection is when a specific phenotype is more advantageous to the local environment, causing increased survivability over other phenotypes. This results in higher allele frequencies for the advantageous phenotype due to the decrease in the less advantageous, ultimately leading to decreased genetic diversity (Carey 2003; Kingsolver and Diamond 2011). Stabilizing, or balancing, selection refers to the favoring of average phenotypes, such as birth weights, to ensure fitness of the organism. While directional is generally considered the primary mode of selection, stabilizing is the most common mode given most organisms have become stable in their local environments (Carey 2003; Kingsolver and Diamond 2011). Stabilizing selection is believed to eventually decrease genetic variation over time (Carey 2003). The last mode, disruptive selection, opposes stabilizing selection where extreme phenotypes are favored over the average with the average having reduced fitness. This is the lease common mode of selection and does not appear to occur in human evolution as it generally leads to speciation (Carey 2003). Examining events of changes in allele frequencies helps to explore the evolutionary histories of organisms across the landscape and adds to the 21 puzzle of how they adapted over time (Wells and Dlugosch 2019). Therefore, it is critical to understand the processes by which allele frequencies are altered over time and space. Calculations of population structure can estimate the genetic relationship between populations using methods like principal components analysis or distance measures to estimate relatedness (similarity/dissimilarity), known as identity-by-state (IBS) and identity-by-descent (IBD) analyses (Salojärvi 2019). IBS refers to identical alleles or DNA sequences that are not identical due to common ancestry while IBD refers to those that are a result of common ancestry (Templeton 2019). Some studies have created computational models for estimating admixture to examining population structure, but these models are more assumptive than PCA or distance analyses (Salojärvi 2019). Rosenberg and colleagues (2002) conducted a global study of population structure using microsatellite data to measure sample relatedness through cluster analyses. The authors reported the number of genetic clusters of populations (k) as being k = 6, but other studies suggest k was assigned by the authors rather than determined by the data. Although many assumed there were six genetic populations following Rosenberg and colleagues’ findings, Serre and Pääblo (2004) reanalyzed subsets of their data suggesting the number genetic clusters was dependent upon the k selected. This phenomenon demonstrates continuity in human genetic variation and fluidity in the number of defined populations depending on which loci are used in the analyses. Therefore, we cannot define discrete human groups; however, allelic frequencies are predictable across geographic space, making it possible to approximate geographic origins for most individuals from population structures (Serre and Pääblo 2004; Madrigal and Barbujani 2007). 22 Population Histories Anthropological geneticists are particularly interested in using genetic information for revealing the origin of anatomically modern humans and other significant human migration events like the African diaspora and Peopling of the Americas (Jackson and Borgelin 2010; Reyes-Centeno et al. 2014; Rasmussen 2015; Skoglund et al. 2015; Spradley 2006; Moreno- Mayar 2018a, 2018b; Montinaro 2021). The advances in molecular genetics coupled with the development of worldwide genetic datasets greatly improve the ability to map human population history (Relethford 2019; Templeton, 2019). We must be cognizant that different types of markers used to recreate human population histories and understand patterning in human variation can create different narratives. For example, combined genome and skeletal data has indicated two dispersals of modern humans from Africa, first into Asia around the Middle Pleistocene and later into northern Eurasia (Reyes- Centeno et al. 2014), while mtDNA and Y chromosome studies suggest a single dispersal from Africa (Reyes-Centeno et al. 2017a). Confoundingly, the type of skeletal information used in population history studies is important to consider. The cranium is one of the most frequently used skeletal elements in these studies due to the abundance of available phenotypic information, but not all cranial features preserve population history signatures to the same degree (Havarti and Weaver 2006; Reyes-Centeno et al. 2017b). Specifically, the temporal and facial bones are found to have the highest degree of preservation based the association between geometric morphometric and Fst distances from SNPs (Reyes-Centeno et al. 2017b). Similarly, dentition, another frequently used skeletal element, has been equated to genetic relationships due to their high heritability and assumed neutral selection. However, comparisons of dental metric and nonmetric data with both SNPs and STRs demonstrate that not all dental phenotypes are 23 selectively neutral (Rathmann et al. 2017). Rathmann and colleagues (2017) pair the two types of genetic data (SNPs and STRs) following recommendations that both neutral genomic data types should be compared to phenotypic variation given differential mutation rates and apportionment of genetic variation (Reyes-Centeno et al. 2016, 2017a, 2017b). Again, this highlights the necessity to carefully select the types of skeletal and genetic data used for examining the complexities of human population histories. While genetic information is incredibly valuable, there remains a need to combine these data with other biological variables capturing phenotypic variation, such as cranial metric and cranial nonmetric traits (Cheverud 1988; Carson 2006; Raff 2019; Relethford 2019). Combining these lines of data allow us to consider genetic variation and phenotypic plasticity simultaneously and examine the interaction of biology and the environment in evolutionary contexts. Phenotypic Plasticity The first use of the term plasticity regarding the human species I was able to locate was in Lasker’s book, The Evolution of Man (1961), where he defined plasticity as the “the capacity of the individual to change in response to his environment’ (Lasker 1961:206). He viewed plasticity as permanent adaptive change occurring during an individual’s lifetime. Early studies of plasticity were overly simplistic and typically examined only one environmental variable per study. Over time, the complexity of these studies increased as our understanding of environmental factors influencing the phenotype improved, as well as our understanding in how changes in phenotype influence evolution, regardless of genetics (West-Eberhard 1989; Fuentes 2010). Variation in phenotypic expressions can result from several environmental causes, such as temperature, humidity, altitude, and nutrition, and can relate to overall body size and ratios, body composition, soft tissue morphology, and skeletal structures (Roberts 2011; Agarwal 2016; 24 Suzuki et al. 2020). Lasker (1969) identified three types of environmental adaptation: 1) natural selection of genotypes (selection), 2) behavioral or physiological responses in adults that are reversible in response to the immediate environment (acclimatization), and 3) phenotype modification during ontogeny with no genetic effects (phenotypic plasticity). However, the role of phenotypic plasticity in adaptation and evolution through natural selection is quite complex and the degree of genetic influence on plasticity is not universally agreed upon (Futuyma 2021; Levis and Pfenning 2021; Pfenning 2021). Until the 21st century, evolutionary research focused heavily on genes with very little attention given to phenotypes. This is shifting in evolutionary biological studies which are revealing the significance of plasticity in human evolution, particularly developmental plasticity (West-Eberhard 2021). Developmental plasticity was coined in 2003 by West-Eberhard. This term refers to permanent behavioral or physiological changes in the phenotype due to internal and/or external environmental factors during ontogeny of a trait (West-Eberhard 2003, 2008). Suzuki and colleagues (2020:404) describe the process of developmental plasticity in three phases: “the organism senses environmental cues (sensors), integrates these inputs (integrators), and modulates effectors that alter developmental trajectories (effectors)”. These sensors can be neural or epigenetic sensors and integrated into the body through the physiological processes by which the central nervous and endocrine systems interact. Effectors can regulate growth or pattern at the cellular and tissue level. In other words, the central nervous system integrates environmental signals which control developmental processes by stimulating signaling pathways in the endocrine system that determine gene expression, growth, and morphogenesis in specific tissues throughout the body (Suzuki et al. 2020). Physiological responses to these intrinsic and extrinsic forces can lead to specific phenotypic expressions which can occur at any time throughout 25 development including embryogenesis, fetal growth, postnatal growth, and adolescence (Fausto- Sterling 2005; Agarwal and Beauchesne 2011, Agarwal 2016). These developmentally plastic responses—defined by Lasker’s (1969) third category (phenotype modification during ontogeny with no genetic effects)—have captured the interest of anthropological researchers aiming to understand how the environment influences lasting biological affects through adaptation. Plasticity, particularly cranial plasticity, has been grossly under researched in anthropology due to inconsistent terminology and definitions of ‘plasticity’ coupled with difficulty in measuring non-genetic adaptations (Schell 2011). Genetic adaptation studies develop relatively straightforward models by measuring generational changes in allele frequencies. However, this is much more difficult for non-genetic adaptations which requires measuring the benefit and success of adaptive responses. The two modes of environmental adaptation also create a problem for the scale of measurement (e.g., individual or population). Acclimatization, or physiological adaptation, is more easily examined on the individual level and may (or may not) be visible between generations due to their reversible nature (Schell 2011). Anthropologists and evolutionary biologists are typically interested in population level analyses for answering larger human evolutionary questions. This leads to an emphasis on extrinsic forces of adaptation, such as climate, rather than intrinsic (e.g., heart conditions, V02 max), because we are interested in what influences the human condition and form. Some researchers believe external environmental changes lead to more prevalent plastic responses than do internal changes (Garland and Kelly 2018), meaning that anthropological focus on extrinsic forces will still provide ample information for the exploration of human evolution and variation (Schell 2011). 26 The general lack of understanding of phenotypic plasticity—how it works, and how it can be measured—has caused a significant unbalance in perspectives concerning human adaptation. Biological and social scientists across disciplines agree that phenotypic plasticity plays a major role in human evolution and adaptation, yet we often rely on genetics to explain evolutionary processes and patterns (West-Eberhard 2003). There continues to be a lack of consensus about the evolutionary significance of plasticity due to debates on whether plasticity is adaptive or non- adaptive and on the heritability of plastic phenotypes (Bonduriansky 2021; Chenard and Duckworth 2021; Futuyma 2021; Schlichting 2021; West-Eberhard 2021). While some believe that plasticity does lead to adaptive responses across generations, some argue that there are also non-adaptive responses (i.e., not beneficial but rather a byproduct of reproductive stress) or that plasticity may even retard adaptation, increasing generational stress (Bonduriansky 2021; Fox et al. 2019). However, other studies argue plasticity may increase rather than hinder evolutionary change via selection mechanisms, particularly when an individual moves into a new environment (Kuzawa and Bragg 2012; Perry et al. 2018). There is clearly a lot that is still not understood about plasticity in the human species. Plasticity studies with carefully constructed research designs to examine intergenerational phenotypic responses as well as studies combining theoretical concepts of plasticity with population structure will begin to illuminate some of these heavily debated concerns. When exploring craniofacial plasticity on a global scale, we are examining plastic structural features established during ontogeny and carried intergenerationally through epigenetics. If structural facial features are correlated to climate variables across space and human groups, then we can assume these phenotypic environmental responses were adaptive and selected for over generations. Our scope is limited given the nature of skeletal collections 27 whereby we have no control or definitive data for all intrinsic and extrinsic forces which may be influencing the structural, behavioral, biochemical, and metabolic response. However, some behavioral factors are also likely captured within the structural features through dietary habits influencing muscular structures or other cultural practices, such as artificial cranial modification. Additional analyses could also be performed to incorporate biochemical information, but these would likely be destructive methods which lead to a host of ethical considerations. Finally, there are recent efforts to utilize more skeletal biological approaches with anthropological research aiming to consider the interactions between metabolic and skeletal systems (Devlin 2015; Ives and Humphrey 2018; Snoddy et al. 2018; Brickley et al. 2020), as well as hormonal-skeletal interactions (Gosman et al. 2011; Western and Bekvalac 2017; Berger et al. 2021). By incorporating these various lines of information, we can lean into a more complete approach to skeletal developmental plasticity which includes considerations of both intrinsic and extrinsic forces in future research. Interdisciplinary research engaging with literature from evolutionary biology, human biology, genetics, ecology, and anthropology and combining theoretical concepts from all fields is key to future research in skeletal phenotypic variation. The next section explores literature from these various disciplines to situate this dissertation within current conversations. Human Phenotypic Variation, Climate, and Genetics Below, I highlight three major areas of research applying evolutionary theoretical concepts, combining population structure and histories with environmental factors to explore human evolution and variation. These areas include human body size and shape, cranial variation and plasticity, and internal and external nasal morphologies. Special attention is given to the 28 significant findings for the relationships between human form and climate to provide necessary background information for the project design presented in the next chapter. Evolution of Body Size and Shape The relationship between climate and body size and shape has been extensively explored in many species, including humans. These studies aimed to answer a range of questions about how and why our bodies, as well as the bodies of other mammals, adapt and evolve. Anthropologists typically center their hypotheses and research questions about patterns of body size and shape around two widely accepted principles: Bergmann’s Rule and Allen’s Rule. These principles describe different thermoregulatory mechanism in the body. Bergmann’s Rule describes an ecogeographic phenomenon where larger animal species are correlated with colder climates and smaller species with warmer climates (Bergmann 1847; Mayr 1956). This ‘rule’ is typically quantified using body size measures, specifically height and weight. Allen’s Rule states that animal species in colder climates tend to have shorter appendages than those in warm climates (Allen 1877) and is studied using limb element length measures. Both rules have been documented in numerous species expressing these conditions latitudinally; organisms increase in size and decrease in limb length as latitude increases. However, body size also increases with altitude (Teplitsky and Millien 2014). The causative forces of these patterns have been heavily debated. Researchers argue body size and shape can be influenced by a number of variables, including temperature for heat conservation, food availability, and gene flow (Cushman et al. 1993; Blackburn et al. 1999; McNab 2010; Yom-Tov and Geffen 2011). Teplitsky and Millien (2014) suggest body size is mostly likely affected by a number of correlated environmental factors. Regardless of the causative forces, these patterns have been documented and widely accepted for the human species in Eurasian and African (Ruff 1994, 2002; Holliday 1997). 29 Anthropological studies have noted the increase in human body size with latitude as an adaptive strategy to combat heat loss in colder environments (Roberts 1978; Ruff 1994, Katzmarzyl and Leonard 1998; Auerbach 2012). General temporal trends are apparent and occur in three dominant evolutionary shifts: 1) increase in body size of early Homo, 2) increase in body size around 50kya, particularly in higher latitudes, and 3) decline in average body size beginning ca. 50kya, again particularly in higher latitudes (Ruff 2002). As a result, modern humans in higher latitudes are closer in size to their Pleistocene ancestors while humans in lower latitudes are smaller. The third shift corresponds with other biometric studies using estimated stature and cranial size in which decrease among modern living humans in higher latitudes (Ruff 2002). However, Teplitsky and Millien (2014) reviewed the Bergmann’s Rule literature to determine if climate warming will result in decreased body sized in a current population. The authors found only weak evidence for any recent decrease in size across the research and very few of those studies examined the evolutionary and adaptive influence on physical traits. General body proportions have changed over time to shift remnants of arboreal body functionality in our early hominin ancestors to the more modern form. Early hominins had longer upper limbs than lower limbs and larger relative body breadth to stature (Ruff 2002). Modern humans show clinal patterns in body form with those in colder climates having shorter limbs and wider bodies than those in warm climates (Roberts 1978; Trinkaus 1981; Ruff 1994; Holliday 1997). Ruff (1991) argued body shape variables have the same adaptive responses as size variables since he believes both relate to heat conservation. Thus, body breadth is strongly correlated with latitude which directly affects body surface area, while stature is not correlated with latitude. 30 Several studies have attempted to determine the causative forces of human body size and shape to understand the evolutionary forces acting on postcranial variation. Teplitsky and colleagues (2008) studied a wild bird population to determine whether genetic or environmental forces were responsible for variation in bird body size. They argue body size is related more to plasticity than evolutionary forces identifying no significant genetic influences on postcranial morphology (Teplitsky et al. 2008). Like the cranial studies described above, Roseman and Auerbach (2015) used a mixed model approach controlling for population structure to explore environmental influences on human body size and shape variables traditionally used in studies to identify examples of Bergmann’s and Allen’s Rules. The authors found that population structure significantly influences human postcranial variation and cannot be reduced to clinal variation, arguing previous studies on the ecogeographic principles are insufficient because they do not account for genetic contributions (Roseman and Auerbach 2015). These studies demonstrate the conflicting perspectives on causative forces for body size variation, which leaves these questions unanswered. The literature for both cranial and postcranial variation needs to be considered to answer questions on human plasticity and evolution; this is particularly true given the complexity and currently ambiguous nature of these larger anthropological questions. Cranial Variation and Plasticity Franz Boas conducted the earliest research on plasticity relating to human cranial variation in anthropology (Boas 1911, 1912). His research stems from previous works by Gould, Baxter, and Bowditch on changes in stature of American-born and European-born citizens, as well as Fishberg’s work comparing stature and cephalic index (cranial breadth x cranial length) between Jewish immigrants in America and Europe (Roberts 2011). These studies examined the implications of people moving into novel environments on their and their offspring’s body form. 31 Boas compared cephalic indices between foreign-born and U.S.-born children, as well as between immigrant parents and their children. These studies found that U.S.-born children have significantly different cephalic indices from their foreign-born parents, leading Boas to conclude that the American environment was influencing size and shape of developing crania, paralleling concepts of developmental plasticity generated decades later (Boas 1911, 1912). Based on his findings, Boas argued fixed races cannot be defined by cranial form given crania change in response to local environments, departing from race theory during that time. However, the study was highly criticized for decades as statistically and biologically invalid. Between 1940 and 1960 subsequent studies examined the influence of immigration on stature, many demonstrating that immigrant groups exhibit greater stature and body size than their counterparts in the countries of origin. This was attributed to high caloric diets and food surpluses in the U.S. rather than an enriched climate or better quality of life (Lasker 1969). Interestingly, the early 21st century saw Boas’s data revisited, and renewed attention was given to cranial plasticity, despite the heavy focus on evolutionary research and human population variation in anthropology. Boas’ data were more recently analyzed by two separate research teams applying modern and more robust statistical procedures to measure the influence of plasticity on cranial form. The two teams drew diametrically opposed conclusions. Gravlee and colleagues (2003a) argue Boas was generally correct. However, they found weaker correlations than Boas, suggesting the magnitude of environmental forces may not be as strong as originally considered. Sparks and Jantz (2003) also found weaker correlations in the dataset, arguing Boas was generally incorrect. The authors go further and state outright that craniometric variation is under more genetic control than environmental plasticity (Sparks and Jantz 2003). Gravlee and colleagues (2003b) address these contradictory findings as differences in research questions, 32 more nuanced understanding of the original study, and general interpretation. Both studies found differences between group cephalic indices suggesting Boas was at least partially correct; the presence of cranial changes between and within groups indicate plasticity. Further, plasticity, by definition, refers to the interaction of genetics and environment to shape bodily form rather than one or the other (Gravlee et al. 2003b). Despite pressures to abandon concepts of biological race, some researchers continued attempts to prove ‘scientifically’ the existence of race, presenting them in evolutionary contexts with discussions of gene flow. The seemingly empirical nature of these types of studies had lasting effects on our discipline and society’s perspective of race, leaving artifacts in current anthropological textbooks and literature still depicting typical cranial forms using 3-group models. The new wave of research aimed to understand the cause of cranial form as it relates to evolutionary processes. These seminal studies did not directly test the relationship of environmental variables on cranial form but rather inferred those relationships based on group differences. Additionally, these studies ignored the facial skeleton which is very likely shaped by plastic responses to the environment. Several studies investigated the complex causations of craniometric variation to determine the evolutionary processes acting upon each cranial measure. Most notably, a series of studies focused on global variation in the Howells craniometric dataset, (Relethford 1994, 2004, 2009; Roseman 2004; Roseman and Weaver 2004; Katz et al. 2016), a dataset comprising 82 measurements for approximately 2,500 individuals from 28 populations (Howells 1973, 1989, 1995). The Howells data have been used for many anthropological studies examining human variation due to the large sample size and the number of populations represented across the world. Relethford (2004) used the Howells dataset to determine the different mechanisms of 33 migration influencing global craniometric variation. Situating his study within the Boas debate on cranial plasticity, he distinguishes between these mechanisms and interprets the role of genetic and environmental influences on cranial form. He argues that as geographic distance increases so too does variation, indicating gene flow is influencing the variation; however, no influence of craniometric variation conceals another. Each study that used Howells data applied different theoretical perspectives to explore the relationship between craniometric and genetic data, shifting emphasis from neutral theory (Relethford 1994) to population structure for assessing the relationship between climate and cranial variation (Roseman 2004; Roseman and Weaver 2004). While greater attention was given to population structures and rates of evolutionary change, their statistical frameworks limited the number of variables contributing to the final models and, thus, limited the analytical power of the methods through this reduced dimensionality. Katz and colleagues (2016) applied a robust quantitative genetic statistical framework, known as Bayesian Sparse Factor Analysis of Genetic Covariance Matrices (BSFG), to explore the relationship between population genetic structure from microsatellites, climate, and cranial measurements from the Howells dataset. Two models were tested: 1) a climate model-with structured random effects and fixed effect for sex, and 2) a reduced model excluding climate effect (i.e., sex effect model). Most of the climate coefficients were negative, including the seven largest magnitude coefficients. These models indicate cranial size and shape are related to climate, particularly the cranial vault, with larger crania associated with colder climates. However, climate coefficients for cranial length measurements are close to zero, indicating less association with climate than cranial width measurements. The authors reported four specific trends: 1) temperature overwhelmed structured random effects at almost all levels of shared 34 history, 2) temperature overwhelmed by structured random effects at almost all levels of shared history, 3) strong population structure effect but weak temperature associations, and 4) weak temperature association and structured effects indicating most variance is residual. Climate is correlated with cranial shape and size with larger body sizes associated with colder temperatures. However, the authors recognize the study is limited to cold temperature variables, likely neglecting many climatic influences in cranial form. This study by Katz and colleagues (2016) is a great model for exploring environmental and genetic influences on cranial form concurrently using robust and appropriate statistical frameworks. However, there remains considerable work to be achieved to truly unravel the evolutionary forces behind human cranial variation, particularly beyond the genome. Literature beyond anthropology, including ecological, quantitative genetics, and clinical literature, should be explored for current methods and theoretical approaches to understanding the relationships between climate, genetics, and facial form given their accessibility to data on living populations and other organisms. Nasal Morphology and Climate Anatomists and geneticists also contributed greatly to the understanding of human facial variation in soft tissue morphology. Most of these studies focus on nasal morphology, an area of the human face directly interacting with the environment as inspired air is modified to protect the health of the respiratory system. Inspired air must be heated and humidified to prevent damage to the lungs, which takes place primarily in the internal nasal passages (Cole 1982). This adaptive response necessitates appropriate changes in the phenotypic expression of the nasal module and, as such, means the nasal complex is highly susceptible to selective forces. Several studies demonstrated a correlation between regional and population differences in internal nasal 35 morphology (e.g., nasal conchae complexity) (Yokley 2009; Noback et al. 2011; Evteev et al. 2014). These studies note populations living in cold-dry climates have higher and narrower nasal cavities than those in hot-humid climates (Noback et al. 2011). Narrow nasal channels coupled with complex nasal turbinates increase surface area between respired air and the mucosa providing ample opportunity for heat and moisture exchange (Maddux et al. 2017). Cold-dry climates create the most respiratory stress of the climatic zones, followed by hot-dry climates; therefore, the greatest variation between groups living between these zones is expected. Hot-wet climates create less respiratory stress making heat reduction the selective pressure in these regions (Maddux et al. 2016); thus, populations in hot-wet regions would have less complex nasal turbinates and wider and lower nasal cavities. This pattern is reflected in Yokley’s (2009) assessment of nasal passages between European and African populations where European populations have adapted to colder climates than African populations. Similar patterning in MMS trait expressions associated with the nasal module would be expected, including nasal aperture width, nasal aperture shape, anterior nasal spine, and nasal bone contour. Researchers demonstrated similar findings in external nasal morphology, as local adaptation plays a role in determining phenotype of the nasal module (Weiner 1954; Maddux et al. 2017; Zaidi et al. 2018). Zaidi and colleagues (2018) identified a significant correlation between nares width and temperature and absolute humidity. The authors report populations living in cold-dry climates exhibit narrower nares than those in hot-humid climates, like the internal nasal morphology above (Zaidi et al. 2017). However, alar base width was not correlated with temperature and only somewhat correlated with absolute humidity, suggesting nares width, and not nose width, may be more associated with adaptation. Maddux and colleagues (2017) tested the relationship of four morphological units of the nose (external pyramid, nasal aperture, 36 internal nasal fossa, and naropharynx) with climatic zones and determined that only the internal nasal fossa was significantly correlated with climatic zones. The authors argue internal nasal structures are more important in respiratory function than the external structures; thus, internal nasal fossa are more influenced by evolutionary and plastic responses leading to population variation. This may indicate that MMS traits, such as Nasal Aperture Width and Nasal Aperture Shape, may not be significantly associated with climate. The dissertation research presented below applies the theoretical frameworks presented above to build upon current knowledge of the relationship between climate, genetics, and human form, specifically in the human craniofacial skeleton. This study aims to recognize the significance of building plastic effects into an evolutionary model for understanding patterns of human craniofacial variation. As presented in the background review, previous knowledge is gathered from human biological and evolutionary studies, as well as anthropology, ecology, and quantitative genetics in attempts to build a comprehensive model for addressing the research goals of this project. 37 MATERIALS AND METHODS This chapter outlines the project design and material and methods used for my dissertation research. The three datasets are described in the Materials section, along with details concerning the collection of these data. Finally, the methods of data analysis are described and outlined by subsections, including data preparation and summary statistics, preliminary analyses, and sections for each of the stated research goals. Materials Three lines of data are used in this study: Macromorphoscopic traits, microsatellite, and climate data. The datasets represent appropriate geographic space to adequately address the research aims. Each line of data was pruned to include only populations and regions with significant representation and geographic proximity across all three lines of data. Macromorphoscopic Traits The highly criticized history of morphological features used to support racial hierarchies in anthropology has deterred many biological anthropologists from using cranial morphology (i.e., nonmetric traits) to study human variation. The quantitative nature of craniometric data was more appealing, which is unfortunate since macroscopic assessments are faster to collect and more useful for fragmentary remains. Traditional methods of ancestry estimation in forensic anthropology viewed these traits as mutually exclusive, considering their presence as distinguishing characteristics for group membership. There was no standardized scoring system or statistical evaluations for group classification. Rather, the practitioner would make ad hoc ancestry assessments based on traits they viewed to be characteristics of a particular population (Klepinger, 2006). This outdated approach used a three-group model to determine the race of unknown individuals, classifying them as European, African, or Asian. Unfortunately, cranial 38 morphology lacked advancement beyond the three-group model for many years. To this end, Hefner (2009, 2018) created the macromorphoscopic (MMS) trait scoring protocol by modifying previously used traits with new standardized scoring systems and the recognition that they are not mutually exclusive to specific populations. This approach provides trait scores on a quasi- continuous scale based on degree of expression of the traits, allowing for variations in trait expression to be captured within and between population, mimicking phenotypic trait expressions. Table 1 provides a list of MMS traits and range of expression scores. Hefner (2009) and others demonstrated systematic patterns of human variation in MMS trait expressions (Hefner and Ousley 2014; Ratliff 2014; Hefner et al. 2015; Klales and Kenyhercz 2015; Plemons et al. 2017). In 2016, Hefner was awarded a National Institute of Justice grant to develop the Macromorphoscopic Databank (MaMD). The goal of this databank was to record global human craniofacial variation in MMS trait scores. 39 Table 1. Macromorphoscopic traits, their abbreviations, and the range of possible score values (character states). Trait Abbreviation Score Anterior Nasal Spine ANS 1-3 Inferior Nasal Aperture INA 1-5 Interorbital Breadth IOB 1-3 Malar Tubercle MT 0-3 Nasal Aperture Shape NAS 1-3 Nasal Aperture Width NAW 1-3 Nasal Bone Contour NBC 0-4 Nasal Bone Shape NBS 1-4 Nasal Overgrowth NO 0-1 Nasofrontal Shape NFS 1-4 Orbital Shape OBS 1-3 Palate Shape PS 1-4 Postbregmatic Depression PBD 0-1 Posterior Zygomatic Tubercle PZT 0-3 Supranasal Suture SPS 0-2 Transverse Palatine Suture TPS 1-4 Zygomaticomaxillary Suture ZS 0-2 The MaMD currently houses trait data for over 8,000 individuals from across the world with 5 levels of population information, including 3-Group ‘Race’, Ancestry, Peer-perceived Ancestry, Geographic Origin, and Population. Each level further refines the geographic region for the individual to improve ancestry estimation and to create more informed models of population variation across geographic space. 3-Group ‘Race’ refers to the traditional broad population assignments (European, African, and Asian); Ancestry further refines the Asian population given the greater geographic region represented in this group (Native American, Eskimo, Asian, Hispanic, and Pacific Island); and Peer-Perceived Ancestry refers to social classifications, such as American White or American Black. Geographic Origin assigns the individual to their country of origin and Population further refines origin for Native American or African populations using tribe information. 40 The detailed geographic information in the MaMD brings MMS data up to current standard procedures available with craniometric analyses and provides the ability to scale population information needed to test their relationship with genetic and climate data. Many biodistance studies rely on patterns of variation in reference materials to situate a study population within a biological and social context, lacking a true understanding of why the variation exists. This project aims to answer the ‘why’ question for MMS trait expression using MMS data obtained from the MaMD. The complete MaMD contains data for 44 populations at the geographic origin level. This study uses a reduced sample comprising 28 geographic groups to best parallel the microsatellite groups (Table 2). Microsatellite Data The microsatellite dataset is compiled by Pemberton and colleagues (2013) from previous global and regional studies. These microsatellite data include 645 loci in 5,795 individuals from the 9 geographic regions presented in Table 2. This dataset is a reduced sample from the original dataset including only populations that overlap in geographic space with individuals in the MaMD. Microsatellite data, or repeated sequences in the genome also referred to as single tandem repeats or STRs, are highly polymorphic, meaning more than one allele can take place at a specific locus, exponentially increasing genetic diversity. This makes microsatellite data highly useful in population studies since related individuals exhibit similar allele frequencies and thus populations that do not exhibit similar allele frequencies are likely genetically unrelated. Conversely, if populations exhibit similar allele frequencies, gene flow likely occurred between those groups. Microsatellite variation is not influenced by natural selection, but rather by gene flow, genetic drift, and mutation (Haasl and Payseur 2013). Therefore, this line of data is useful for defining population differences due to natural selection by controlling for other evolutionary 41 mechanisms. Population structure for these groups has been extensively explored for the Pemberton dataset (Rosenberg et al. 2002; Ramachandran et al. 2005; Rosenberg et al. 2005; Rosenberg et al. 2006; Wang et al. 2007, 2008; Friedlander et al. 2008; Kopelman et al. 2009; Tishkoff et al. 2009; Pemberton et al. 2012; Pemberton et al. 2013). Population information is compared across the MaMD and Pemberton datasets to identify overlap in geographic regions between the two lines of data. Both datasets are reduced to only those groups overlapping in geographic space. The Regional Populations in Table 2 assigned as the metapopulations for all lines of data are defined by Pemberton et al. (2013). These larger scale regional assignments are used for building all statistical models for the project (Table 2). Climate Data Climate information includes coldest month averages for temperature effects and driest months coupled with monthly vapor pressure averages for absolute humidity data. These maximum and minimum values are also used to create ratios of temperature and absolute humidity variation experienced by the regional populations. These data are obtained from NOAA using weather stations in close proximity to the MaMD populations (weather stations provided in Table 2). Average monthly temperature data for 1986 to 2015 are downloaded directly from NOAA reports. This 30-year criteria is the longest timeframe available for modern climate reports and is selected here to more closely reflect temperatures experienced by the skeletal populations. Absolute humidity is calculated using ideal gas law with a constant of 2.1667 gK/J (Vaisala 2013). This calculated as A = C · Pw/T, where A is absolute humidity, C is the constant of 2.1667, Pw is vapor pressure, and T is temperature in Kelvins. Monthly temperature values from NOAA are converted from Celsius to Kelvin. Average monthly vapor pressure data for the same timeframe available for temperature (1986-2015) are obtained from Climate Research Unit 42 (CRU) TS 4.05 (https://crudata.uea.ac.uk/cru/data/hrg/cru_ts_4.05/). Absolute humidity is calculated for each month using the formula above per region. Lowest monthly humidity values are selected for this project to represent driest climatic experience for each regional population. Dry climates have been documented to cause responses in other bony features, such as internal nasal complexities; therefore, if there are correlations between humidity and MMS expression, we can expect to see this through driest month averages (Maddux et al. 2016, 2017). Finally, models test the relationship between the degree of climate variation experienced on trait expressions. To examine the influence of temperature variation on trait expression, the difference is calculated between the number of month averages above 65 degrees and the number of month averages below, modeled from calculations of cooling and heating degree days. Degree days are used for developmental models in ecological studies (Moore and Remais 2014); however, these data are only available for the U.S. datasets in NOAA and could not be applied to the project dataset. However, this model aims to generally reflect climate theories using 65 degrees as a rough sectioning point given this temperature is found as a temperature change that causes physiological response (Moore and Remais 2014). Finally, humidity variation values are calculated by the difference in the range of highest and lowest month averages. These difference values give a rough estimate in variation between the regional groups with greater differences representing greater climatic annual fluctuations in humidity. 43 Table 2. Summary of MMS and genetic populations with sample sizes along with weather station locations. Regional populations refer to the metapopulations assignments that encompass all three lines of data for that region. Regional Weather Populations MMS Local Populations n Genetic Local Populations n Stations Kenya (11), Rwanda (6), Ethiopia (110), Kenya (862), Eastern Africa Somalia (8), Tanzania (8), 38 1986 Garissa, Kenya Rwanda (16), Tanzania (998) Uganda (5) East Asia China (59), Japan (15) 74 China (344), Japan (58) 402 Shanghai, China Jyvaskyla, Europe Finland (165) 165 Orkney Islands (32), Russia (50) 82 Finland Southwest US (101), Colombia (476), Guatemala Santa Marta, Latin America 395 790 Colombia (244), Peru (50) (24), Mexico (250), Peru (40) Colombia Australia (43), Fiji (13), French Polynesia (6), Indonesia (40), Malaysia (5) Australia (20), New Guinea (98), Malacca, Pacific Island New Zealand-Māori (7), 222 198 Polynesia (80) Malaysia New Zealand (27), Papua New Guinea (46), Philippine Islands (35) Udon Thani, SE Asia Thailand (756) 756 Cambodia (22) 22 Thailand Aleut (168), Barrow (8), Siberia 185 Siberia (76) 76 Nome, Alaska St. Lawrence (9) Malawi (19), Cape Town, South Africa 45 Namibia (14), South Africa (76) 114 South Africa (26) South Africa Ghana (82), Ivory Coast (10), West Africa Gambia (3), Nigeria (60) 63 424 Minna, Nigeria Nigeria (332) Methods Data Preparation and Summary Statistics As mentioned above, all three lines of data are pruned for best fit models with significant geographic overlap while also maximizing spatial coverage and skeletal sample sizes. Due to missing data, ten of the 17 MMS traits are used for this study to maintain coverage rather than removing entire populations from the dataset. These traits include ANS, INA, IOB, MT, NAW, NBC, NO PBD, PZT, and ZS. These traits are dichotomized as “0” (absent) or “1” (present). Remaining missing trait data at the individual level is imputed using the ‘mice’ package (Van Buuren and Groothuis-Oudshoorn 2011) in R (R Core Team 2021) to maintain necessary sample sizes for robust analyses. Dichotomized MMS trait data were imputed across the dataset using 44 the predictive mean matching method with 50 iterations to best fit true trait distributions (Figure 1). This imputation method for MMS trait data has been validated by Kenyhercz and colleagues (2019) and sectioning points are provided in Table 3. Table 3. Sectioning points for dichotomizing MMS traits. Trait Sectioning Point Anterior Nasal Spine 1-2 Inferior Nasal Aperture 3-4 Interorbital Breadth 1-2 Malar Tubercle 1-2 Nasal Aperture Width 2-3 Nasal Bone Contour 2-3 Nasal Overgrowth NA Postbregmatic Depression NA Posterior Zygomatic Tubercle 0-1 Zygomaticomaxillary Suture 0-1 Figure 1. Distributions of MMS imputation iterations using mice R package (Blue line is true data; Red lines are imputation iterations). 45 Summary statistics of MMS trait frequencies by regional population are calculated in R and provided in the Results chapter. Additionally, given the MMS data are skeletal data, the population samples span a wide time range with both archaeological and contemporary collections. Therefore, samples are designated as archaeological (“0”) or contemporary (“1”). Contemporary samples include Southeast Asia and Latin America while the remaining samples are archaeological. Preliminary Analyses Genetic distance matrices are generated using delta mu squared distances from population means of microsatellite repeat lengths. The distance matrices are calculated in R using the ‘dist’ function (R Core Team 2021). Matrix A was generated from delta mu squared distances using the ‘distmat’ function from the ‘pracma’ package in R (Borchers 2021) and imported into MATLAB. No further preliminary or summary assessments were performed to explore within data relationships of the genetic data given these data have been extensively examined (Pemberton et al. 2013). However, MMS traits have limited exploration in these types of analyses and, thus, are examined in greater detail. Tests of Mean Measure of Divergence are performed using the ‘AnthropMMD’ package (Santos 2018) in R. This package provides more detailed information for the distance analyses, including the hierarchical clustering plots and significance matrices. These distance matrices provide a means of preliminary exploratory analyses of relationships within and between the genetic and MMS data, before incorporating fixed effects into the model. Finally, a mantel test is performed on the genetic and MMS distance matrices to verify the correlation between these data (Mantel 1967; Cheverud 1988) using the ‘vegan’ package (Oksanen et al. 2020). The mantel test is important for determining if MMS traits can be used as 46 a proxy for genetic relationships like findings for craniometric and that MMS variation has some degree of evolutionary backing. Validating these relationships is an important first step before delving into the environmental forces on craniofacial form and justifying the need to control for population structure in the subsequent statistical models. BSFG Model Setup and Genetic Output Information Samples were collected after a burn in of 500,000 and thinning rate of 1,000 for a total of 1,000 samples for every fixed effect tested. BSFG models were performed using The BSFG model generates some useful genetic information that is worth discussing before delving into the climatic influences of MMS. While heritability and genetic correlations have been validated for the use cranial nonmetric traits and craniometrics as genetic proxies in skeletal studies (Relethford 1994, 2002, 2004; Carson 2006; Martínez-Abadías et al. 2009), MMS has only been briefly explored in one study which suggest good genetic correspondence (Reyes-Centeno and Hefner, 2021). The BSFG provides a bit more detailed information on the relationship between genetics and MMS. The model first determines the number of latent traits, or factors, for the sample. These latent traits represent unobserved traits that underlie variation in the observed traits, in this instance, MMS traits. BSFG assumes these latent traits are present due high modularity in phenotypes from gene networks and developmental pathways. Therefore, the latent traits are used to calculate phenotypic covariances of the MMS traits. Any latent traits that are responsible for less than 1% of phenotypic variation are removed from the sample. Broad sense heritability for each of these are calculated, which are used to determine P, phenotypic covariance, as well as narrow sense heritabilities, which are used to estimate G, genetic covariances (Runcie and 47 Mukherjee 2013). In the Results chapter, each hypothesis will provide the total number of latent traits for the model, narrow and broad sense heritabilities, and matrices P and G. Data Analyses for Climate Effects BSFG models are performed on each of the climate variables. The influences of these variables are measured using the BSFG model described above to accommodate the high dimensionality. The BSFG mixed model (Y=XB+ZU+E) incorporates a series of matrices for the genetic and phenotypic described above. Variable B is a 3*10 matrix which includes an intercept, temporal distance measures, and climate coefficient for the 10 MMS traits for each set of climate analyses. Variable U is a matrix of structured random effects for the sampled populations consisting of the row covariance matrix (H), the genetic covariance matrix (A), and a normal distribution with a mean of 0. A is calculated using the delta-mu squared genetic distance measure, a measure used to calculate interpopulation distances between microsatellite populations for assessing shared population history and time of group divergence (Goldstein et al., 1995a; Goldstein et al., 1995b; Cooper et al., 1999). H uses population variances and covariances assuming the groups had common ancestors but evolved independently. The model combines matrices for disentangled phenotypic covariance (P) and entangled genetic covariance (A). Variable X is a known design matrix for B, as well as Z for U. Lastly, E is an error matrix consisting of column (I) and row (R) matrices and a matrix normal distribution with a mean of 0 (sd=1), where I is the identity matrix for all individuals in the model and R is the residual covariance matrices for the MMS traits. The error matrix illuminates relationships within the data and builds consideration of trait correlations into the model. All individuals from the same population are assigned the same climate and population structure values and all contemporary samples are assigned the same time effect. 48 The BSFG models are performed in MATLAB using an adaptive Gibbs sampler to calculate posterior densities of model parameters (Runice and Mukherjee, 2013) with some modifications provided by Katz and colleagues (2016). Output provides coefficients calculating the relationship between climate and MMS traits with those near 0 having neutral impacts and those further from 0 (negative or positive) being correlated with climate with information for magnitude and directionality of these fixed effects. The posterior means for each climate sample iterations are further explored in R for directionality and importance of climate predictors through additional statistical frameworks, such as linear mixed models, and graphical representation. Contrasts of variance are calculated to compare MMS population relatedness before and after removal of population structure and examine impact of considering the various climate variables on craniofacial variance. 49 RESULTS Results for the summary and preliminary analyses are presented below followed by results for each research goal. Discussions and interpretations of these findings are carried out in the following chapter. Macromorphoscopic Trait Data As discussed, MMS data were dichotomized and then imputed. Figure 2 depicts the distribution of missing data by MMS trait. Hierarchical clustering of these data demonstrates nearly identical group linkages between the imputed and non-imputed data (Figure 3) suggesting the application of this method did not result in significant changes in group relationships. While this method has been validated and preliminary assessments of group relationships demonstrate reliability of imputed data, we should be cognizant of the patterns of missing data while interpreting the results. Specifically, there is a notable high imputation rate for NO. MMS traits are dichotomized following breakpoints provided in Chapter 3. Figure 4 provide a correlation matrix plot for MMS traits. Table 4 provides a summary of MMS trait frequencies for dichotomized data by regional population with “0” being absent and “1” indicating frequencies of present traits. 50 Figure 2. Histogram of missing data distribution for MMS traits. The pattern of missing is depicted in the right plot with gray representing missing and green representing present data. Figure 3. Hierarchical clustering of regional populations for non-imputed and imputed MMS data. 51 Figure 4. Correlation plot of MMS traits. 52 Table 4. MMS trait frequencies for dichotomized scores by Regional Population. n ANS INA IOB MT NAW NBC NO PBD PZT ZS 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 East Africa 38 0.84 0.16 0.92 0.08 0.08 0.92 0.40 0.61 0.55 0.48 0.74 0.26 0.71 0.29 0.45 0.55 0.03 0.97 0.16 0.68 East Asia 74 0.00 1.00 0.97 0.03 0.42 0.58 0.43 0.57 0.91 0.10 0.85 0.15 0.73 0.27 0.89 0.11 0.00 1.00 0.10 0.91 Europe 165 0.32 0.68 0.90 0.10 0.43 0.57 0.52 0.48 0.98 0.02 0.83 0.17 0.76 0.24 0.66 0.35 0.05 0.95 0.07 0.93 Latin America 395 0.13 0.87 0.58 0.42 0.28 0.72 0.49 0.51 0.83 0.18 0.61 0.39 0.24 0.76 0.65 0.35 0.02 0.98 0.26 0.74 Pacific Island 222 0.55 0.46 0.93 0.07 0.19 0.81 0.64 0.36 0.75 0.25 0.82 0.18 0.79 0.21 0.78 0.22 0.06 0.94 0.54 0.46 SE Asia 756 0.37 0.63 0.87 0.13 0.17 0.83 0.51 0.50 0.85 0.15 0.82 0.18 0.87 0.13 0.79 0.21 0.04 0.96 0.44 0.56 Siberia 185 0.41 0.59 0.96 0.04 0.21 0.79 0.83 0.17 0.86 0.14 0.87 0.13 0.82 0.18 0.83 0.17 0.15 0.85 0.62 0.38 South Africa 45 0.73 0.27 0.91 0.09 0.02 0.98 0.44 0.56 0.44 0.56 0.89 0.11 0.80 0.20 0.38 0.62 0.02 0.98 0.58 0.42 West Africa 63 0.71 0.29 0.94 0.06 0.00 1.00 0.32 0.68 0.35 0.65 0.86 0.14 0.73 0.27 0.29 0.71 0.05 0.95 0.56 0.44 53 A population distance matrix for MMS data with significance indicated in the lower triangle for each group distance measure is provided in Table 5. Almost all groups demonstrate significant distances with the exception the three Africa groups. These findings are reasonable given the African groups are in closer geographic proximity than the 6 other regional populations and likely experienced a higher rate of gene flow between them. However, small samples sizes for the African populations may make it difficult to detect significant craniofacial variation for discriminating between these groups. The 2D multidimensional scaling plot with ordinal scaling in Figure 5 depicts group similarities and distance relationships between the regional populations from MMS data. Notice the three regional African populations all plot right of the 0 along the x- axis with East and West Africa plotting in the bottom right quadrant and South Africa in the top right. Latin America is significantly different from all other groups as it is the only group to plot in the bottom left quadrant. East Asia, Europe, Siberia, and Southeast Asia all fall within the top left quadrant; however, all groups display significant distances between them, particularly East Asia and Siberia which show high distance values from all other groups. Table 5. MMS population distance matrix. Upper triangle provides Mahalanobis Mean Measure of Divergence distances while lower triangle indicates significant population distances (*) and non-significance (NS). East East Latin Pacific South West Africa Asia Europe America Island SE Asia Siberia Africa Africa East Africa NA 0.675 0.387 0.483 0.143 0.201 0.332 0.010 0.025 East Asia * NA 0.115 0.333 0.396 0.270 0.457 0.810 0.855 Europe * * NA 0.233 0.244 0.175 0.314 0.556 0.617 Latin America * * * NA 0.387 0.328 0.465 0.615 0.612 Pacific Island * * * * NA 0.029 0.044 0.149 0.237 SE Asia * * * * * NA 0.097 0.229 0.317 Siberia * * * * * * NA 0.303 0.433 South Africa NS * * * * * * NA 0.000 West Africa NS * * * * * * NS NA 54 Figure 5. 2D multidimensional scaling plot of MMS population distances using ordinal scaling. Summary Genetic Data Genetic distance is calculated for the microsatellite dataset using average repeat lengths for the same regional population assignments as for MMS data. The table below provides delta- mu squared distances in the upper triangle and pairwise estimates for the genetic matrix, A, in the lower triangle (Table 6). The three groups from Africa have small genetic distances indicating high group relatedness, mimicking patterns of MMS distances. The association between MMS and genetic distances are depicted below applying a Procrustes transformation to scale the datasets for appropriate comparison in Cartesian space (Figure 6). This plot demonstrates the multivariate relationship between the genetic and MMS trait data, expressed in two dimensional space. Figure 6 demonstrates fair correspondence between the two lines of data. Visually, there appears to be three clusters in the genetic dataset: 1) East, West, and South Africa, 2) Europe, Pacific Island, and Latin America, and 3) Southeast Asia, Siberia, and East Asia. There seems to be only two distinguishable clusters for MMS data: 1) East, West, and 55 South Africa and 2) Pacific Island, Southeast Asia, Siberia, Latin America, and Europe. Unlike the genetic relationships, East Asia plots alone for MMS data showing exaggerated distance from all other groups. A mantel test detects positive and significant correlations between the MMS and genetic population distance matrices from group means (r=0.48, p=0.01), indicating that MMS traits are reasonable proxies for genetic population structure. This also justifies the need to control for genetic effects on MMS trait expression when examining the traits for climatic influences. Table 6. Genetic distance matrix for regional populations. Upper triangle provides delta- mu squared Euclidean distances from group means. Lower triangle provides pairwise estimates for the genetic matrix, A. East East Latin Pacific South West Africa Asia Europe America Island SE Asia Siberia Africa Africa East Africa 0.000 7.476 6.743 10.029 8.354 9.235 8.079 2.412 1.631 East Asia 0.465 0.000 5.474 4.902 5.244 2.659 1.789 11.340 11.146 Europe 0.518 0.608 0.000 6.064 7.972 7.264 4.446 10.913 10.627 Latin America 0.282 0.649 0.566 0.000 9.507 7.357 4.581 13.977 13.737 Pacific Island 0.402 0.625 0.430 0.320 0.000 6.716 6.008 12.415 11.674 SE Asia 0.339 0.810 0.480 0.474 0.520 0.000 4.635 13.279 12.873 Siberia 0.422 0.872 0.682 0.672 0.570 0.668 0.000 12.144 11.973 South Africa 0.827 0.189 0.219 0.000 0.112 0.050 0.131 0.000 2.101 West Africa 0.883 0.203 0.240 0.017 0.165 0.079 0.143 0.850 0.000 56 Figure 6. Procrustes plot of genetic and MMS distance matrices demonstrating group relatedness within and between the datasets. Climate data are obtained from NOAA (https://www.ncdc.noaa.gov/cdo-web/datasets) and CRU TS 4.05 (Harris et al. 2020). Table 7 provides summary climate data for each variable used in this study by regional population. The average coldest month temperatures are used to assess temperature effects on MMS expression (column 1, Table 7). The difference between the number of month averages above and below 65 degrees Fahrenheit are provided under Temperature Ratio. For example, all month averages for East Africa are above 65 degrees resulting in a difference of 12 (n above 65°– n below 65°). This means that higher values represent very little annual climatic fluctuations while ratios of 0 have high climatic variability. As previously stated, absolute humidity values were calculated using monthly temperature and vapor pressure data and ideal gas law with a constant of 2.1667. The lowest month absolute humidity values were selected for each region to test the influence of dry climates on MMS trait 57 expression, resulting in the values provided in Table 7. The final column of this table provides the difference of the range of monthly absolute humidity values to test climatic variation in humidity on craniofacial form. Here, higher values represent greater annual climatic fluctuations while lower values represent little variation. Regional climate values are assigned to every individual within the population and imported into MATLAB to calculate climate coefficients for MMS traits. Climate coefficients for all explored variables are provided below in hypothesis-specific contexts. Table 7. Summary climate data used for BSFG analyses. Regional Population Cold, Celsius Absolute Humidity Temperature Ratio Humidity Range East Africa 26.6 13.63 12 0.70 South Africa 12.4 6.10 2 6.57 West Africa 25.3 9.13 12 10.20 East Asia 6 2.26 0 9.13 SE Asia 22.6 15.54 12 5.92 Europe -8.1 2.46 12 7.41 Pacific Island 26.8 20.33 12 1.16 Latin America 16.2 17.99 4 1.58 Siberia -13.6 0.97 12 8.61 Genetic Effects on MMS Variation The genetic data summarized here reflects posterior means and exploration of data obtained from the BSFG model for temperature data. These data vary slightly between each BSFG model given the values are averages of the 1000 sample iterations performed, but these values are a good representation of the overall genetic information. The BSFG model calculates 18 factors for the humidity model, while the humidity BSFG calculates 17. Narrow sense heritability, representing additive genetic effects, ranges from 0.06 (low heritability) to 0.51 (moderate heritable). Broad sense heritabilities, incorporating genetic and phenotypic variation, 58 range from approximately 0 (no heritability) to 0.95 (high heritability). However, most of the broad sense heritabilities for the factors have low to moderate measures. Narrow sense heritability values for the factors are used to estimate the genetic covariances (G) (Table 8) and broad sense heritabilities estimate phenotypic covariance (P) (Table 9). Notably, IOB negatively covaries with almost all other traits in both matrices except for MT. This means that as IOB becomes more greatly expressed (i.e., interorbital breadth is wider), so too does MT (i.e., tubercle is larger), but all other traits are less expressed. Following IOB, NAW negatively covaries with INA, MT, and PBD. Interestingly, NAW and IOB are both traits capturing facial breadth. A heatmap of trait covariances incorporating temperature and genetic effects provides a graphical representation of matrices G and P (Figure 7). Patterns of small covariances between ANS, INA, and NAW suggest some level of modularity in the nasal complex. Table 8. Genetic variance-covariance matrix (G) among traits calculated from narrow sense heritability estimates of factors. ANS INA IOB MT NAW NBC NO PBD PZT ZS ANS 0.01107 0.00132 -0.00044 -0.00013 0.00100 0.00079 0.00096 0.00073 0.00020 0.00090 INA 0.00132 0.01330 -0.00010 0.00050 -0.00027 0.00017 0.00078 0.00062 0.00001 0.00069 IOB -0.00044 -0.00010 0.01714 0.00066 -0.00135 -0.00180 -0.00053 -0.00021 -0.00033 -0.00078 MT -0.00013 0.00050 0.00066 0.01010 -0.00045 -0.00032 0.00018 0.00022 0.00020 0.00008 NAW 0.00100 -0.00027 -0.00135 -0.00045 0.01030 0.00120 0.00022 -0.00046 0.00003 0.00082 NBC 0.00079 0.00017 -0.00180 -0.00032 0.00120 0.01695 0.00033 0.00011 0.00023 0.00010 NO 0.00096 0.00078 -0.00053 0.00018 0.00022 0.00033 0.01625 0.00009 0.00023 0.00048 PBD 0.00073 0.00062 -0.00021 0.00022 -0.00046 0.00011 0.00009 0.02009 0.00035 0.00056 PZT 0.00020 0.00001 -0.00033 0.00020 0.00003 0.00023 0.00023 0.00035 0.00414 0.00067 ZS 0.00090 0.00069 -0.00078 0.00008 0.00082 0.00010 0.00048 0.00056 0.00067 0.00985 59 Table 9. Phenotypic (MMS) variance-covariance matrix (P) calculated among traits from broad sense heritability estimates of factors. ANS INA IOB MT NAW NBC NO PBD PZT ZS ANS 0.20838 0.01777 -0.00527 -0.00148 0.01857 0.01129 0.01185 0.00777 0.00411 0.01337 INA 0.01777 0.13360 -0.00183 0.00745 -0.00522 0.00608 0.01128 0.00935 -0.00084 0.00870 IOB -0.00527 -0.00183 0.17917 0.01056 -0.01971 -0.03204 -0.00954 -0.00383 -0.00441 -0.01175 MT -0.00148 0.00745 0.01056 0.24728 -0.00691 -0.00465 0.00324 0.00344 0.00242 -0.00002 NAW 0.01857 -0.00522 -0.01971 -0.00691 0.14497 0.01804 0.00195 -0.00947 0.00191 0.01424 NBC 0.01129 0.00608 -0.03204 -0.00465 0.01804 0.17816 0.00662 0.00268 0.00300 0.00136 NO 0.01185 0.01128 -0.00954 0.00324 0.00195 0.00662 0.17001 0.00309 0.00335 0.00568 PBD 0.00777 0.00935 -0.00383 0.00344 -0.00947 0.00268 0.00309 0.20493 0.00452 0.00672 PZT 0.00411 -0.00084 -0.00441 0.00242 0.00191 0.00300 0.00335 0.00452 0.04715 0.00994 ZS 0.01337 0.00870 -0.01175 -0.00002 0.01424 0.00136 0.00568 0.00672 0.00994 0.22330 Figure 7. Heatmap of traitwise covariance matrix for MMS traits with inclusion of both temperature and genetic effects. Yellow indicates high covariance and blue indicates low. Broad sense heritabilities (h2) of MMS traits are obtained from matrix P (Table 10). The MMS h2 values range from 0.05 to 0.24. Traits MT, ZS, ANS, and PBD have the highest h2 measures with over 20% of variance being due to genetics. Unfortunately, the BSFG does not 60 provide population specific values. Therefore, these h2 values are aggregated measures of population heritabilities and should be evaluated for population-specific heritabilities in the future. Despite these low h2 values, a mantel test for the correlation between matrices G and P indicates good correspondence for these data (r = 0.62, p = 0.00). This suggests that a more in- depth analysis of h2 is needed. Notably, the traits with lower heritability measures are all part of the nasal complex (INA, NAW, NBC, and NO), as well as IOB, a facial breadth trait. These traits capturing the nasal complex and facial breadth have similar heritabilities ranging between 0.14 and 0.16. ANS is the only nasal trait not included in this trend. Table 10. Heritability measures of MMS traits. Trait h2 ANS 0.22 INA 0.14 IOB 0.16 MT 0.24 NAW 0.15 NBC 0.17 NO 0.16 PBD 0.21 PZT 0.05 ZS 0.23 Temperature Effects on Craniofacial Morphology Research aim 1 explores coldest month averages to determine the correlation of temperature with MMS trait expression. Figure 8 is a forest plot of climate coefficients for coldest month averages with 95% confidence intervals. Coefficients are pulled directly from the β matrix in MATLAB output files and confidence intervals are calculated from the 1000 sample iterations for fixed effects for each MMS trait. The plot is sorted by coefficient magnitudes to show the relationships between MMS expression and their correlation to cold temperatures. 61 Those furthest from 0 demonstrate climatic influence on the traits and positive and negative values represent the directionality of the correlation. All traits exhibit positive correlations with cold climate except for NAW, PBD, and IOB, but none of these are reliably nonzero. ZS and NBC have the greatest positive correlations with ZS being reliably nonzero. Despite ZS having one of the highest heritability measures, this trait exhibits the greatest positive temperature effect. This may indicate that cold climates overpower the genetic effects on ZS expression with increased sutural complexity in colder regions. The magnitudes of these temperature effects are small and only one trait is reliably nonzero. However, noting preliminary patterns in these data is important for more refined analyses in the future. Figure 8. Climate coefficients (black dots) with 95% confidence intervals for coldest month averages. Values right of the red line have positive directionality, and those to the left have negative directionality. 62 A comprehensive model is used to assess contrast of variance in craniofacial morphology between groups incorporating temperature difference and phenotypic and residual error variance, while removing effects of population structure. Figure 9 depicts linear relationships between pairwise temperature differences between regional populations and contrast of variances. A Procrustes plot demonstrates the difference between group relatedness between MMS Mean Measure of Divergence distance and the contrast of variance model (Figure 10). The plot shows that when genetic variance in removed for the temperature effects model many groups cluster very tightly indicating that temperature alone does account from group dissimilarities. Figure 9. Comparison of temperature effects on contrasts of variance for ANS (left) and INA (right). The X-axis plots temperature difference for pairwise regional populations and Y-axis in contrast of variance. ANS is the only trait demonstrating a strong correlation with temperature when population structure is removed. The remainder of the traits resembles INA with weak correlation. 63 Figure 10. Procrustes plot of group relatedness for MMS Mean Measure of Divergence Distances (red) and contrasts of variance (black) for temperature effects, phenotypic and residual error variance, and population structure. Humidity Effects on Craniofacial Morphology Research aim 2 explores the effects of humidity on MMS trait expression using driest month averages (Table 7). Figure 11 depicts humidity coefficients for driest month averages with 95% confidence intervals. Like temperature effects, coefficients are obtained from the β matrix and confidence intervals are calculated from the 1000 sample iterations. The plot shows the magnitude and directionality of humidity effects on each MMS trait. All traits exhibit positive correlations with humid climate except for MT, NAW, NBC, and PZT with MT directly at 0. IOB and NO have the greatest positive correlations with IOB being the only trait reliably nonzero. Interestingly, IOB has one of the higher heritabilities measures, yet this suggest strong humidity effects on the trait variance. This may be a result of muddled information from dichotomizing these data. INA and ANS are traits present at the 64 inferior portion of the nasal aperture which seem to correlate (Figure 11), and both can be seem with some degree of positive influence by humidity. This means that these traits have greater expression in wetter climates. Conversely, NAW and NBC which capture shape of the nose (width and taller bridges of the nose) are negatively correlated meaning these traits are less expressed in drier climates. However, these are weak magnitudes which may be overpowered by genetics. Figure 11. Climate coefficients (black dots) with 95% confidence intervals for driest month averages. Values right of the red line have positive directionality, and those to the left have negative directionality. Again, contrast of variance in craniofacial morphology is assessed between groups incorporating humidity difference and phenotypic and residual error variance, removing population structure. The Procrustes plot (Figure 12) demonstrates the difference between group relatedness between MMS Mean Measure of Divergence distance and the contrast of variance 65 model for humidity. Here, the plot depicts how removing genetics effects from the humidity model influence population relatedness via craniofacial phenotypic variance comparing contrast of variance to distances from trait expression means. There are clear group distances between the regional population when incorporating humidity effects while removing population structure and, thus, absolute humidity may play a key role craniofacial variation due to selection pressures. There appears to be greater divergence between the groups using this comprehensive approach to measuring phenotypic variance, more so than the temperature effect model (Figure 12) Figure 12. Procrustes plot of group relatedness for MMS Mean Measure of Divergence Distances (red) and contrasts of variance (black) for humidity effects, phenotypic and residual error variance, and population structure. Effects of Climatic Fluctuations on Craniofacial Morphology The final aim of this research is to examine if greater annual climate variation shapes craniofacial morphology. This is tested using two climate variables, annual ratio of months above and below 18 degrees Celsius (65°F) and range of annual absolute humidity. The range of annual temperature (Figure 13) and humidity (14) are depicted below to show within and 66 between climate variation for each region. Siberia, Europe, and East Asia have the widest distribution of monthly temperature ranges, yet Siberia and Europe fall completely below the 18- degree breaking point. East Asia, Latin America, and South Africa fall almost equally on either side of the breaking point making their variance values the highest. The remainder of the groups East and West Africa, Southeast Asia, and Latin America all fall completely above the break line indicating these regions experience high temperatures year-round (Figure 13). Regarding absolute humidity, the regions demonstrate greater annual variability than seen for temperature. East Asia, Europe, Siberia, and South Africa experience dry to moderate humidity throughout the year, while Southeast Asia experiences moderate to high humidity. East Africa, Latin America, and Pacific Island have narrow annual humidity distributions with East Asia having moderately high humidity and Pacific Island and Latin America having high year-round humidity. Figure 13. Plot of annual monthly temperature averages for (1986-2015) for each regional population. The x-axis is temperature in degrees Celsius from coldest to warmest. The red line indicates 18-degree breaking point for ratio calculations. 67 Figure 14. Plot of annual humidity temperature averages for (1986-2015) for each regional population. The x-axis is absolute humidity values from driest to wettest. The BSFG model calculated 14 factors for estimating matrices G, P, and R. Figure 15 depicts coefficients for temperature ratios with 95% confidence intervals. MT, PZT, IOB, NBC, and ZS show positive directionality with greater climate variance, while INA, ANS, and PBD exhibit negative correlations. NAW and NO are 0 with neither negative nor positive correlations; however, none of the traits are reliably nonzero. Again, INA and ANS group together as in the previous two models suggesting modularity of the nasal complex features. NAW and NO are also close in magnitude and directionality for this effect model. While modularity is suggested in response to climatic fluctuations, the dichotomized structure of MMS data diminishes the ability to illuminate whether climate variability leads to more intermediate expression of the traits. This may be why coefficients for this model are small and no traits are reliably nonzero. The effects of climate variance is also examined with annual humidity ranges. Figure 16 plots magnitude and directionality coefficients for humidity ranges with 95% confidence 68 intervals. INA, ZS, PZT, NBC, NAW, and MT exhibit positive directionality demonstrating greater trait expression with greater annual variability in humidity. INA is reliably positive indicating a nasal sill is more prominent under these environmental conditions. INA also exhibits positive correlation in previous humidity analyses with greater nasal sills in more wet environments. Interestingly, while IOB is reliably nonzero in previous humidity analyses (wide breadth in wetter environments, this trait has slight negative directionality when humidity is more variable. This shift in IOB may be indicative of trends in more intermediate traits as IOB becomes smaller with greater degrees in annual humidity fluctuations. Directionality of temperature ratio and humidity range coefficients demonstrate similar patterns between the two models. The most notable differences are INA and IOB where IOB is positively correlated with temperature variance and slightly negative with humidity variance. This pattern is reversed for INA where INA is negatively correlated with temperature variance and reliably positively correlated with humidity variance. 69 Figure 15. Climate coefficients (black dots) with 95% confidence intervals for annual temperature ratios. Values right of the red line have positive directionality, and those to the left have negative directionality. 70 Figure 16. Climate coefficients (black dots) with 95% confidence intervals for annual ranges of absolute humidity. Values right of the red line have positive directionality, and those to the left have negative directionality. Interaction of Genetics and Climate on Craniofacial Morphology Phenotypic variance is due to environmental and genetic factors, but the magnitude of these factors varies depending on the populations and variables under investigation. Therefore, it is important to establish these relationships when estimating biodistance in both archaeological and forensic contexts. Without rooting our methods deep in evolutionary theory, we increase the likelihood of misinterpreting our population histories or diminishing our ability to identify unknown peoples. The following discussion synthesizes the magnitude of genetic and climatic effects on MMS trait expression. Genetic Effects Gross level assessments of these data were conducted as precursors to explorations of the various climate variables. There were presumptions that MMS traits have good correspondence 71 with genetic data given previous preliminary research (Reyes-Centeno and Hefner 2021) and craniometric (Relethford 1994, 2002; González-José et al. 2004; Roseman 2004) and cranial nonmetric findings (Carson 2006; Martínez-Abadías 2009). Therefore, the preliminary nature of these analyses could not verify neutral evolution influencing MMS trait manifestations. Mantel tests for population distances from group means and for comparison between matrices P and G from the BSFG results indicate significant correspondence between MMS and genetic data. The Procrustes plot (Figure 6) and the matrices (Tables 5 and 6) for these data show significant clustering of East, West, and South Africa within and between the two datasets. This reflects the close geographic proximity between these three groups resulting in greater genetic relatedness than other regional groups in this study. All other groups seem to correspond well in cartesian space except for the distances between Siberia and East Asia. These two regions are closely related genetically with Siberia plotting between East Asia and Southeast Asia, but Siberia plots between Europe, Latin America, and Southeast Asia for MMS distances. East Asia appears to be phenotypically distinct from all other groups. This may be a product of small sample size in the East Asian MMS data. Heritability estimates of factors and MMS traits demonstrate low to moderate heritabilities for these craniofacial morphological features. However, given the correlation between genetic and phenotypic matrices, these heritabilities are likely not reliable measures. Heritability is measured from genetic and phenotypic variance. When we dichotomize MMS traits, we are greatly reducing the variance amongst each trait (i.e., phenotypic variance values). The use of group means to calculate matrices P and G likely reduces model accuracy compared to complete datasets which would more accurately calculate these variance-covariance matrices 72 for heritability estimates. Despite these shortcomings, the genetic effects results serve as a rough estimation to explore these data. The most notable finding from these initial phases of genetic effects is that traits associated with the nasal complex (INA, NAW, NBC, and NO) and facial breadth (IOB) have the lowest heritability estimates. This suggests these features likely respond to selection pressure more than neutral evolutionary forces (migration, genetic drift, or mutation). These traits are of particular interest throughout the assessment of climate effects in the series of research aims. Conversely, this could indicate features with higher heritabilities should be of greater concern in modeling biodistance and estimating ancestry, but these would need to be explored in much greater detail for definitive claims. Given these values are more relational rather than accurate measures in this early stage of research, selective pressures may still be overshadowed by genetic effects in terms of their evolutionary significance and, thus, continue to serve as proxies for genetic relatedness. Current genetic estimates should be considered as the minimum influence of population structure on phenotypic variance due to the inherent limitations of these data (e.g., temporal inconsistencies, unmatched datasets, dichotomized traits, and group means) (Roseman 2004). Climate Effects Climate effects are discussed by their estimated heritabilities to grasp some degree of the ratio of genetic and climate effects on each MMS trait. First, we review traits with higher heritabilities, including ANS, ZS, MT, and PBD. Theoretically, we should expect these more heritable traits to exhibit lower correlations with climate variables due to reduced selective pressure with higher neutral genetic effects. 73 High Heritability Traits ANS ANS is the only feature located within the nasal complex with a higher heritability measure, yet ANS positively covaries with all other nasal complex traits. ANS exhibits zero correlation with cold climates but does present a weak positive correlation with humidity. This may indicate that higher noses (greater projection from the face) are more prominent in colder climates, consistent with previous findings of nasal morphology and structure in cold-dry climate (Noback et al. 2011). Greater nasal projection increases surface-volume ratio of the internal nose, thereby increasing the ability to exchange moisture and warmth with the mucosa (Franciscus and Long 1991; Cole 2000; Churchill et al. 2004; Inthavong et al. 2007, 2019). However, when climate is more variable, ANS decreases in size (Figures 15 and 16). The nose becomes less pronounced, decreasing surface-volume ratio, and improving the ability of the nasal module to respond to changes in the environment. Cold-dry and hot-wet environments result in the greatest thermoregulatory stress whereby cold-dry environments (Noback et al. 2011; Maddux et al. 2017) cause significant respiratory stress and hot-wet environments lead to greater risk of hyperthermia (Frye & Kamon 1983; Havenith 2005; Taylor 2006a, b; Maddux et al. 2017). The nasal complex must respond to increased respiratory stress in cold-dry environment to prevent damage to the lungs (Noback et al. 2011; Maddux et al. 2017). Therefore, when climates have greater annual fluctuations, the respiratory system undergoes less continual stress reducing the need for the nasal complex to respond to the environmental pressures. While temperature effects of cold weather climate we not detected, humidity captures some degree of temperature effects given absolute humidity is directly correlated with 74 temperature. Temperature alone may not be enough environmental stress to necessitate a plastic response during development, but humidity and temperature combined could result in these measurable changes. The genetic influence of ANS may be significant enough to mask the temperature effects alone. MT MT has the highest heritability estimate of all MMS traits. Climate coefficients for temperature, absolute humidity, and annual range of absolute humidity were all 0. The only climate correlation with MT was positive directionality when tested for annual temperature variability (i.e., larger MT in more variable temperatures). Given this correlation is not reliably nonzero and all other climate coefficients has no directionality, we can conclude that MT expression is likely a reflection of population structure. The MT serves as a muscle attachment site for the masseter muscles used for mastication (Hefner and Linde 2018). Therefore, the MT serves no functional purpose directly related to climate. This trait is located just lateral to ZS, both of which are incorporated into the zygomaticomaxillary complex and demonstrate the highest heritabilities for MMS traits. It appears that trait expression within the zygomaticomaxillary complex are due to neutral evolution. ZS ZS is interesting because not only does this trait exhibit the second highest heritability, but it is also the only trait reliably positively correlated with cold climates. ZS is also mildly positively correlated with humidity variation yet lacks any degree of correlation between temperature ratio values of absolute humidity. These patterns indicate that greater sutural complexity corresponds with colder climates and greater annual fluctuations in absolute 75 humidity. However, sutural complexity is not affected by humidity or annual temperature variability. Cranial sutures function to hold together various cranial elements and absorbs and transmits mechanical stress. Sutures with greater complexity are capable of withstanding and distributing sutural dynamic stress from surrounding bones (Zhang and Yang 2015). The zygomaticomaxillary suture connects the maxilla and zygomatic bones of the midfacial skeleton and help to manage mechanical stress from the adjacent masseter muscle (Hefner and Linde 2018). Therefore, greater expression of ZS should be suited to deal with greater masseteric stress. The form-function relationship of ZS appears to have little direct association with temperature or humidity, supporting greater genetic effects due to neutral evolutionary processes rather than climate effects. The correspondence between ZS expression and temperature and humidity variability may be due to patterns of subsistence across the regions, as well as temporal differences between the skeletal samples that may represent great contrasts in subsistence patterns and available resources based on the populations local environments at that time. Katz and colleagues (2017) found that the transition to softer diets with transition to agriculture did change global patterns of cranial size and shape; however, diet effects were not significant enough to overwhelm genetic forces of evolution. Expression of ZS is likely more strongly controlled by neutral evolutionary processes rather than temperature and humidity, but other spatially distributed environmental factors should not be excluded. PBD High heritability of PBD suggests that the expression of this trait is also controlled by neutral evolutionary processes. PBD does exhibit negative directionality in all three phases of analysis suggesting PBD is present more often in hot-wet environments. These results are 76 misleading, because only the three groups from Africa express a PBD. As discussed earlier, these groups are highly related genetically and phenotypically due to close geographic proximity. If these climatic patterns were true, we would expect to also see the trait present in groups from other hot-wet environments who are not genetically closely related. We can reasonably conclude that PBD is due to genetic forces with little to no climatic influence. The functionality and origins of PBD are still largely unknown. This trait appears at location of the anterior fontanelle present at birth. The most accepted causation of this trait to date suggests that PBD occurs when intrinsic and extrinsic stressors occurring during ontogeny resulting in reduced ability to ossify the fontanelle (Wood 2012; Hefner and Linde 2018). More research is needed to understand the functionality of this trait before we can conclude on the evolutionary forces. For example, PBD may be more present in hot climates to help release heat from the brain and reduce the risk hyperthermia. Studies find other cranial features that help to cool the brain for reducing hyperthermia (Mariak et al. 1999; Jessen 2001; White et al. 2011), so it seems unreasonable to rule out this functionality given so little is known about PBD. Low Heritability Traits The remainder of traits presented in this section are those with relatively low heritabilities. First, we review four features incorporated into the nasal complex, INA, NAW, NBC, and NO, followed by IOB and PZT. The nasal complex is of particular interest given the wealth of research demonstrating the significance of the nose to regulating inspired air in humans and how this shapes our internal and external nasal morphology (Weiner 1954; Farkas et al. 1986; Franciscus and Long 1991; Churchill et al. 2004; Ferrario et al. 2007; Inthavong et al. 2007; Bigham et al. 2009; Yokley 2009; Keck and Lindemann 2010; Noback et al. 2011; Guo et al. 2014; Shaffer et al. 2016; Zaidi et al. 2018) 77 INA Unfortunately, when reviewing these data, it became apparent that all regional populations have high frequencies of “0” for INA with 87% or greater. Latin America is the only exception with only 58% absence. Therefore, these data are largely uninformative and the sectioning point for dichotomizing INA need to be reevaluated. INA makes up the nasal anterior aspect of the nasal aperture floor. This internal aspect of this region of the nose is lined with mucosa to aid in regulating inspired air. INA represents part of the nasal complex and one of the first points of contact with inspired air. The anatomy and functionality of INA coupled with the fact that other traits associated with the nasal module exhibit some level of climatic influence, it is likely that this trait responds to selective processes regardless of population structure. This need to be reexamined with modified sectioning points or the full suite of trait expressions. NAW NAW was predicted to have strong climatic effects due to influence of nasal width on overall surface-volume ratio of the nose. Despite the limited evidence of genetics on this trait and previous research correlating nasal breadth and climate, all climate coefficients had no significant directionality. Interestingly, Roseman (2004) and Katz et al. (2016) did not find nasal breadth measures to be significantly correlated with cold climate, yet nasal height was significant. The correspondence between nasal height in previous studies and nasal projection via ANS may indicate that these two features work together to increase nasal surface-volume ratio rather more so than nasal width. It should be noted however that NAW can be subjective as expression is scored relative to overall facial breadth. Wide facial breadths may result in observers to underscore which skew frequencies toward “0” scores across many regional 78 populations (Table 4) and vice versa for narrow crania. Roseman (2004) and Katz et al. (2016) found wide cranial breadth measures to be associated with cold climates. Consequently, these high associations in cranial breadth may cause macroscopic assessments of relative nasal widths to appear more intermediate, information that gets lost when dichotomizing the traits. More accurate assessments of NAW will require refinement of these data with raw trait scores. NBC NBC exhibits positive directionality indicating lower nasal bone curvature with cold climates. No other climate variables determined directionality for NBC. However, again, I believe this to be a result of dichotomized data losing a great deal of variation. Frequencies of NBC by regional population (Table 4) show similar frequencies of NBC across all groups demonstrating absence of the trait (i.e., low rounded versus high nasal bridge). Hefner and Linde (2018) provide frequencies by group for raw trait data which show great diversity between these regions. Raw frequency trait data shows higher rates of expression for high nasal bridges in colder populations, like European and East Asian (Siberian not specified), corroborating studies on the association of greater projection of the nose for cold-dry environments. As with INA and NAW, this trait needs to be reexamined with full consideration of variation through raw trait scores. NO NO is positively correlated with both temperature and humidity meaning NO is more often expressed in hot-wet climates and absent in cold-dry environments. There is no indication that NO expression responds to climatic fluctuations for neither temperature nor humidity. NO is scored as either present or absent, so there is no concerns about the loss of variation with dichotomizing this trait. 79 NO occurs from ossification of cartilage anchored to the nasal bones. This ossification can occur during growth and development or during adulthood (Hefner and Linde 2018). The lower frequencies of this trait in cold-dry environments may indicate that ossification around the nose is reduced during growth and development to improve the surface-volume ratio of the nose. IOB IOB negatively covaries with almost all other traits to suggest that as other features are less expressed, IOB increases in expression. Only MT covaries positively with IOB which may be due to some degree of modularity in the midfacial skeleton. IOB demonstrates negative directionality with temperature meaning that IOB is narrower in cold populations. Examination of temperature ratio values estimate positive directionality meaning that IOB exhibits greater expression in more variable climates. These patterns are consistent with studies mentioned throughout stating narrow and tall nasal passages are associated with cold-dry environments. Maddux et al. (2017). Measurements of interorbital breadth (dacryon-dacryon) are also used to capture upper fossa morphology of the internal nasal fossa (Franciscus, 1995, 2003; Yokley 2009; Noback et al. 2011; Maddux et al. 2017). This region of the nose includes nasal turbinates which provide the greatest mucosal space for regulating inspired air. Narrower passages increase force of inspired air through the internal nasal fossa improving the ability to warm and moisten the air. Likewise, IOB exhibits reliably nonzero positive directionality for absolute humidity with wider narrow breadths associated with dry climates. Evidence of wider IOBs with temperature variability suggests the potential for more intermediate traits with greater annual fluctuations due to reduced respiratory pressure from continuous cold climates. This is consistent with the 80 findings for ANS suggesting the nose is tall and narrow in cold-dry climates with reduction of the features in more variable environments. PZT Finally, PZT has the lowest heritability estimates of all MMS traits. The trait exhibits positive directionality for temperature and temperature and humidity annual variability measures, while showing negative directionality for humidity. There appears to be greater magnitude for the temperature variables than for humidity. These patterns indicate that PZT is less expressed in cold environments. The PZT anchors the temporalis muscle which is important for mastication. Therefore, expression of this trait may be a response of more indirect climate forces or other environmental factors, like the discussion for MT. Frequencies for trait presence are high for all regional groups except for Siberia; thus, this group is driving the climatic effects for PZT. 81 DISCUSSION AND CONCLUSION This chapter comprehensively discusses the findings of the analyses described above to present the overall interactions of genetics and climate on craniofacial morphology. I contextualize these findings to similar research, discuss the broader evolutionary significance of the findings, and present limitations of this project with suggested future avenues of research. While this study aimed to measure the influence climate and genetics on MMS traits, this was the first time exploring these data through an evolutionary lens, and the finding are simply preliminary assessments of general global patterns of variance and relationships between the various lines of data. Therefore, the findings raise more questions and open more avenues of research more so than answering definitively the questions sought to answer. Nevertheless, patterns can be discussed to guide future research for delving deeper into these evolutionary questions on phenotypic response to our genes and our environments. Contextualizing Causative Forces of MMS Trait Expression To contextualize the findings above, we situate these results within previous research starting with preliminary works by Reyes-Centeno and Hefner (2021) on the genetic relationships with MMS traits. We then briefly explore research on the heritability of nonmetric traits followed similar mixed models examining the genetic and climate effects on cranial form using cranial measurements. MMS Genetic Relatedness The study on neutral evolutionary processes on MMS traits explored seven populations (African American, White American, Chinese, Japanese, Columbian, Mexican, and Papuan) estimated similar correlation between MMS and genetic distances compared to those provided above (Reyes-Centeno and Hefner 2021). Reyes-Centeno and Hefner performed a mantel test 82 between distance matrices calculating significant correspondence (r = 0.64, p = 0.009) which reflects the mantel performed on matrices P and G in this studying (r = 0.62, p = 0.00). These two studies validate that MMS traits serve as reliable proxies for population structure and history and, thus, are reliable for biodistance analyses which aim to measure genetic relatedness from phenotypic variation. The authors of this research also noted traits they believe to respond to genetic drift or natural selection. Those noted to evolve greatest under neutrality include ZS, PZT, and ANS (Reyes-Centeno and Hefner 2021). Heritability estimations for ZS and ANS in this study corroborate these conclusions. PZT did not correspond to previous findings, but these data could be potentially problematic for this study and need to be explored in greater detail. PBD, a trait found likely to evolve under neutrality in this study, plots more toward selection for Reyes- Centeno and Hefner; however, trait utility for PBD in their analyses hover around 0. Those noted as most likely being due to selective processes include MT and NAW (Reyes-Centeno and Hefner 2021). MT was found to have higher heritabilities in this study, but there also appears to be some degree of evidence for climatic correlations. Therefore, evolutionary forces of this trait are not fully disentangled. Evolutionary processes of NAW are still believed to be due to selection as suggested by both studies, but NAW needs to be further assessed using raw trait data to capture true phenotypic variance. IOB and NO, two traits found to have stronger climate effects here, also tend toward selection in Reyes-Centeno and Hefner (2021). Overall, there seems to good correspondence between predictions of evolutionary mechanisms for many of traits between the two studies. It should be noted that the studies used the same overall genetic (microsatellite) and MMS datasets, but we used different population 83 subsets which will create some disagreement. Thus far, it seems that we can reasonably conclude that ZS and ANS evolve primarily due to neutral evolutionary processes while IOB evolve largely through selection processes. Cranial Metric Mixed Models for Genetics and Climate This research largely builds on a series of works testing the proportion of neutral evolutionary versus selection processes on human crania using regional craniometric data (Roseman 2004; Katz et al. 2016). The results of this research need to be situated within these previous studies to illuminate patterns of concordance or discordance of evolutionary forces on the facial skeleton. Many patterns were mentioned when interpreting results in the previous section of this chapter, Interaction of Genetics and Climate on Craniofacial Morphology, and are only briefly summarized here for comprehension. Both previous studies found cold climate to be significantly associated with many cranial measurements illustrating that cranial form does not respond solely to genetic factors. Roseman (2004) notes similar correlation via Mantel’s test between genetic and phenotypic matrices (r = 0.495, p = 0.006) to MMS. Consequently, MMS can be considered as equally useful proxies for genetic relatedness as craniometric data. Controlling for population structure, Roseman (2004) found a series of measurements shaped by natural selection, such as cranial breadth measures, nasal height, and measures of midfacial prognathism, all of which increase in cold environments. These measurements are features known to be responses for thermoregulation (Beals 1983, 1984; Franciscus 1991; Krovitz 2003; Roseman 2004). Katz and colleagues (2016) built upon Roseman also removing population structure to examine cold climate effects on the same craniometric dataset using a new statistical approach, BSFG. This study reported consistent patterns of associations between those traits significant in 84 Roseman (2004). Consistencies between these findings and those in this dissertation are the increasing facial prognathism in cold climates, seen in trends for ANS. I suspect that NBC will corroborate this pattern with high nasal bridges associated with cold environments when using original trait score frequencies, which is apparent in frequency tables provided in Hefner and Linde (2018). Another consistent finding is the lack of significance in nasal breadth measures and NAW. Roseman (2004) does indicate that nasal height and nasal breadth negatively covary suggesting that breadth does play an overall role in nasal surface-volume ratios, but the correspondence between nasal breath and climate may be masked by genetics. Finally, IOB does appear to be correlated with cold temperatures, but this is not reported as meaningful through dacryon-dacryon measurements in either Roseman (2004) or Katz et al. (2016). It should be noted though that Maddux et al. (2017) found interorbital breadth to be an important measure for understanding selective pressures on the internal nasal fossa in cold-dry environments. Overall, there does appear to be some consistency between the facial macroscopic morphology and suites of metrics for cranial size and shape. Cranial Nonmetric Trait Heritabilities Although traditional cranial nonmetric traits and MMS traits have not been directly linked, nonmetric trait are separate series of cranial morphological features used to examine and measure population history and biodistance (Berry and Berry 1967; Lane and Sublett 1972; Ossenberg 1976; Milne et al. 1983; Hauser and De Stefano 1989; Saunders 1989; Hanken and Hall 1993; Ishida and Dodo 1993; Hallgrimsson et al. 2004; Hanihara et al., 2003; Sutter and Mertz 2004). Cranial modularity during growth and development make it likely that these two lines of macroscopic analyses of cranial variation overlap considerably. Therefore, we can use previous models of genetic and environmental effects on traditional nonmetric traits to infer 85 about MMS evolvability. Discussing correspondence between these biodistance methods can also serve as a starting point for testing relationships across the two series of traits. Nonmetric traits have been used as genetic proxies for population, operating under the assumption that these traits are inherited and have little environmental influence (Berry and Berry 1967; Sjøvold 1973; Ossenberg 1976; Saunders 1989). Most studies of nonmetric trait heritability used dichotomized scores for presence or absence (Carson 2006), like this dissertation research. Early works on cranial nonmetric trait heritability in skeletal assemblages estimated a wide range of h2 estimates (0.008 to 0.954) (Sjøvold 1973), as well as generally high rates of heritability across studies (Cheverud and Buikstra 1981; Richtsmeier and McGrath 1986). Carson (2006) calculated heritability of nonmetric traits using pedigree data for an Austrian skeletal collection. These heritabilities were calculated using maximum likelihood estimates incorporating influence of sex and birth year. She estimated h2 values between 0 and 1 with most traits ranging from h2 = 0 - 0.4. One trait that appears to correspond between these two macroscopic approaches is zygomaxillary tubercle size (traditional) and MT (MMS). Estimates for zygomaxillary tubercle size is approximately h2 = 0.341 for quasi-continuous scoring and h2 = 0.492 using a dichotomized approach (Carson 2006). This indicates that MT (h2 = 0.24) and zygomaxillary tubercle are amongst higher heritabilities in both series of morphological traits. Both methods indicate low to moderate heritability. Most notably, Carson (2006) found that heritability estimates all decreased significantly when traits are dichotomized; thus, given we see low levels of heritability in MMS traits when translated into presence/absence scores, we can expect these levels to increase when they are assessed using the original quasi-continuous scores. Obviously, environmental factors also a play role in variation of cranial nonmetric traits given they develop during ontogeny, as do MMS traits. These factors can include climate, 86 nutrition, and other deformities and pathological conditions that may influence the development of other features (Ossenberg 1970; Trinkhaus 1978; Saunders 1989; O’Loughlin 2004). These effects are often examined through asymmetry of trait expressions (Trinkhaus 1978; Tagaya 2021), but the effects of climate are still largely unknown. However, discussions of these environmental factors on cranial nonmetric variation suggests that MMS evolutionary models should explore the influences of environmental factors beyond climate, such as nutritional stress and natural and artificial cranial modifications. Limitations and Future Directions The major limitations of this study include dichotomizing of the traits, loss of variation due to regional population assignments, unmatched genetic and phenotypic datasets, and temporal variation between skeletal collections. These limitations are outlined here with suggestions for how to resolve these concerns in future research. Dichotomized Traits As mentioned throughout this chapter, dichotomizing MMS traits from their original quasi-continuous scale reduce the ability to capture the ranges in trait expression. Traits were dichotomized for this project to simplify analyses given this was the first attempt in examining these data in this way. Scoring the traits as present/absent diminishes the trait variance ultimately used to calculate heritability. Thus, heritability estimates were not accurate values and could only be used as relative information. This method also made it impossible to examine actual trait intermediacy (e.g., a score of “2” on a scale of 1-3, like IOB and ANS) when examining the influence of annual climate fluctuations. The impact of climate variability had to be assessed in relation the temperature and humidity effects, such as the skewed of IOB to wider expression when climate is more variable than in cold environments. 87 Future research should rely on original trait scores. Software like SOLAR for MLE heritability estimates are capable of handling ordinal data while also measuring the influence of confounding factors (Carson 2006). The BSFG model in MATLAB can also be used for ordinal data. The use of these ordinal data are feasible and will more accurately estimate phenotypic variances of MMS traits. Regional Population Assignments Large scale regional population assignments were used to correspond with assignments from similar research for more appropriate comparisons, comprehend these data more easily on a large scale, and reliably overlap geographic regions between MMS and microsatellite data. However, grouping various populations on this large scale likely aggregates a range of variation and climates. For example, you may have two groups that are in close geographic proximity, but one population lives in high elevations while the other lives in low elevations. This would likely result in different selective pressures leading to some degree of phenotypic variation. If selection is of particular interest, more geographically refined models need be applied in the future. Unmatched Datasets Another major limitation of this study is the lack of unmatched datasets. Microsatellite data is very useful as proxies for genetic variation within and between populations. The nature of skeletal collections makes it difficult to obtain genetic information specific to your sample given these analyses are destructive; thus, these global genetic datasets provide a means for modeling general genetic trends for skeletal data. Nevertheless, there will be some degree of discordance between genetic and skeletal data that results in somewhat skewed correlations between genetic and phenotypic distances due to sample biases. Patterns of gene flow is greatly affected by assortative mating practices and other geographic, social, and political barriers which make it 88 nearly impossible to select the same sample without obtaining genetic information directly from the skeletal individuals. Future research should include more refined studies on skeletal samples with either pedigree or matched genetic information to (1) accurately calculate heritability of MMS traits and (2) accurately control for genetic effects in climate models. Alternatively, geographic distances can be used as genetic proxies for large scale studies. Temporal Variation between Skeletal Collections A final limitation of this study is the temporal variation between the MMS populations, as well as between the MMS and genetic populations. All MMS populations are from archaeological contexts, except those from Southwest America, Colombia, and Thailand which are contemporary (19th and 20th century birth years). The microsatellite data are also obtained from contemporary humans. As groups of people move across the landscape and interact, their population structures and phenotypic expressions may be altered. Additionally, over time, populations will experience new genetic mutations and adaptations to their local environments, which can lead to changes in phenotypic expression (West-Eberhard 2008; Relethford 2019). Thus, we expect secular change in MMS trait expression across groups, as well as different rates of change due to population histories. This likely means some level of discordance between the MMS and genetic data, as well across MMS groups. Future work should aim to use only contemporary phenotypic data to best reflect modern population structures and available genetic datasets. Additionally, reconstruction of paleoclimates would improve results for archaeological skeletal materials with more accurate climate models for the material being examined. Incorporating paleoclimate reconstructions along with hominin fossil data, such as cranial and genome data from Neandertal and Denisovan specimens, expands the ability to observe 89 evolutionary processes over time and how remnants of genetic admixture from these hominin ancestors play a role in modern phenotypic expression. For example, Neanderthals found in Siberia are documented as having wide crania and greater facial prognathism (Krovitz 2003; Gregory et al. 2017). While the genes responsible for these traits are reduced in anatomically modern humans (Gokhman, et al. 2020), these patterns of phenotypic expression exist in modern Siberian populations. This may suggest Neanderthal admixture continues to influence cranial form today given there is evidence of greater admixture documented in these regions (Simonti et al. 2016). However, the cause of these changes in genetic expression are still unknown and the geographic patterning of trait expressions make it difficult to postulate about the evolutionary mechanisms. Applying the BSFG modeling to examine evolutionary forces of phenotypic changes between hominin fossils and anatomically modern humans would help to understand the reduction of these traits over time, as well as how admixture from these early ancestors influence cranial form today. Broader Evolutionary Significance The evolutionary mechanisms behind craniometric have been explored rather extensively (Relethford and Harpending 1994; Relethford 2002, 2004a, 2004b; Gonzalez-Jose et al 2004; Roseman 2004; Roseman and Weaver 2004, 2007; Harvati and Weaver 2006; Manica et al. 2007; von Cramon-Taubadel and Lycett 2008; Betti et al. 2009; Hubbe et al. 2009; Smith 2009; von Cramon-Taubadel 2009; von Cramon-Taubadel and Weaver 2009; Reyes-Centeno et al. 2014; Katz et al. 2016). Craniometric studies use a full suite of traits that capture overall size and shape of the human cranium, while MMS traits focus on macroscopic morphological assessments primarily on the midfacial skeleton with only one trait on the cranial vault. Most significant craniometric features significantly associated with climate are vault measures, making 90 MMS a unique approach to exploring global patterns of phenotypic variation and evolution. MMS traits are particularly important for expanding our understanding of natural selection whereby a large portion of cranial evolutionary research has centered around the neutral evolutionary processes. The wealth of research demonstrating that the middle face is highly responsive to climate due to respiratory stress emphasizes the importance of exploring the proportion of genetics and environments on MMS traits, given the heavy focus on the midfacial skeleton, as well as their utility in addressing these larger questions. Additionally, many craniometric studies rely on the Howells dataset due to global coverage of these data. The MaMD provides a new global dataset to compare to spatially patterned craniometric variation, as well as offers unique population datasets for biodistance and evolutionary research. Evolutionary research has moved well beyond Darwinian theory which oversimplifies evolution-by-phenotype processes and led to reductionism in Western biological sciences, or the idea that everything can be explained using a genetic model (Weiss and Buchanon 2004). The scientific community generally accepts that evolution is a complex relationship between genotypes and phenotypes due to modularity, the role of non-coding DNA, polymorphism, and many other behaviors of genes; and, moreover, that nongenetic forces, like the environment, also play an important role (Weiss and Buchanon 2004; Noble 2015). Therefore, we need to model these complex relationships in our evolutionary and human variation studies of skeletal remains and refrain from falling back to overly simplistic Darwinian models. The genetic data selected for this research helps to tease apart different evolutionary mechanisms. Microsatellite data is ideal for detecting gene flow and genetic drift (i.e., neutral forces of evolution) but are not good agents for detecting natural selection processes (Haasl et al. 2013). Therefore, when using microsatellite data to test genetic relatedness of MMS expression, we are marginalizing error in 91 misinterpreting evolutionary mechanisms. The geographically patterned nature of both selective environmental pressures and genetics due to neutral forces means that removing population structure should reveal selective pressure of climate (Roseman 2004). Theoretically, MMS traits with lower correspondence to microsatellite data are likely not due to neutrality and, thus, have greater response to natural selection and plasticity. MMS trait data combined with microsatellite data and methods for removing population structure effects in climate models furthers our efforts to understand the role of natural selection in craniofacial morphology. Another unique perspective of this study building on our understanding of evolution in craniofacial morphology is the exploration of influence of annual climate fluctuations on selective pressures. Katz et al. (2016) and Roseman (2004) noted that Siberia was driving many of their significant associations between craniometrics and cold climates. Other studies also indicate that climate correlations are significant only when models include populations from very cold regions (Harvati and Weaver 2006; Hubbe et al. 2009; Betti et al. 2010). This supports thermoregulatory studies stating that cold-dry environments cause greatest respiratory stress and signal phenotypic responses to improve survivability while hot-wet environments reduce selective pressure (Walker et al. 1961; Franciscus & Long 1991; Noback et al. 2011; Lieberman 2015; Maddux et al. 2017). This reduced pressure is evident in patterns of skewed trait expression when climates are more variable. For example, ANS and IOB, which shape the overall nose morphology (projection of the nose and breadth of internal nasal fossa), demonstrate high and narrow noses in cold environments. However, in more variable climate regions, the nose becomes less pronounced from the face and wider like a result of less selective pressure. Incorporating these various environment experiences across regional populations also aids in our ability to disentangle selection from neutral evolution. 92 Finally, the MMS method is frequently used in forensic anthropological research and casework. This project continues to root the method deep withing evolutionary theory improving the reliability and utility of MMS traits. Forensic anthropology often relies on patterns of variation to answer questions about the unknown, yet we often become very detached from our evolutionary roots. Improving identification efforts in forensic investigations will require more than a name change from “race” to “ancestry” to “population affinity”. It requires our methods to be situated back within evolutionary theory to truly capture what it is we aim to answer. Population affinity estimates aim to decipher to whom the unidentified individual was most closely related at the most refined population possible. To do this, we aim to capture genetic relatedness through phenotypic variation. This genetic relatedness is shaped by our parent’s migration history, sociopolitical boundaries, language barriers, geographic barriers, and various microcultural and behavioral actions. The accuracy of our methods rely on our ability to properly relate phenotypic and genetic variation while understanding the confounding factors such as natural selection. Therefore, establishing evolutionary foundations of our methods and the features we are using to measure phenotypic variation is critical to moving toward more evolutionarily informed estimates of population affinity for unidentified human remains. 93 APPENDIX 94 The macromorphoscopic trait populations used in this dissertation are described below. These data were all obtained from the Macromorphoscopic Databank. The information for each regional population provided below includes the local populations, the physical locations where collections are currently housed, whether the collections are either archaeological or contemporary, who collected the macromorphoscopic trait data, and, when known, how the collections were obtained by the current institution. Eastern Africa The Eastern African sample includes individuals from Kenya (n=11), Rwanda (n=6), Somalia (n=8), Tanzania (n=8), and Uganda (n=5), all of which are from archaeological contexts. These data were collected by Amber Plemons from individuals housed at the Natural History Museum in London. East Asia The East Asian population comprises individuals from China (n=59) and Japan (n=15) housed at the National Museum of Natural History (NMNH), Smithsonian Institution in Washington, DC. These individuals are from archaeological contexts. Chinese individuals are from a cannery in Karluk, Alaska and donated by Hrdlicka, and Japanese individuals are from Tokyo and donated by Tokyo Imperial University. Data were collected by Dr. Joseph Hefner. Europe The European sample consists of individuals from Pedersöre, Finland (n=165) and are currently housed at the Finnish Museum of Natural History in Helsinki, Finland. These individuals were associated with the Pedersöre Church from the medieval time period to the 1800’s. Data were collected by Amber Plemons. 95 Latin America The Latin American sample comprises individuals from the Southwest United States (Southwest Hispanic) (n=101), Colombia (n=244), and Peru (n=50). The Southwest Hispanic population includes contemporary individuals housed at Texas State University and the Pima County Medical Examiner’s Office. These data were collected by Amber Plemons and Dr. Caitlin Vogelberg. The Colombian individuals are also contemporary, and data were collected by Timisay Monsalve. Individuals from Peru are from archaeological contexts and were collected by Nicole Jastremski. Pacific Island The Pacific Island sample includes individuals from Australia (43), Fiji (13), French Polynesia (6), Indonesia (40), Malaysia (5) New Zealand - Maori (7), New Zealand (27), Papau New Guinea (46), and the Philippine Islands (35). These data are currently housed at the University of Pennsylvania Museum of Archaeology and Anthropology and the Smithsonian NMNH. Data were collected by Melanie Ratliff. Southeast Asia The Southeast Asian sample comprises individuals from Thailand (n=756) housed at Khon Kaen University in Khon Kaen, Thailand. These are contemporary individuals who donated their body to the university for research, primarily residents of Northeast Thailand. Data were collected by Dr. Joseph Hefner, Amber Plemons, and Kelly Kamnikar. Siberia The Siberian sample includes Aleut (n=168), Barrow (n=8), and St. Lawrence Inuit (n=9). These individuals are from archaeological contexts (pre- and protohistoric) and are housed 96 at the NMNH, Smithsonian Institution in Washington, DC. Data were collected by Dr. Joseph Hefner. South Africa The South African sample includes individuals from Malawi (n=19) and South Africa (n=26) and are from archaeological contexts. These data were collected by Amber Plemons from individuals housed at the Natural History Museum in London. West Africa The West African sample includes individuals from Gambia (n=3) and Nigeria (n=60) and are from archaeological contexts. These data were collected by Amber Plemons from individuals housed at the Natural History Museum in London. 97 LITERATURE CITED 98 LITERATURE CITED Humidity Conversion Formulas: Calculation formulas for humidity. 2013. Finland: Vaisala Oyj. Agarwal S, and P Beauchesne. 2011. It is not carved in bone: development and plasticity of the aged skeleton. In Social bioarchaeology, edited by S. Agarwal and B. Glencross. New York: Wiley-Blackwell. Agarwal SC. 2016. Bone morphologies and histories: Life course approaches in bioarchaeology. American Journal of Physical Anthropology 159 (S61):130-149. Allen JA. 1877. The influence of physical conditions in the genesis of species. Radical Review 1:108-140. Auerbach BM. 2012. Skeletal variation among early Holocene North American humans: implications for origins and diversity in the Americas. American Journal of Physical Anthropology 149 (4):525. Beals KL, CL Smith, and SM Dodd. 1983. Climate and the evolution of brachycephalization. American Journal of Physical Anthropology 62 (4):425-37. Beals KL, CL Smith, and SM Dodd. 1984. Brain size, cranial morphology, climate, and time machines. Current Anthropology 25 (3):301-330. Berger SM, JS Griffin, and SC Dent. 2021. Phenotypes and pathways: Working toward an integrated skeletal biology in biological anthropology. American Journal of Human Biology 33 (2):e23450. Bergmann C. 1847. Ueber die Verh€altnisse der W€arme€okonomie der Thiere zu ihrer Gr€osse. Gottinger Studien 3:595-708. Berry CA, and RJ Berry. 1967. Epigenetic variation in the human cranium. Journal of Anatomy 101 (Pt 2):361-379. Bethard JD, and EA DiGangi. 2020. Letter to the Editor—Moving Beyond a Lost Cause: Forensic Anthropology and Ancestry Estimates in the United States. Journal of Forensic Sciences 65 (5):1791-1792. Betti L, F Balloux, W Amos, T Hanihara, and A Manica. 2009. Distance from Africa, not climate, explains within-population phenotypic diversity in humans. Proc R Soc Lond B Biol Sci 276:809-814. Betti L, F Balloux, T Hanihara, and A Manica. 2010. The relative role of drift and selection in shaping the human skull. American Journal of Physical Anthropology 141:76-82. 99 Bigham AW, X Mao, R Mei, T Brutsaert, MJ Wilson, CG Julian, EJ Parra, JM Akey, LG Moore, and MD Shriver. 2009. Identifying positive selection candidate loci for high-altitude adaptation in Andean populations. Human Genomics 4 (2):79-90. Blackburn TM, KJ Gaston, and N Loder. 1999. Geographic Gradients in Body Size: A Clarification of Bergmann's Rule. Diversity and Distributions 5 (4):165-174. Blois JL, JW Williams, MC Fitzpatrick, ST Jackson, and S Ferrier. 2013. Space can substitute for time in predicting climate-change effects on biodiversity. Proceedings of the National Academy of Sciences of the United States of America 110 (23):9374-9379. Boas F. 1911. Changes in Bodily Form of Descendants of Immigrants. Washington, DC: Government Printing Office. Boas F. 1912. Changes in Bodily Form of Descendants of Immigrants. New York: Columbia University Press. Boas F. 1940. Race, Language, and Culture. New York: Macmillan. Bonduriansky R. 2021. Plasticity across Generations. In Phenotypic Plasticity & Evolution: Causes, Consequences, Controversies, edited by D. Pfenning. Boca Raton, FL: CRC Press. Borchers HW. 2021. Practical Numerical Math Functions. Buikstra JE, SR Frankenberg, and LW Konigsberg. 1990. Skeletal biological distance studies in American Physical Anthropology: Recent trends. American Journal of Physical Anthropology 82 (1):1-7. Calafell F and MHD Larmuseau. 2017. The Y chromosome as the most popular marker in genetic genealogy benefits interdisciplinary research. Human Genetics 136 (5):559-573. Carey G. 2003. Human Genetics for the Social Sciences. New York: Sage Publications, Inc. Carson EA. 2006. Maximum-Likelihood Variance Components Analysis of Heritabilities of Cranial Nonmetric Traits. Human Biology 78 (4):383-402. Chaix R, F Austerlitz, T Khegay, S Jacquesson, MF Hammer, E Heyer, and L Quintana-Murci. 2004. The Genetic or Mythical Ancestry of Descent Groups: Lessons from the Y Chromosome. American Journal of Human Genetics 75 (6):1113-1116. Chenard KC, and RA Duckworth. 2021. The Special Case of Behavioral Plasticity? In Phenotypic Plasticity & Evolution: Causes, Consequences, Controversies, edited by D. Pfenning. Boca Raton, FL: CRC Press. Cheverud JM. 1988. A Comparison of Genetic and Phenotypic Correlations. Evolution 42 (5):958-968. 100 Churchill SE, LL Shackelford, JN Georgi, and MT Black. 2004. Morphological variation and airflow dynamics in the human nose. American Journal of Human Biology 16 (6):625-638. Claus J. 2001. Selective Brain Cooling in Mammals and Birds. The Japanese Journal of Physiology 51 (3):291. Cole P. 1982. Modification of inspired air. In The nose: Upper airway physiology and the atmospheric environment, edited by D. Proctor and I. Andersen. New York: Elsevier Biomedical Press. Cole P. 2000. Biophysics of Nasal Airflow: A Review. American Journal of Rhinology 14 (4):245-250. Cooper G, W Amos, R Bellamy, MR Siddiqui, A Frodsham, AVS Hill, and DC Rubinsztein. 1999. An Empirical Exploration of the (Δμ) 2 Genetic Distance for 213 Human Microsatellite Markers. The American Journal of Human Genetics 65 (4):1125-1133. Crawford MH. 2007. Foundations of Anthropological Genetics. In Anthropological Genetics: Theory, Methods, and Applications, edited by M Crawford. Cambridge: Cambridge University Press. Crawford MH. 2018. Anthropological genetics. In The International Encyclopedia of Biological Anthropology. Crawford MH. 2019. A Companion to Anthropological Genetics, edited by DH O'Rourke. Newark, United States: John Wiley & Sons, Incorporated. Cushman JH, and JH Lawton. 1993. Latitudinal Patterns in European Ant Assemblages: Variation in Species Richness and Body Size. Oecologia 95 (1):30-37. Devlin MJ. 2015. The “Skinny” on brown fat, obesity, and bone. American Journal of Physical Anthropology 156 (S59):98-115. DiGangi EA, and JD Bethard. 2021. Uncloaking a Lost Cause: Decolonizing ancestry estimation in the United States. American Journal of Physical Anthropology 175 (2):422-436. Dumache R, V Ciocan, C Muresan, and A Enache. 2016. Molecular DNA Analysis in Forensic Identification. Clinical laboratory 62 (1-2):245-248. Dumache R, V Ciocan, C Muresan, and A Enache. 2016. Molecular Genetics and its Applications in Forensic Sciences. In Forensic Analysis - From Death to Justice, edited by BSK Shetty and JR Padubidri: IntechOpen. Edgar H and M Pilloud. 2021. A Reassessment of Assessing Race: "Ancestry" Estimation and Its Implications for Forensic Anthropology and Beyond. Forensic Anthropology 4 (4):67-72. 101 Evteev A, AL Cardini, I Morozova, and P O'Higgins. 2014. Extreme climate, rather than population history, explains mid-facial morphology of Northern Asians. American Journal of Physical Anthropology 153 (3):449. Farkas LG, JC Kolar, and IR Munro. 1986. Geography of the nose: A morphometric study. Aesthetic Plastic Surgery 10 (1):191-223. Fausto-Sterling A. 2005. The Bare Bones of Sex: Part 1—Sex and Gender. Signs 30 (2):1491-1527. Ferrario VF, CS, CE Poggio, and JH Schmitz. 1997. Three-dimensional study of growth and development of the nose. Cleft Palate Craniofacial Journal 34 (4):309-17. Fox Rebecca J., Jennifer M. Donelson, Celia Schunter, Timothy Ravasi, and Juan D. Gaitán- Espitia. 2019. Beyond buying time: the role of plasticity in phenotypic adaptation to rapid environmental change. Philosophical Transactions of the Royal Society B: Biological Sciences 374 (1768):20180174. Franciscus RG. 1995. Later Pleistocene nasofacial variation in western Eurasia and Africa and modern human origins. Ph.D. thesis, University of New Mexico, Albuquerque. Franciscus RG. 2003. Internal nasal floor configuration in Homo with special reference to the evolution of neandertal facial form. Journal of Human Evolution 44:701-729. Franciscus RG, and JC Long. 1991. Variation in human nasal height and breadth. American Journal of Physical Anthropology 85 (4):419-27. Friedlaender JS, FR Friedlaender, FA Reed, KK Kidd, JR Kidd, GK Chambers, RA Lea, JH Loo, G Koki, JA Hodgson, DA Merriwether, and JL Weber. 2008. The genetic structure of Pacific Islanders. PLoS Genet 4 (1):e19. Frye AJ, and E Kamon. 1983. Sweating efficiency in acclimated men and women exercising in humid and dry heat. Journal of Applied Physiology Respir Environ Exerc Physiol 54 (4):972-7. Fuentes A. 2010. The New Biological Anthropology: Bringing Washburn’s New Physical Anthropology Into 2010 and Beyond—The2008 AAPA Luncheon Lecture. Yearbook of Physical Anthropology 53 (2):2-12. Futuyma DJ. 2021. How Does Phenotypic Plasticity Fit into Evolutionary Theory? In Phenotypic Plasticity & Evolution: Causes, Consequences, Controversies, edited by D. Pfenning. Boca Raton, FL: CRC Press. Garland T. Jr and SA Kelly. 2018. Phenotypic plasticity and experimental evolution. The Journal of Experimental Biology 209:2344-2361. Giles E, S Wyber, and RJ Walsh. 1970. Microevolution in New Guinea: Additional Evidence for Genetic Drift. Archaeology & Physical Anthropology in Oceania 5 (1):60-72. 102 Gokcumen O. 2018. The Year In Genetic Anthropology: New Lands, New Technologies, New Questions. American Anthropologist 120 (2):266-277. Goldstein DB, AR Linares, LL Cavalli-Sforza, and MW Feldman. 1995a. An evaluation of genetic distances for use with microsatellite loci. Genetics 139 (1):463-471. Goldstein DB, AR Linares, LL Cavalli-Sforza, and MW Feldman. 1995b. Genetic absolute dating based on microsatellites and the origin of modern humans. Proceedings of the National Academy of Sciences 92 (15):6723-6727. González-José R, S Van Der Molen, E González-Pérez, and M Hernández. 2004. Patterns of phenotypic covariation and correlation in modern humans as viewed from morphological integration. American Journal Physical Anthropology 123 (1):69-77. Gosman JH, SD Stout, and CS Larsen. 2011. Skeletal biology over the life span: A view from the surfaces. American Journal of Physical Anthropology 146 (S53):86-98. Gravlee CC, HR Bernard, and WR Leonard. 2003. Boas's "Changes in Bodily Form": The Immigrant Study, Cranial Plasticity, and Boas's Physical Anthropology. American Anthropologist 105 (2):326-332. Gravlee CC, HR Bernard, and WR Leonard. 2003. Heredity, Environment, and Cranial Form: A Reanalysis of Boas's Immigrant Data. American Anthropologist 105 (1):125-138. Guo J, J Tan, Y Yang, H Zhou, S Hu, A Hashan, N Bahaxar, S Xu, TD Weaver, L Jin, M Stoneking, and K Tang. 2014. Variation and signatures of selection on the human face. Journal of Human Evolution 75:143-52. Haasl RJ, and BA Payseur. 2012. Microsatellites as Targets of Natural Selection. Molecular Biology and Evolution 30 (2):285-298. Hallgrímsson B, BO Donnabháin, GB Walters, DM Cooper, D Gudbjartsson, and K Stefánsson. 2004. Composition of the founding population of Iceland: biological distance and morphological variation in early historic Atlantic Europe. American Journals of Physical Anthropology 124 (3):257-74. Hamilton M. 2009. Population Genetics. Hoboken, New Jersey: Wiley-Blackwell. Hanihara T, H Ishida, and Y Dodo. 2003. Characterization of biological diversity through analysis of discrete cranial traits. American Journal of Physical Anthropology 121 (3):241-251. Hanken J, and BK Hall. 1993. The Skull. Chicago: University of Chicago Press. Harris I, TJ Osborn, P Jones, and D Lister. 2020. Version 4 of the CRU TS monthly high- resolution gridded multivariate climate dataset. Scientific Data 7 (1):109. 103 Harvati K, and TD Weaver. 2006. Human cranial anatomy and the differential preservation of population history and climate signatures. The Anatomical Record Part A: Discoveries in Molecular, Cellular, and Evolutionary Biology 288A (12):1225-1233. Hauser G, and G De Stefano. 1989. Epigenetic Variants of the Human Skull. Schweizerbart: Stuttgart. Havenith G. 2005. Temperature Regulation, Heat Balance and Climatic Stress. In Extreme Weather Events and Public Health Responses, edited by W Kirch, R Bertollini and B Menne. Berlin, Heidelberg: Springer Berlin Heidelberg. Hefner JT. 2009. Cranial Nonmetric Variation and Estimating Ancestry. Journal of Forensic Sciences 54 (5):985-995. Hefner JT. 2018. The macromorphoscopic databank. American Journal of Physical Anthropology 166 (4):994-1004. Hefner JT, and KC Linde. 2018. The Macromorphoscopic Databank. San Diego: Elsevier, Academic Press. Hefner JT, and SD Ousley. 2014. Statistical Classification Methods for Estimating Ancestry Using Morphoscopic Traits. Journal of Forensic Sciences 59 (4):883-890. Hefner JT, MA Pilloud, CJ Black, and BE Anderson. 2015. Morphoscopic Trait Expression in Hispanic Populations. Journal of Forensic Sciences 60 (5):1135-1139. Hefner JT, MA Pilloud, J Buikstra, and C Vogelsberg. 2016. A Brief History of Biological Distance Analysis. In Biological Distance Analysis: Forensic and Bioarchaeological Perspectives, edited by MA Pilloud and JT Hefner. London: Academic Press. Holliday TW. 1997. Postcranial evidence of cold adaptation in European Neandertals. American Journal of Physical Anthropology 104 (2):245. Holton NE, TR Yokley, and A Figueroa. 2012. Nasal septal and craniofacial form in European- and African-derived populations. Journal of Anatomy 221 (3):263-274. Honnay O. 2013. Genetic Drift. In Brenner's Encyclopedia of Genetics (Second Edition), edited by S. Maloy and K. Hughes. San Diego: Academic Press. Howells WW. 1973. Cranial variation in man: a study by multivariate analysis of patterns of difference among recent human populations. Papers of the Peabody Museum No. 67. Cambridge, MA: Peabody Museum. Howells WW. 1989. Skull shapes and the map: craniometric analyses in the dispersion of modern Homo. Papers of the Peabody Museum No. 79. Cambridge, MA: Peabody Museum. Howells WW. 1995. Who’s who in skulls: ethnic identification of crania from measurements. Papers of the Peabody Museum No. 82. Cambridge, MA: Peabody Museum. 104 Howells WW. 1996. Howells’ craniometric data on the internet. Am J Physical Anthropology 101:441-442. Hubbe M, T Hanihara, and K Harvati. 2009. Climate signatures in the morphological differentiation of worldwide modern human populations. Anatomical record 292 (11):1720- 1733. Inthavong K, J Ma, Y Shang, J Dong, ASR Chetty, J Tu, and D Frank-Ito. 2019. Geometry and airflow dynamics analysis in the nasal cavity during inhalation. Clinical Biomechanics 66:97- 106. Inthavong K, ZF Tian, and J Tu. 2007. CFD Simulations on the Heating Capability in a Human Nasal Cavity. 16th Australasian Fluid Mechanics Conference Gold Coast, Australia. Ishida H, and Y Dodo. 1993. Nonmetric cranial variation and the population affinities of the pacific peoples. American Journal of Physical Anthropology 90:49-57. Ives R, and L Humphrey. 2018. Endochondral growth disruption during vitamin D deficiency rickets in a mid-19th century series from Bethnal Green, London, UK. American Journal of Physical Anthropology 167 (3):585-601. Jackson FLC, and Borgelin, LFJ. 2010. How Genetics Can Provide Detail to the Transatlantic African Diaspora. In The African Diaspora and the Disciplines, edited by T Olaniyan, JH Sweet. Bloomington: Indiana University Press. Katz DC, MN Grote, and TD Weaver. 2016. A mixed model for the relationship between climate and human cranial form. American Journal of Physical Anthropology 160 (4):593. Katzmarzyk PT, and WR Leonard. 1998. Climatic influences on human body size and proportions: ecological adaptations and secular trends. American Journal of Physical Anthropology 106 (4):483. Keck T, and J Lindemann. 2010. Numerical simulation and nasal air-conditioning. GMS Current Topics Otorhinolaryngology Head Neck Surg 9:Doc08. Kelly M. 2019. Adaptation to climate change through genetic accommodation and assimilation of plastic phenotypes. Philosophical Transactions of the Royal Society B: Biological Sciences 374 (1768):20180176. Kenyhercz M, NV Passalacqua, and JT Hefner. 2019. Missing Data Imputation Using Morphoscopic Traits and Their Performance in the Estimation of Ancestry. Forensic Anthropology 2 (3):178-188. Kingsolver JG, and SE Diamond. 2011. Phenotypic Selection in Natural Populations: What Limits Directional Selection? The American Naturalist 177 (3):346-357. Klales AR, and MW Kenyhercz. 2015. Morphological Assessment of Ancestry using Cranial Macromorphoscopics. Journal of Forensic Sciences 60 (1):13-20. 105 Klepinger LL. 2006. Fundamentals of Forensic Anthropology. Edited by M. Cartmill and K. Brown. Hoboken, New Jersey: John Wiley & Sons, Inc. Kopelman NM, L Stone, C Wang, D Gefel, MW Feldman, J Hillel, and NA Rosenberg. 2009. Genomic microsatellites identify shared Jewish ancestry intermediate between Middle Eastern and European populations. BMC Genetics 10:80. Kottack CP. 2002. Cultural Anthropology. 9th ed. New York: McGraw Hill. Kowalczyk M, E Zawadzka, D Szewczuk, M Gryzińska, and A Jakubczak. 2018. Molecular markers used in forensic genetics. Medicine, science, and the law 58 (4):201-209. Krovitz GE. 2003. Shape and growth differences between Neandertals and modern humans: Grounds for a species-level distinction? In Patterns of Growth and Development in the Genus Homo, edited by J Thompson, G Krovitz and A Nelson. Kuzawa CW, and Bragg JM. 2012. Plasticity in Human Life History Strategy: Implications for Contemporary Human Variation and the Evolution of Genus Homo. Current Anthropology 53 (6):369-382. Lane RA, and AJ Sublett. 1972. Osteology of Social Organization: Residence Pattern. American Antiquity 37 (2):186-201. Lasker GW. 1961. The Evolution of Man. New York: Holt Rinehart & Winston. Lasker GW. 1969. Human biological adaptability. The ecological approach in physical anthropology. Science 166 (2913):1480-1486. Levis NA, and DW Pfenning. 2021. Innovation and Diversification Via Plasticity-Led Evolution. In Phenotypic Plasticity & Evolution: Causes, Consequences, Controversies, edited by D Pfenning. Boca Raton, FL: CRC Press. Lewontin RC. 1967. An estimate of average heterozygosity in man. American Journal of Human Genetics 19 (5):681-685. Lewontin RC. 1972. The Apportionment of Human Diversity. In Evolutionary Biology: Volume 6, edited by T Dobzhansky, MK Hecht and WC Steere. New York, NY: Springer US. Li YC, AB Korol, T Fahima, and E Nevo. 2004. Microsatellites within genes: structure, function, and evolution. Molecular Biology and Evolution 21 (6):991-1007. Lieberman DE. 2015. Human locomotion and heat loss: An evolutionary perspective. Comprehensive Physiology 5:99-117. Livingston FB. 1971. Malaria and Human Polymorphisms. Annual Review of Genetics (5):33- 64. 106 Maddux SD, TR Yokley, RG Franciscus, and BM Svoma. 2016. Absolute humidity and the human nose: A reanalysis of climate zones and their influence on nasal form and function. American Journal of Physical Anthropology. Maddux SD, LN Butaric, TR Yokley, and RG Franciscus. 2017. Ecogeographic variation across morphofunctional units of the human nose. American Journal of Physical Anthropology 162 (1):103-119. Madrigal L, and G Barbujani. 2007. Partitioning of Genetic Variation in Human Populations and the Concept of Race. In Anthropological Genetics: Theory, Methods, and Applications, edited by MH Crawford. Cambridge: Cambridge University Press. Manica A, W Amos, F Balloux, and T Hanihara. 2007. The effect of ancient population bottlenecks on human phenotypic variation. Nature 448:346-348. Mantel N. 1967. The detection of disease clustering and a generalized regression approach. Cancer Research 27 (2):209-20. Mariak Z, MD White, J Lewko, T Lyson, and P Piekarski. 1999. Direct cooling of the human brain by heat loss from the upper respiratory tract. Journal of Applied Physiology 87 (5):1609- 13. Markel AL. 2014. Physiological genetics. Russian Journal of Genetics: Applied Research 4 (4):301-310. Marks J. 2009. What is the Viewpoint of Hemoglobin, and Does It Matter? History and Philosophy of the Life Sciences 31 (2):241-262. Marks J. 2012. The Origins of Anthropological Genetics. Current Anthropology 53 (S5):S161- S172. Martínez-Abadías N, M Esparza, T Sjøvold, R González-José, M Santos, and M Hernández. 2009. Heritability of human cranial dimensions: comparing the evolvability of different cranial regions. Journal of Anatomy 214 (1):19-35. Mayr E. 1956. Geographical character gradients and climatic adaptation. Evolution 10:105-108. McNab BK. 2010. Geographic and temporal correlations of mammalian size reconsidered: a resource rule. Oecologia 164 (1):13-23. Mielke JH, and AG Fix. 2007. The Confluence of Anthropological Genetics and Anthropological Demography. In Anthropological Genetics: Theory, Methods, and Applications, edited by MH Crawford. Cambridge: Cambridge University Press. Milgroom MG. 2017. Population Biology of Plant Pathogens: Genetics, Ecology, and Evolution. St. Paul, Minnesota: The American Pathological Society. 107 Milne N, LH Schmitt, and L Freedman. 1983. Discrete trait variation in Western Australian Aboriginal skulls. Journal of Human Evolution 12 (2):157-168. Montinaro F, V Pankratov, B Yelmen, L Pagani, and M Mondal. 2021. Revisiting the out of Africa event with a deep-learning approach. The American Journal of Human Genetics 108 (11):2037-2051. Moore JL, and JV Remais. 2014. Developmental models for estimating ecological responses to environmental variability: structural, parametric, and experimental issues. Acta Biotheoretica 62 (1):69-90. Moore JL, and JV Remais. 2014. Developmental models for estimating ecological responses to environmental variability: structural, parametric, and experimental issues. Acta Biotheoretica 62 (1):69-90. Moreno-Mayar JV, BA Potter, L Vinner, M Steinrücken, S Rasmussen, J Terhorst, JA Kamm, A Albrechtsen, AS Malaspinas, M Sikora, JD Reuther, JD Irish, RS Malhi, L Orlando, YS Song, R Nielsen, DJ Meltzer, and E Willerslev. 2018. Terminal Pleistocene Alaskan genome reveals first founding population of Native Americans. Nature 553 (7687):203-207. Moreno-Mayar JV, Lasse V, P de Barros Damgaard, C de la Fuente, J Chan, PS Jeffrey, EA Morten, T Vimala, F Racimo, T Pinotti, S Rasmussen, A Margaryan, MI Orbegozo, D Mylopotamitaki, M Wooller, C Bataille, L Becerra-Valdivia, D Chivall, D Comeskey, T Devièse, KG Donald, L George, H Harry, V Alexandersen, C Primeau, J Erlandson, C Rodrigues-Carvalho, S Reis, QRB Murilo, J Cybulski, C Vullo, F Morello, M Vilar, S Wells, K Gregersen, LH Kasper, N Lynnerup, M Mirazón Lahr, K Kjær, A Strauss, M Alfonso-Durruty, A Salas, H Schroeder, T Higham, SM Ripan, TR Jeffrey, L Souza, RS Fabricio, AS Malaspinas, M Sikora, R Nielsen, SS Yun, JM David, and E Willerslev. 2018. Early human dispersals within the Americas. Science 362 (6419):eaav2621. Nesheva DV. 2014. Aspects of ancient mitochondrial DNA analysis in different populations for understanding human evolution. Balkan Journal of Medical Genetics: BJMG 17 (1):5-14. Noback ML, K Harvati, and F Spoor. 2011. Climate-related variation of the human nasal cavity. American Journal of Physical Anthropology 145 (4):599. Nobel D, Hoppeler HH. 2015. Evolution beyond neo-Darwinism: a new conceptual framework. Journal of Experimental Biology 218 (1): 7–13. Oksanen J, F Blanchet, M Friendly, R Kindt, P Legendre, D McGlinn, PR Minchin, RB O'Hara, GL Simpson, P Solymos, MHH Stevens, E Szoecs, H Wagner. 2020. vegan: Community Ecology Package. O'Loughlin VD. 2004. Effects of different kinds of cranial deformation on the incidence of wormian bones. American Journal of Physical Anthropology 123 (2):146-155. Omer G. 2017. Evolution, Function, and Deconstructing Histories: A New Generation of Anthropological Genetics. Human Biology 89 (1):5-6. 108 O'Rourke DH. 2003. Anthropological Genetics in the Genomic Era: A Look Back and Ahead. American Anthropologist 105 (1):101-109. Ossenberg NS. 1970. The influence of artificial cranial deformation on discontinuous morphological traits. American Journal of Physical Anthropology 33 (3):357-371. Ossenberg NS. 1976. Within and between race distances in population studies based on discrete traits of the human skull. American Journal of Physical Anthropology 45 (3):701-715. Pemberton TJ, FY Li, EK Hanson, NU Mehta, S Choi, J Ballantyne, JW Belmont, NA Rosenberg, C Tyler-Smith, and PI Patel. 2012. Impact of restricted marital practices on genetic variation in an endogamous Gujarati group. American Journal of Physical Anthropology 149 (1):92-103. Pemberton TJ, M DeGiorgio, and NA Rosenberg. 2013. Population structure in a comprehensive genomic data set on human microsatellite variation. G3 3 (5):891-907. Perry BW, DR Schield, and TA Castoe. 2018. Evolution: Plasticity versus Selection, or Plasticity and Selection? Current Biology 28 (18):2970-2977. Pfenning DW. 2021. Key Questions about Phenotypic Plasticity. In Phenotypic Plasticity & Evolution: Causes, Consequences, Controversies, edited by D Pfenning. Boca Raton, FL: CRC Press. Pietrusewsky M. 2017. Biological Distance in Bioarchaeology and Human Osteology. In Encyclopedia of Global Archaeology, edited by C Smith. Cham: Springer International Publishing. Plemons AM, JT Hefner, and KR Kamnikar. 2018. Refining Asian Ancestry Classifications via Cranial Macromorphoscopic Traits. Paper at American Association of Physical Anthropology, at Austin, TX. Potts R. 1998. Variability selection in hominid evolution. Evolutionary Anthropology 7 (3):81- 96. Potts R. 2013. Hominin evolution in settings of strong environmental variability. Quaternary Science Reviews 73:1-13. Potts R, and JT Faith. 2015. Alternating high and low climate variability: The context of natural selection and speciation in Plio-Pleistocene hominin evolution. Journal of Human Evolution 87:5-20. Raff JA. 2019. Ancient DNA and Bioarchaeology. In A Companion to Anthropological Genetics, edited by O. R. DH. Hoboken, New Jersey: John Wiley & Sons, Inc. Ramachandran S, O Deshpande, CC Roseman, NA Rosenberg, MW Feldman, and LL Cavalli- Sforza. 2005. Support from the relationship of genetic and geographic distance in human 109 populations for a serial founder effect originating in Africa. Proceedings of the National Academy of Science of the US 102 (44):15942-7. Rasmussen M, M Sikora, A Albrechtsen, TS Korneliussen, JV Moreno-Mayar, GD Poznik, CPE Zollikofer, MS Ponce de León, ME Allentoft, I Moltke, H Jónsson, C Valdiosera, RS Malhi, L Orlando, CD Bustamante, TW Stafford, DJ Meltzer, R Nielsen, and E Willerslev. 2015. The ancestry and affiliations of Kennewick Man. Nature 523 (7561):455-458. Rathmann H, and H Reyes-Centeno. 2020. Testing the utility of dental morphological trait combinations for inferring human neutral genetic variation. Proceedings of the National Academy of Sciences of the United States 117 (20):10769. Rathmann H, H Reyes-Centeno, S Ghirotto, N Creanza, T Hanihara, and K Harvati. 2017. Reconstructing human population history from dental phenotypes. Scientific Reports 7 (1):12495. Ratliff MD. 2014. Evaluating Morphoscopic Trait Frequencies of Southeast Asians and Pacific Islanders. Missoula, MT University of Montana 2014-01-01T08:00:00Z. Relethford JH. 1994. Craniometric variation among modern human populations. American Journal of Physical Anthropology 95 (1):53. Relethford JH. 2002. Apportionment of global human genetic diversity based on craniometrics and skin color. American Journal of Physical Anthropology 118 (4):393-398. Relethford JH. 2004a. Boas and Beyond: Migration and Craniometric Variation. American Journal of Human Biology 16 (4):379-386. Relethford JH. 2004b. Global patterns of isolation by distance based on genetic and morphological data. Human Biology 76:499-513. Relethford JH. 2009. Race and global patterns of phenotypic variation. American Journal of Physical Anthropology 139 (1):16. Relethford JH. 2019. Human Population Structure and History. In A Companion to Anthropological Genetics, edited by DH O’Rourke. Hoboken, New Jersey: John Wiley & Sons, Inc. Relethford JH, and HC Harpending. 1994. Craniometric variation, genetic theory, and modern human origins. American Journal of Physical Anthropology 95 (3):249-70. Reyes-Centeno H, S Ghirotto, F Détroit, D Grimaud-Hervé, G Barbujani, and K Harvati. 2014. Genomic and cranial phenotype data support multiple modern human dispersals from Africa and a southern route into Asia. Proceedings of the National Academy of Sciences of the United States of America 111 (20):7248-7253. 110 Reyes-Centeno H, S Ghirotto, and K Harvati. 2017b. Genomic validation of the differential preservation of population history in modern human cranial anatomy. American Journal of Physical Anthropology 162 (1):170-179. Reyes-Centeno H, K Harvati, and G Jäger. 2016. Tracking modern human population history from linguistic and cranial phenotype. Scientific Reports 6 (1):36645. Reyes-Centeno H, and JT Hefner. 2021. Evolution of cranial macromorphscopic trait variation in modern humans. Annual Human Biology Association Conference. Virtual. Reyes-Centeno H, H Rathmann, T Hanihara, and K Harvati. 2017a. Testing Modern Human Out- of-Africa Dispersal Models Using Dental Nonmetric Data. Current Anthropology 58 (S17):S406-S417. Roberts DF. 1978. Climate and Human Variability. 2nd ed. Menlo Park, C: Benjamin- Cummings. Roberts DF. 1995. The pervasiveness of plasticity. In Human Variability and Plasticity, edited by B Bogin and CGN Mascie-Taylor. Cambridge: Cambridge University Press. Roberts DF. 2011. The pervasiveness of plasticity. In Human Variability and Plasticity, edited by CGN Mascie-Taylor and B Bogin. Roseman CC. 2004. Detecting interregionally diversifying natural selection on modern human cranial form by using matched molecular and morphometric data. National Academy of Sciences. Roseman CC, and TD Weaver. 2004. Multivariate apportionment of global human craniometric diversity. American Journal of Physical Anthropology 125 (3):257. Roseman CC, and TD Weaver. 2007. Molecules versus morphology? Not for the human cranium. Bioessays 29:1185-1188. Roseman CC, and BM Auerbach. 2015. Ecogeography, genetics, and the evolution of human body form. Journal of Human Evolution 78:80-90. Rosenberg NA, S Mahajan, C Gonzalez-Quevedo, MG Blum, L Nino-Rosales, V Ninis, P Das, M Hegde, L Molinari, G Zapata, JL Weber, JW Belmont, and PI Patel. 2006. Low levels of genetic divergence across geographically and linguistically diverse populations from India. PLoS Genetics 2 (12):e215. Rosenberg NA, S Mahajan, S Ramachandran, C Zhao, JK Pritchard, and MW Feldman. 2005. Clines, clusters, and the effect of study design on the inference of human population structure. PLoS Genetics 1 (6):e70. Rosenberg NA, JK Pritchard, JL Weber, HM Cann, KK Kidd, LA Zhivotovsky, and MW Feldman. 2002. Genetic structure of human populations. Science 298 (5602):2381-5. 111 Ross AH., and MA Pilloud. 2021. The need to incorporate human variation and evolutionary theory in forensic anthropology: A call for reform. American Journal of Physical Anthropology 176 (4):672-683. Rubicz RC, PE Melton, and MH Crawford. 2007. Molecular Markers in Anthropological Genetic Studies. In Anthropological Genetics: Theory, Methods, and Applications, edited by MH Crawford. Cambridge: Cambridge University Press. Ruff CB. 1994. Morphological adaptation to climate in modern and fossil hominids. Yearbook of Physical Anthropology 37:65. Ruff CB. 2002. Variation in Human Body Size and Shape. Annual Review of Anthropology 31:211-232. Salojärvi J. 2019. Computational Tools for Population Genomics. In Population Genomics: Concepts, Approaches, and Applications, edited by R. Op. Cham, Switzerland: Springer Nature Switzerland. Sanchez-Faddeev H, J Pijpe, D van Bodegom, T van der Hulle, KJ van der Gaag, UK Eriksson, T Spear, RGJ Westendorp, and P de Knijff. 2013. Ancestral Stories of Ghanaian Bimoba Reflect Millennia-Old Genetic Lineages. PLOS ONE 8 (6):e65690. Santos F. 2018. AnthropMMD: an R package with a graphical user interface for the mean measure of divergence. American Journal of Physical Anthropology 165 (1):200-205. Saunders SR. 1989. Nonmetric Skeletal Variation. In Reconstruction of Life from the Skeleton, edited by M Iscan and K Kennedy. New York: Liss, AR. Schell LM. 2011. Human plasticity: history and future research. In Human Variability and Plasticity, edited by CGN Mascie-Taylor and B Bogin. Boca Raton, FL: CRC Press Schlichting CD. 2021. Plasticity and Evolutionary Theory: Where We Are and Where We Should Be Going. In Phenotypic Plasticity & Evolution: Causes, Consequences, Controversies, edited by CGN Mascie-Taylor and B Bogin. Boca Raton, FL: CRC Press. Serre D, and S Pääbo. 2004. Evidence for gradients of human genetic diversity within and among continents. Genome Research 14 (9):1679-85. Shaffer JR, E Orlova, MK Lee, EJ Leslie, ZD Raffensperger, CL Heike, ML Cunningham, JT Hecht, CH Kau, NL Nidey, LM Moreno, GL Wehby, JC Murray, CA Laurie, CC Laurie, J Cole, T Ferrara, S Santorico, O Klein, W Mio, E Feingold, B Hallgrimsson, RA Spritz, ML Marazita, and SM Weinberg. 2016. Genome-Wide Association Study Reveals Multiple Loci Influencing Normal Human Facial Morphology. PLoS Genetics 12 (8):e1006149. Sjøvold T. 1973. The occurrence of minor non-metrical skeletal variants in the skeleton and their quantitative treatment for population comparisons. Homo 24:204-233. 112 Skoglund P, A Mallick, MC Bortolini, N Chennagiri, T Hünemeier, ML Petzl-Erler, FM Salzano, N Patterson, and D Reich. 2015. Genetic evidence for two founding populations of the Americas. Nature 525 (7567):104-108. Smith HF. 2009. Which cranial regions reflect molecular distances reliably in humans? Evidence from three-dimensional morphology. American Journal of Human Biology 21:36-47. Snoddy AME, HR Buckley, GE Elliott, VG Standen, BT Arriaza, and SE Halcrow. 2018. Macroscopic features of scurvy in human skeletal remains: A literature synthesis and diagnostic guide. American Journal of Physical Anthropology 167 (4):876-895. Sparks CS, and RL Jantz. 2002. A reassessment of human cranial plasticity: Boas revisited. Proceedings of the National Academy of Sciences of the United States of America 99 (23):14636-14639. Sparks CS, and RL Jantz. 2003. Changing Times, Changing Faces: Franz Boas's Immigrant Study in Modern Perspective. American Anthropologist 105 (2):333-337. Spradley MK. 2006. Biological anthropological aspects of the African diaspora: geographic origins, secular trends, and plastic versus genetic influences utilizing craniometric data: The University of Tennessee. Spradley MK. 2021. Use of craniometric data to facilitate migrant identifications at the United States/Mexico border. American Journal of Physical Anthropology 175 (2):486-496. Stoneking M. 2016. An Introduction to Molecular Anthropology. Hoboken, New Jersey: John Wiley & Sons, Inc. Stull KE, EJ Bartelink, AR Klales, GE Berg, MW Kenyhercz, EN L'Abbé, MC Go, K MCormick, and C Mariscal. 2021. Commentary on: JD Bethard, EA DiGangi. Letter to the Editor—Moving beyond a lost cause: Forensic anthropology and ancestry estimates in the United States. Journal of Forensic Sciences 66 (1):417-420. doi: 10.1111/1556‐4029.14513. Sutter R, and L Mertz. 2004. Nonmetric cranial trait variation and prehistoric biocultural change in Azapa Valley, Chile. American Journal of Physical Anthropology 123:130-145. Suzuki Y, KZ McKenna, and HF Nijhout. 2020. Regulation of phenotypic plasticity from the perspective of evolutionary developmental biology. In Phenotypic Switching, edited by H Levine, MK Jolly, PK and V Nanjundiah. Academic Press. Tagaya A. A general model for symmetry and asymmetry of nonmetric traits and congenital anomalies and tumors: reviving the proposals sacrificed to false myths. Anthropological Science. Taylor NAS. 2006a. Ethnic differences in thermoregulation: Genotypic versus phenotypic heat adaptation. Journal of Thermal Biology 31 (1):90-104. Taylor NAS. 2006b. Challenges to Temperature Regulation When Working in Hot Environments. Industrial Health 44 (3):331-344. 113 Templeton AR. 2019. Human Population Genetics and Genomics. London: Academic Press. Teplitsky C, and V Millien. 2014. Climate warming and Bergmann's rule through time: is there any evidence? Evolutionary Applications 7 (1):156-168. Teplitsky C, JA Mills, JS Alho, JW Yarrall, and J Merilä. 2008. Bergmann's Rule and Climate Change Revisited: Disentangling Environmental and Genetic Responses in a Wild Bird Population. Proceedings of the National Academy of Sciences of the United States of America 105 (36):13492-13496. Tishkoff SA, FA Reed, FR Friedlaender, C Ehret, A Ranciaro, A Froment, JB Hirbo, AA Awomoyi, JM Bodo, O Doumbo, M Ibrahim, AT Juma, MJ Kotze, G Lema, JH Moore, H Mortensen, TB Nyambo, SA Omar, K Powell, GS Pretorius, MW Smith, MA Thera, C Wambebe, JL Weber, and SM Williams. 2009. The genetic structure and history of Africans and African Americans. Science 324 (5930):1035-44. Trinkaus E. 1978. Bilateral asymmetry of human skeletal non-metric traits. American Journal of Physical Anthropology 49 (3):315-8. Trinkaus E. 1981. Neandertal limb proportions and cold adaptation. In Aspects of Human Evolution, edited by C Stringer. London: Taylor and Francis. Underhill PA, P Shen, AA Lin, L Jin, G Passarino, WH Yang, E Kauffman, B Bonné-Tamir, J Bertranpetit, P Francalacci, M Ibrahim, T Jenkins, JR Kidd, SQ Mehdi, MT Seielstad, RS Wells, A Piazza, RW Davis, MW Feldman, LL Cavalli-Sforza, and PJ Oefner. 2000. Y chromosome sequence variation and the history of human populations. National Genet 26 (3):358-61. Underhill PA, and T Kivisild. 2007. Use of Y Chromosome and Mitochondrial DNA Population Structure in Tracing Human Migrations. Annual Review of Genetics 41 (1):539-564. van Buuren S, and K Groothuis-Oudshoorn. 2011. Multivariate Imputation by Chained Equations. Journal of Statistical Software 45 (3):1-67. Vieira, MLC, L Santini, AL Diniz, and C de Freitas Munhoz. 2016. Microsatellite markers: what they mean and why they are so useful. Genetics and Molecular Biology 39 (3):312-328. von Cramon-Taubadel N. 2009. Congruence of individual cranial bone morphology and neutral molecular affinity patterns in modern humans. American Journal of Physical Anthropology 140:205-215. von Cramon-Taubadel N, and SJ Lycett. 2008. Brief communication: human cranial variation fits iterative founder effect model with African origin. American Journal of Physical Anthropology 136:108-113. von Cramon-Taubadel N, and TD Weaver. 2009. Insights from a quantitative genetic approach to human morphological evolution. Evolutionary Anthropology 18:237-240. 114 Walker JEC, RE Wells, and EW Merrill. 1961. Heat and water exchange in the respiratory tract. American Journal of Medicine 30:259-267. Wang S, CM Lewis, M Jakobsson, S Ramachandran, N Ray, G Bedoya, W Rojas, MV Parra, JA Molina, C Gallo, G Mazzotti, G Poletti, K Hill, AM Hurtado, D Labuda, W Klitz, R Barrantes, MC Bortolini, FM Salzano, ML Petzl-Erler, LT Tsuneto, E Llop, F Rothhammer, L Excoffier, MW Feldman, NA Rosenberg, and A Ruiz-Linares. 2007. Genetic variation and population structure in native Americans. PLoS Genetics 3 (11):e185. Wang S, N Ray, W Rojas, MV Parra, G Bedoya, C Gallo, G Poletti, G Mazzotti, K Hill, AM Hurtado, B Camrena, H Nicolini, W Klitz, R Barrantes, JA Molina, NB Freimer, MC Bortolini, FM Salzano, ML Petzl-Erler, LT Tsuneto, JE Dipierri, EL Alfaro, G Bailliet, NO Bianchi, E Llop, F Rothhammer, L Excoffier, and A Ruiz-Linares. 2008. Geographic patterns of genome admixture in Latin American Mestizos. PLoS Genetics 4 (3):e1000037. Weiss KM, Buchanon AV. 2004. Genetics and the Logic of Evolution. Hoboken, New Jersey: John Wiley & Sons, Inc. Weiner JS. 1954. Nose shape and climate. American Journal of Physical Anthropology 12 (4):615-618. Wells SR, and KM Dlugosch. 2019. Population Genomics of Colonization and Invasion. In Population Genomics: Concepts, Approaches, and Applications, edited by O Rojora. Cham, Switzerland: Springer Nature Switzerland. West-Eberhard MJ. 1989. Phenotypic Plasticity and the Origins of Diversity. Annual Review of Ecology and Systematics 20:249-278. West-Eberhard MJ. 2003. Developmental Plasticity and Evolution: Oxford University Press. West-Eberhard MJ. 2008. Phenotypic Plasticity. In Encyclopedia of Ecology, edited by SE Jørgensen and BD Fath. Oxford: Academic Press. West-Eberhard MJ. 2021. Forward: A perspective on 'plasticity'. In Phenotypic Plasticity & Evolution: Context, Causes, Consequences, edited by D Pfenning: CRC Press. Western AG, and JJ Bekvalac. 2017. Hyperostosis frontalis interna in female historic skeletal populations: Age, sex hormones and the impact of industrialization. American Journal of Physical Anthropology 162 (3):501-515. White MD, JG Greiner, and PLL McDonald. 2011. Point: humans do demonstrate selective brain cooling during hyperthermia. Journal of Applied Physiology 110 (2):569-571. Winburn AP, and B Algee‐Hewitt. 2021. Evaluating population affinity estimates in forensic anthropology: Insights from the forensic anthropology database for assessing methods accuracy (FADAMA). Journal of Forensic Sciences 66 (4):1210-1219. 115 Wood C. 2012. The Influence of Growth and Development in the Expression of Human Morphological Variation, University of Toronto, Toronto. Yokley TR. 2009. Ecogeographic variation in human nasal passages. American Journal of Physical Anthropology 138 (1):11. Yom-Tov Y, and E Geffen. 2011. Recent spatial and temporal changes in body size of terrestrial vertebrates: probable causes and pitfalls. Biological Reviews 86 (2):531-541. Zaidi AA, BC Mattern, P Claes, B McEvoy, C Hughes, and MD Shriver. 2018. Correction: Investigating the case of human nose shape and climate adaptation. PLoS Genetics 14 (1):e1007207-e1007207. Zhang ZQ, and JL Yang. 2015. Biomechanical Dynamics of Cranial Sutures during Simulated Impulsive Loading. Applied Bionics and Biomechanics 2015:59684. 116