Morphological analyses

We examined morphological relationships among body parts and facial features using published anthropometric data sets. We focused our analyses on the ANSUR II data set because it provides a large, consistent database of individual level facial and body measurements. We analyzed the linear anthropometric measurements (Supplementary Tables 1 and 2). In our analysis, we considered four groups of service members based on their sex and racial identity: African American females (n=181, mean height=64.29±0.17 inches), African American males (n=457, mean height=69.12±0.13 inches), European American females (n=204, mean height=64.27±0.15 inches), and European American males (n=1,168, mean height=69.32±0.08 inches). Compared with the general civilian population, the individuals measured in the ANSUR II data set tend to be taller and have lower levels of body fat55. Neither of these factors should influence our results or conclusions because our comparisons use facial and body measurements from the same individuals. Identity signalling predicts that traits used for recognition will have greater variance and be less correlated with each other compared with non-recognition traits in the same group of individuals.

Using the ANSUR II data set, we tested two predictions of the identity-signalling hypothesis. First, we considered the levels of variation in each trait by calculating the coefficient of variation—by dividing the standard deviation of each trait by the mean. Coefficients of variation provide a scale-free method for comparing variation across samples that differ in average size as is the case for human morphological data. Second, we considered the correlations among traits by calculating the inter-trait Pearson’s correlations for all pairwise combination of traits within each class of traits. To compare the distribution of correlation coefficients between bodies and faces, we recorded the correlation coefficients significant at the P<0.05 level. For any pair of traits which did not show a significant correlation at P<0.05, we recorded the correlation as 0. Pearon’s correlation test is sensitive to the sample size such that correlations are more likely to be significant when larger samples are used. Therefore, comparisons between the different groups considered should be made with caution because of differences in sample size. For example, differences in the percentage of significant pairwise comparisons between males and females likely reflect differences in samples sizes. Within a group, however, the same individuals were measured for both facial and body traits, providing a direct comparison of the relative degree of correlation among traits.

Selection of genomic regions for analysis

Face-associated SNPs were taken from two recent genome-wide association studies of normal facial morphology. Paternoster et al.33 conducted a discovery phase association study where they examined the relationship between facial characteristics and more than 2.5 million imputed SNPs in a sample of 2,185 15 year olds from the Avon Longitudinal Study of Parents and their Children56. Only subjects who genetically clustered with the CEU HapMap population were included in their analysis. The study identified 30 loci associated with facial morphology at P<5 × 10−7, which we examined in our study. Liu et al.32 examined the relationship between facial morphology and more than 2.5 million SNPs in a discovery phase association study of 5,388 adults. The samples in the Liu et al.32 study came individuals of European ancestry living in the Netherlands, Australia, Canada, Germany and the United Kingdom. They identified 29 loci associated with facial morphology at P<5 × 10−7, which we examined in this study. None of the SNPs identified by the two studies overlapped, providing a total of 59 loci for investigation. In both studies, multiple linked SNPs were often identified in association with a particular phenotype. When more than one SNP was associated with a trait, we chose the SNP with the smallest P-value within a 1-MB region of a chromosome from the association study. The SNPs identified for further examination from the two studies include one from each of 59 loci distributed throughout the autosomes. The SNPs are largely intergenic (95%) though a few occur within introns (5%). None were located in coding regions.

We compared face-associated genomic regions with two sets of control regions. First, we examined SNPs associated with height taken from the genome-wide association study catalog of the National Human Genome Research Institute ( www.genome.gov) on 25 April 2013. In order to prevent multiple sampling of any regions, we only considered SNPs that were separated by more than 4 kb. In the instances where multiple nearby SNPs had been associated with height, we chose the SNP that had been associated with the smallest P-value as reported in the genome-wide association study catalog. We excluded six SNPs associated with height that fall within the HLA region, though including the SNPs in our analyses does not alter our pattern of results. This produced a total of 365 loci associated with variation in height. Like faces, height is a composite character that depends on the morphology of numerous different bones. Additionally, both height and facial morphology have complex genetic bases with many loci of small effect contributing overall phenotypic variation43. SNPs associated with height are predominantly located in intronic regions (54%) and intergenic regions (36%) with a smaller percent found near the 3′ and 5′ end of genes (6%) or exons (3%). Second, we considered the genome-wide patterns of diversity by examining 5,000 2 kb intergenic regions. We identified putatively neutral intergenic regions in Europeans using the Neutral Region Explorer webserver42. The same set of intergenic regions was used for all populations.

We analyzed regions surrounding the SNPs identified by the association studies at a set window size. The causative mutations underlying the traits are not known, though are likely to be located near the SNPs identified through genome-wide association studies57. The a priori best choice for a window size is not clear, though the patterns of elevated nucleotide diversity we observer are seen over a range of window sizes (Supplementary Fig. 1). We chose to analyze 1 kb both up and downstream of the SNPs, providing windows of 2 kb. Smaller window sizes show marked increased variance in the summary statistics across loci (Supplementary Fig. 1), though this variance levels off at window sizes of 2 kb or greater.

Summary statistics

We calculated summary statistics for each population using binary SNP and indel data from 1000 Genome Project Phase 1 variants. Nine non-admixed populations originating from Europe (CEU: Utah residents with Northern and Western European Ancestry; GBR: British from England and Scotland; FIN: Finnish from Finland; TSI: Toscani from Italia), Asia (CHB: Han Chinese in Beijing, China, CHS: Southern Han Chinese, JPT: Japanese from Tokyo, Japan) and Africa (LWK: Luhya from Webuye, Kenya; YRI: Yoruba from Ibadan, Nigeria) were considered in our study. We downloaded the population data to Galaxy58 using the Table Browser function of the UCSC genome browser. We filtered the data based on the sets of 2-kb windows for face, height and intergenic regions to produce three files for each population, which we subsequently examined using custom macros in Excel.

We used folded site frequency spectra to examine the distribution of minor allele frequencies among SNPs found within each of the demarcated regions. The expected distribution of allele frequencies at loci underlying a polygenic trait under negative frequency-dependent selection is unclear and will depend on the exact form of selection and the genetic architecture of the trait59. Nonetheless, frequency-dependent selection is expected to maintain alleles in a population, on average, longer than expected for neutral alleles60. Thus, the distribution of allele frequencies should differ from that expected in a stationary population at mutation drift equilibrium. In particular, we expect fewer rare alleles under a scenario of frequency-dependent selection. Spectra were compared using the raw counts of SNPs with each minor allele frequency using a MWU-test. To graph the folded site frequency spectra, we binned data into ranges of minor allele frequencies.

In addition to the aggregated site frequency spectrum analysis, we also considered the distribution of multiple summary statistics of genetic diversity across the loci considered within our study. We calculated the following summary statistics for each 2-kb window: π, π corrected for human–macaque divergence, Watterson’s θ, Tajima’s D and Fu and Li’s D*. Both π and θ are estimators of the neutral mutation parameter, 4 N e μ. π is based on the number of pairwise differences among sequences within a sample and θ is based on the proportion of segregating sites. Loci under frequency-dependent selection are expected to show elevated values for π because frequency-dependent selection maintains alleles over longer periods of time. Older alleles accumulate mutations and therefore show higher levels of pairwise sequence divergence. We also examined the distribution of π corrected for the rate of divergence between humans and macaques. Different regions of the genome are known to experience differences in rates of mutation61. Loci with higher mutation rates will show elevated levels of π. The rate of divergence between humans and macaques provides a means of estimating the relative differences in mutation rates among loci62. The maintenance of multiple alleles in a population under frequency-dependent selection is also expected to lead to higher estimates of θ. Tajima’s D is the normalized difference between π and θ. Tajima’s D takes on positive values when there is an excess of intermediate-frequency variants and negative values when there is an excess of rare variants. Fu and Li’s D* is based on the number of nucleotide variants observed only once in a sample63. Negative measures of Fu and Li’s D* indicate an excess of singletons. Loci under frequency-dependent selection are expected to have a relatively smaller number of singletons and therefore more positive values of Fu and Li’s D*.

We calculated the summary statistics using the allele frequencies given in the phase 1 variant files from the 1000 Genomes project. The short indels recorded in the data set were considered in the same manner as SNPs. Human–macaque divergence data were estimated using the LastZ alignment of the two reference genomes. Only regions with alignments between the two species’ genomes were considered in the analysis of π corrected for divergence with macaques (faces=58 regions, height=356 regions, neutral=4,873 regions) and for subsequent analyses using this statistic. For the aligned regions, the average alignment lengths were 1,832.21±4.49 sites out of 2,000. We compared the distribution of summary statistics for face-associated loci with the distributions for the two control data sets using one-tailed MWU-tests.

Patterns of diversity across populations

We asked whether the same loci showed elevated diversity in different populations. To do this, we identified loci in each population for which π/divergence and Tajima’s D were above the 95th percentile as determined from the empirical distribution of intergenic regions examined within that population. We then asked whether or not a disproportionate number of loci with elevated diversity was shared between continents for face regions compared with the intergenic regions examined. For the nine loci showing elevated diversity in at least one population, we investigated the patterns of haplotype sharing across populations. We examined the sequences in the 2-kb window used for previous analyses. For those coordinates, we downloaded a combined PED file including CHB, FIN and YRI from the 1000 Genomes project site (browser.1000genomes.org). We converted the PED files to fasta format using PGD Spider64. This procedure produced a fasta file containing the polymorphic sites found within the examined loci. Using the ‘pegas’ package in R65, we created haplotype networks for each of the loci.

Sliding window analyses

To examine the extent to which selection for identity signalling has been shared or divergent across continents, we conducted a sliding window analysis of the regions identified as having elevated diversity in at least one population. We calculated π and Fst for 5-kb windows every 1 kb for a total of 200 kb. π was calculated for one representative population for each continent (FIN, CHB and YRI). We estimated levels of differentiation between populations using Hudson’s Fst following ref. 66 as it produces unbiased estimates of Fst and is less sensitive to sample size and rare variants than other estimates of Fst such as Weir and Hill67. We estimated Fst for each set of SNPs considered by calculating a ratio of averages rather than an average of ratios, as the former is less sensitive to the presence of rare variants in a sample66.

The SNPs identified in association with facial morphology are not found in coding sequences, so they are likely to influence gene regulation or splicing in some manner. For the two loci examined in greater detail, we used the UCSC genome browser to identify polymorphic sites in ENCODE regulatory regions. We focused on three ENCODE tracks in the UCSC browser68: H3K27Ac marks, DNase sensitivity clusters and transcription factor-binding sites. The H3K27Ac marks show regions for which there is CHIP-seq based evidence of enrichment for the H3K27Ac histone mark. H3K27 acetylation is associated with enhanced transcription. DNase sensitivity clusters show regions sensitive to DNase as assessed across 125 cell types. Promoters and other regulatory regions tend to be DNase sensitive. The transcription factor track shows regions with evidence of transcription factor-binding sites.

Gene trees

For the two 5-kb loci examined, we constructed maximum likelihood gene trees with 10 sequences each from the FIN, CHB and YRI 1000 Genomes populations for a total of thirty sequences. Additionally, we included the human and chimpanzee reference sequences as well as sequences for Denisovans69 and Neanderthals ( http://cdna.eva.mpg.de/neandertal/altai/AltaiNeandertal/bam/). We downloaded the alignment of the human and chimpanzee reference sequences from Ensembl. Denisovan sequences were downloaded using the Table Browser function of the UCSC Genome Browser. The draft Altai Neaderthal sequences were downloaded for the relevant chromosomes from the Department of Evolutionary Genetics at the Max Planck Institute’s website. We constructed individual sequences for the 1000 Genomes, Denisovan and Neanderthal by manually altering the human reference sequence in accordance with the data found in the respective VCF files using Mega 5.2.1 (ref. 70)70. For the phased 1000 Genomes data, we selected one chromosome per individual sample. For the Neanderthal and Denisovan sequences, we included all of the sites that differed from the human reference to make a single sequence. After removing sites with gaps in the alignment, we constructed a maximum likelihood tree using a general time reversible model with a gamma distribution of invariant sites.