D1080–D1089 Nucleic Acids Research, 2019, Vol. 47, Database issue Published online 18 October 2018 doi: 10.1093/nar/gky959 PopHumanScan: the online catalog of human genome adaptation Jesús Murga-Moreno † , Marta Coronado-Zamora † , Alejandra Bodelón, Antonio Barbadilla * and Sònia Casillas * Institut de Biotecnologia i de Biomedicina and Departament de Genètica i de Microbiologia, Universitat Autònoma de Barcelona, 08193 Bellaterra, Barcelona, Spain Downloaded from https://academic.oup.com/nar/article-abstract/47/D1/D1080/5134333 by guest on 05 June 2020 Received August 15, 2018; Revised September 21, 2018; Editorial Decision October 03, 2018; Accepted October 04, 2018 ABSTRACT cas and Eurasia during the last ice age, limiting exposure to sunlight as we moved to higher latitudes or contact with Since the migrations that led humans to colonize new pathogens. Part of the incorporated genetic innovations Earth, our species has faced frequent adaptive chal- may have been introgressed from archaic hominins that left lenges that have left signatures in the landscape of Africa before us, including Neanderthals and Denisovans, genetic variation and that we can identify in our to- with whom we encountered and interbred before they got day’s genomes. Here, we (i) perform an outlier ap- extinct. Around 1–6% of any modern non-African human proach on eight different population genetic statis- genome can be traced back to the genomes of these archaic tics for 22 non-admixed human populations of the populations (2). Another dramatic change occurred within Phase III of the 1000 Genomes Project to detect se- the past 10 000 years coinciding with the transition from a lective sweeps at different historical ages, as well hunting-gathering lifestyle to farming. Selection pressures as events of recurrent positive selection in the hu- for adapting to large settlements and new diets favored ge- netic variants associated with innate immune response, fatty man lineage; and (ii) create PopHumanScan, an on- acid metabolic efficiency, and lactose tolerance, among oth- line catalog that compiles and annotates all candi- ers (3). date regions under selection to facilitate their valida- These selection pressures left signatures in the landscape tion and thoroughly analysis. Well-known examples of genetic variation that can be identified in our today’s of human genetic adaptation published elsewhere genomes (4). Starting from single-locus studies to the first are included in the catalog, as well as hundreds of large-scale catalogs of genetic variation (5–8), dozens of other attractive candidates that will require further targets of positive selection have been identified, providing investigation. Designed as a collaborative database, important insights into recent human evolutionary history PopHumanScan aims to become a central reposi- (3,9,10). Even though genome-wide HapMap genotyping tory to share information, guide future studies and data is able to disentangle the effects of demography and help advance our understanding of how selection selection better than single-locus approaches, it still has the problem of ascertainment bias, which may alter the site fre- has modeled our genomes as a response to changes quency spectrum (SFS) of analyzed single nucleotide poly- in the environment or lifestyle of human popula- morphisms (SNPs) (11). The availability of the most com- tions. PopHumanScan is open and freely available prehensive worldwide nucleotide variation dataset so far at https://pophumanscan.uab.cat. from the 1000 Genomes Project (1000GP) (12,13), based on whole-genome re-sequencing, provides the human lineage with an abundant, ascertained variation dataset on which INTRODUCTION to test molecular population genetics hypotheses and even- Since the split with chimpanzees, and especially since the tually pinpoint targets of positive selection in one or more migrations that led humans to colonize almost every sin- human populations that escape from the background evo- gle place on Earth, our species has faced frequent environ- lutionary dynamics of genetic variation (14). mental and social changes that have shaped the variation To gain deeper understanding of how environmental and patterns of our genomes through the action of natural se- social challenges have shaped our genomes through the ac- lection (1). These environmental challenges include, for ex- tion of natural selection, here we (i) perform a genome- ample, extreme cold temperatures in much of the Ameri- wide scan of selection on the latest version of the 1000GP * To whom correspondence should be addressed. Tel: +34 93 586 8958; Fax: +34 93 581 2011; Email: sonia.casillas@uab.cat Correspondence may also be addressed to Antonio Barbadilla. Tel: +34 93 586 8941; Fax: +34 93 581 2011; Email: antonio.barbadilla@uab.cat † The authors wish it to be known that, in their opinion, the first two authors should be regarded as joint First Authors.  C The Author(s) 2018. Published by Oxford University Press on behalf of Nucleic Acids Research. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted reuse, distribution, and reproduction in any medium, provided the original work is properly cited. Nucleic Acids Research, 2019, Vol. 47, Database issue D1081 data by surveying distinctive signatures of genomic varia- the 1000GP (13) (Supplementary Table S1), mapped to tion left by different selective events, and (ii) create an on- GRCh37/hg19. Specifically, values for seven different neu- line catalog of all candidate genomic regions under selec- trality tests have been obtained for each population in 186 tion to facilitate their validation and thorough analysis. As 549 10-kb non-overlapping sliding windows along the auto- far as we are concerned, dbPSHP (15) is the only previ- somes and the X chromosome (Figure 1). In addition, the ous online database that compiles putative positively se- McDonald and Kreitman test (MKT) (19), as well as the lected loci in human evolution. In their case, regions were proportion of substitutions that are adaptive (α) (20,21), extracted from curated publications based on genotyping were calculated on the protein-coding genes overlapping data––instead of whole-genome re-sequencing data––of the the candidate regions under selection identified with the HapMap III (8) and the 1000GP Pilot 1 (12), and the other seven statistics and that showed some variability in last update is reported as far as May 2014. On the other both polymorphism and divergence, according to gene an- hand, after the publication of the 1000 Genomes Selection notations from GENCODE release 27 (22) and PopHu- Downloaded from https://academic.oup.com/nar/article-abstract/47/D1/D1080/5134333 by guest on 05 June 2020 Browser 1.0 (16), Pybus et al. developed a machine-learning man polymorphism and divergence data (Figure 1). MKT- framework––Hierarchical Boosting––that combines the re- derived calculations were performed using the R pack- sults of multiple tests for detecting positive selection to clas- age iMKT (https://github.com/BGD-UAB/iMKT; last ac- sify genomic regions into different selection regimes (17). cessed: February 2018). In total, eight different neutrality They analyzed within-species polymorphism data for three tests were performed. They identify different types of sig- populations of the 1000GP Phase I (12), and the resulting natures that remain visible in the genomic sequences for a scores were made available as UCSC Hub Tracks. Here, we timescale ranging from few thousands of years after a single perform an outlier approach on the greatest number of sweep selection event to several million years for the case of population genetic statistics and sampled populations avail- recurrent selection (4,23). able so far. This genome-wide scan of selection is able to detect sweeps at different historical ages, as well as evi- Linkage Disequilibrium (LD). Selection signatures in the dence of recurrent selection in the human lineage since the LD, e.g. long haplotypes, were detected using two com- split between our species and chimpanzees. Results have plementary measures: iHS (24) and XP-EHH (9). iHS been made available in a collaborative, online database, Po- has good power to detect selective sweeps with haplotypes pHumanScan, which is aimed at compiling and annotat- at moderate frequency (50–80%), while XP-EHH is more ing adaptation events along the human evolutionary his- powerful for detecting selective sweeps when the selected tory. Well-known examples of human genetic adaptation haplotype has a frequency >80%. In the case of XP-EHH, published elsewhere are included in the catalog, as well as which analyzes pairs of populations, only pairs CEU-YRI, hundreds of other attractive candidates that will require CEU-CHB and YRI-CHB were considered, and the popu- a more thoroughly analysis. PopHumanScan graphically lation showing the evidence of selection was identified with represents each signature of selection within the empirical the locus-specific branch length method (LSBL) (25). These distributions of the corresponding DNA diversity statistic long-range haplotypes persist for relatively short periods of across populations. It also provides structural and func- type, typically up to ∼30 kya, and thus this signature allows tional annotations of the region, links to external databases us to identify recent selection events only (4). and cross-references to 268 publications. Site Frequency Spectrum (SFS). Five statistics have been POPHUMANSCAN ANALYSIS PIPELINE considered: four are based on both the allele frequency spec- trum and the levels of variability––Fay and Wu’s H (26), Fu We have designed and implemented a custom pipeline (Fig- and Li’s D and F (27), and Tajima’s D (28)––, and the other ure 1) to perform a genome-wide scan of selection. Specif- one is based on population differentiation––FST (29,30)––. ically, the pipeline processes eight different neutrality tests Fay and Wu’s H detects an excess of high-frequency derived calculated either in sliding windows along the genome or SNPs, compatible with an incomplete sweep or recombina- for each protein-coding gene, for 22 non-admixed human tion breaking swept linked SNPs. Fu and Li’s D and F and populations. The genomic regions identified should show Tajima’s D assess the lack or excess of rare alleles (or single- signatures that are compatible with natural selection hav- tons). Lack of rare alleles is compatible with balancing se- ing driven the evolution of the region at one or different lection, while an excess is normally explained by either posi- timescales, from recent selective sweeps to recurrent selec- tive or weakly deleterious selection. FST detects population- tion since the split between our species and chimpanzees. specific selective events that changed the genetic composi- These candidate regions under selection are further charac- tion of the affected population. It analyzes pairs of popula- terized with structural and functional annotations of that tions. The pairs CEU-YRI, CEU-CHB and YRI-CHB were particular region. Furthermore, 268 articles reporting ev- considered, and the population showing the evidence of se- idences of natural selection in genomic regions and genes lection was identified with the LSBL (25). SFS signatures using different statistical methods have been manually cu- can persist in the genomes for a longer period than LD, and rated and cross-referenced to the candidate regions detected thus selective events identified by shifts in the SFS might with our pipeline. have occurred up to ∼80 kya (4). Pre-processing of the PopHuman data Protein changes. Recurrent selection since the split be- Population genomic data were retrieved from PopHuman tween our species and chimpanzees (<6 mya) is detected (18) for 22 non-admixed populations of the Phase III of using a test based on comparisons of polymorphism and D1082 Nucleic Acids Research, 2019, Vol. 47, Database issue Downloaded from https://academic.oup.com/nar/article-abstract/47/D1/D1080/5134333 by guest on 05 June 2020 Figure 1. PopHumanScan pipeline. Starting from population genomic data retrieved from PopHuman, 8 different neutrality tests are analyzed in 22 non-admixed human populations (or 3 population pairs). Tests are color-coded depending on the type of signature they are able to detect: Linkage Disequilibrium (LD), Site Frequency Spectrum (SFS) and Protein Changes. The significance of each test is assessed either with a Fisher’s exact test or a rank score, for each of the 22 populations (or 3 population pairs) independently, and independently for autosomes and the X chromosome. Finally, candidate regions under selection are structurally and functionally annotated, and cross-referenced with 268 publications. Nucleic Acids Research, 2019, Vol. 47, Database issue D1083 divergence––MKT (19)––, and the result of the test is sum- they were completely located in DAC Blacklisted regions marized with the estimator α (20,21). For this calculation, (i.e. regions of the reference genome which are troublesome we used an MKT-based methodology that corrects for the for high throughput sequencing aligners) or partially over- presence of non-synonymous slightly deleterious segregat- lapped genomic gaps, as obtained from the UCSC (32). ing sites in order to avoid underestimating α (31). Therefore, a total of 2859 regions were finally considered. Genome-wide scan of selection Structural and functional annotations For the parameter α of MKT, evidence of positive selec- The final 2859 candidate regions under selection were struc- tion for protein-coding genes was inferred when α > 0 and turally and functionally characterized according to 15 dif- the Fisher’s Exact Test for the 2 × 2 MKT contingency ferent annotations categorized into five groups, extracted table was significant (P-value < 0.05) (Figure 1). Because from the UCSC (32) and two publicly available databases Downloaded from https://academic.oup.com/nar/article-abstract/47/D1/D1080/5134333 by guest on 05 June 2020 the other seven selection statistics have not been associated (33,34) (Figure 1). with a simple parametric distribution, candidate windows under selection were identified as the most extreme values Sequencing. (i) Mappability was assessed as the percent- (within the 0.05% tail) in the corresponding empirical dis- age of bases in the region that do not present any trouble- tribution. These empirical distributions were performed in- some to high-throughput sequencing aligners according to dependently for each of the 22 populations (or three pop- the DAC Blacklisted regions of the UCSC (35). (ii) Distance ulation pairs in the case of XP-EHH and FST ), and inde- to closest GAP was computed as the distance (in Mb) to the pendently for the autosomes and the X chromosome (to ac- closest gap (32). count for different demographic histories and the different effective population size of the autosomes compared to the Regulation. (iii) CpG Islands (36), (iv) Vista Enhancers X chromosome; chromosome Y was not analyzed). In to- (37), (v) Transcription Factor Binding Sites (TFBSs) tal, 116 empirical distributions were obtained for autoso- (32) and (vi) ORegAnno Regulatory Elements (38) were mal regions, and 91 for the X chromosome (data of iHS computed as the total number of elements contained in the and XP-EHH was not available for the X chromosome in region, as well as the number of SNPs contained in the over- PopHuman) (Figure 1). lapping elements. From the initial 186 549 10-kb non-overlapping windows from PopHuman for each population and statistic, those Comparative genomics. Evolutionary conservation of the containing <5 segregating sites were discarded (<0.2%) regions was assessed by considering the results of three (Supplementary Figure S1). Then, an empirical P-value was different algorithms––phastCons, PhyloP and GERP––on assigned to each of the remaining windows for each of the the multiple alignments of the genomes of 100 vertebrate 116 combinations of population (or population pair) and species (39). (vii) PhyloP Evolutionary Conservation and statistic, separately for the autosomes as a whole and the (viii) GERP Constrained Elements were assessed as the per- X chromosome. Specifically, for each window i in a pop- centage of bases in the region that have a score >2 for the ulation (or population pair), p is the quantile of that win- given statistic (i.e. constrained sites) (32). (ix) phastCons dow for statistic j, that is, its empirical P-value. In the case Evolutionary Conservation was calculated as the percentage of Tajima’s D, Fu and Li’s D and F, and Fay and Wu’s H, of bases that overlap phastCons conserved elements (32). two-tailed P-values were calculated. Once the significance for each individual 10-kb window in the genome was as- Structural variation. (x) InvFEST Inversions (34), (xi) sessed, a candidate region under selection was defined as DGV Structural Variants (40), (xii) RepeatMasker (41), (xiii) being a contiguous genomic region containing at least one Segmental Duplications (42) and (xiv) TRF Simple Tandem 10-kb significant window (P-value < 0.0005) and spanning Repeats (43) were assessed as the percentage of bases in the adjacent windows with P-values < 0.005 (Supplementary region that overlap these genomic elements. Figure S2). In addition, this region may span stretches <20- kb of contiguous nucleotides not analyzed in PopHuman Archaic introgression. (xv) Archaic introgression was as- (i.e. because they contain non-accessible bases according to sessed as the percentage of bases in the region that over- the Pilot-style Accessibility Mask of the 1000GP (13,18)). lap either Neanderthal or Denisova introgressed haplotypes This outlier approach was designed to face the unique fea- (33). tures and limitations of our PopHuman source data and to be highly conservative defining candidate regions under se- Published references lection. We expect that it likely results in an enriched set of genomic regions that have been targets of natural selection A total of 268 publications from 1954 to 2018 report- along the human evolutionary history (11), and refer to the ing either specific loci or multiple regions from a genome- outlier regions as candidate regions showing signatures of wide scan of selection in the human genome were cross- selection. referenced with our final 2859 candidate regions under se- Once candidate selected regions (or genes) were assigned lection (Figure 1 and Supplementary Table S2). Of these, for the 22 populations (or three population pairs) and eight 132 publications were directly extracted from the dbPSHP statistics, they were collapsed according to their coordinates database (15), while the other 136 were manually curated into a joint set of 2879 candidate regions under selection here. Exhaustive information from the main text and/or genome-wide. Of these, 20 regions were removed because supplementary figures and tables was extracted for each D1084 Nucleic Acids Research, 2019, Vol. 47, Database issue Downloaded from https://academic.oup.com/nar/article-abstract/47/D1/D1080/5134333 by guest on 05 June 2020 Figure 2. Simplified representation of the PopHumanScan interface. The main PopHumanScan table is displayed to the left, while the complete report for a particular candidate region under selection is displayed to the right. reported loci, including the genomic coordinates, affected (v) color-coded dots depicting the meta-population(s) that population(s), statistic(s), type of selection and PubMed ID. show signatures of selection (i.e. Europe (EUR), Africa Genomic coordinates were lifted over to GRCh37/hg19 us- (AFR), South-Asia (SAS) and/or East-Asia (EAS)) ing the LiftOver tool of the UCSC (32), or deduced from and (vi) the source that contributed the candidate region un- protein-coding gene location, if necessary. der selection. At the time of writing, all 2859 regions came uniquely from our genome-wide selection scan (i.e. source OVERVIEW OF THE POPHUMANSCAN ONLINE CAT- labeled as PopHumanScan), but additional data sources by ALOG contributors from the scientific community are expected once PopHumanScan is published (see next section). By In addition to the exhaustive genome-wide selection scan clicking the ⊕ icon at the beginning of each row, detailed that has been performed, we have also created PopHuman- information of the particular candidate region under se- Scan, a collaborative, online database that is aimed at com- lection is displayed, including the values for all significant piling and annotating adaptation events along the human statistics in all target populations (or population pairs) and evolutionary history (Figure 2). PopHumanScan reports an overview of the main structural and functional annota- each evidence of selection with the empirical distributions tions and cross-referenced publications (i.e. non-gray but- of the corresponding DNA diversity statistic across the hu- tons represent overlapping annotations or cross-referenced man genome and among populations, structural and func- publications), as well as access to the complete report for tional annotations of the region, links to external databases, the corresponding candidate region under selection. At last, as well as cross-references to 268 publications. several filters are available at the bottom of the page to nar- row the search. Implementation Complete report. A complete report for each candidate re- PopHumanScan is currently running under Apache on a gion can be accessed from the main table (Figure 2, right). CentOS 7.2 Linux x64 server with 16 Intel Xeon 2.4 GHz The first section of the report displays all the structural and processors and 32 GB RAM. It is mainly built on PHP functional annotations of the region, together with links to as backend framework. It also includes AJAX for specific external databases: (i) PopHuman (18), which complements file requests and MySQL for data storage. The client-side is the population genomics information; (ii) HaploReg (44), build on JavaScript and uses several JavaScript libraries, in- which allows the exploration of evolutionary conservation, cluding jQuery, the jQuery plugin DataTables and Plotly.js, expression eQTSs, epigenomic data and regulatory annota- as well as a custom Bootstrap 4 framework. tions; and (iii) Ensembl (45), which allows the exploration of the LD of the region, among others. The second section lists The PopHumanScan catalog all the genes contained in or partially overlapping the region (if any). For each encoded gene, a short description of the Main table. All 2859 candidate regions under selection are gene and associated Gene Ontology terms for the Biologi- displayed as rows in an interactive table (Figure 2, left). cal Process classification (46) are provided, along with links The information displayed in each row includes: (i) the ge- to external databases: Ensembl (45), NCBI (47), Uniprot nomic coordinates of the candidate locus, (ii) genes con- (48), UCSC (32), Expression Atlas (49), OMIM (50), Open tained in or partially overlapping the region (if any), (iii) Targets (51) and HumanMine (52). The third section con- the most extreme value for each of the eight statistics con- tains cross-referenced publications that support the selec- sidered (i.e. most extreme value in any 10-kb window in- tion evidence found in the region. The fourth section con- cluded in the region, for any population or population pair; tains an interactive graph showing recombination rate val- green (positive) and red (negative) values are outliers (P < ues in cM/Mb along the chromosome in which the region is 0.0005) in the corresponding empirical distribution), (iv) located, calculated from the recombination map by Bhérer color-coded dots depicting different types of selection sig- et al. (53) and extracted from PopHuman (18). The specific natures (i.e. LD, SFS and/or Protein Changes), location of the candidate region under selection is indicated Nucleic Acids Research, 2019, Vol. 47, Database issue D1085 Downloaded from https://academic.oup.com/nar/article-abstract/47/D1/D1080/5134333 by guest on 05 June 2020 Figure 3. Summary of the contents of PopHumanScan. (A) Number of candidate regions under selection unique and shared among the four meta- populations: Europe (EUR), Africa (AFR), Sourth-Asia (SAS) and East-Asia (EAS). (B) Number of candidate regions under selection unique and shared among the three different signature types: Linkage Disequilibrium (LD), Site Frequency Spectrum (SFS) and Protein Changes. (C) Number of candidate regions under selection overlapping different structural and functional annotations. (D) Number of candidate regions under selection cross-referenced with 0, 1, 2, 3, 4 or ≥5 published papers. with dashed vertical lines, and the solid horizontal line rep- well as instructions on how to contribute to it. Interestingly, resents the average recombination rate value in the candi- it also contains a complete tutorial introducing to the usage date region. At last, in the fifth section boxplots show the of the database through a step-by-step worked example. distribution of each significant statistic in all the popula- tions (or population pairs). Highlighted values correspond CONTENTS OF POPHUMANSCAN to those in the candidate region, and those in red are out- liers of the empirical distribution (P < 0.0005). At the time of writing, the PopHumanScan database con- tains 2859 candidate regions under selection derived from the genome-wide selection scan pipeline presented here. Re- Utilities and support resources gions are distributed homogeneously along the autosomes Contributing to PopHumanScan. PopHumanScan has and the X chromosome (Table 1 and Supplementary Fig- been devised as a collaborative database. In order to in- ure S3). Of these, 1453 regions (50.8%) overlap GENCODE corporate information contributed by the scientific commu- protein-coding genes, and 1986 regions (69.5%) are cross- nity, two password-protected tools have been implemented. referenced with at least one publication (Table 1 and Figure The first one allows users to add additional candidate re- 3). gions under selection in the catalog. All contributed regions will be subjected to manual curation and clearly labeled Selection signatures in meta-populations with a data source tag. The second tool allows manually cross-referencing candidate regions already present in the The total number of candidate regions showing signatures database. of selection in the four meta-populations is: 831 (29.1%) in EUR, of which 413 (49.7%) overlap protein-coding genes; Help and tutorial. This section documents the data used 1090 (38.1%) in AFR, of which 580 (53.2%) overlap protein- and the procedures implemented in PopHumanScan, as coding genes; 791 (27.7%) in SAS, of which 401 (50.7%) D1086 Nucleic Acids Research, 2019, Vol. 47, Database issue Table 1. Summary of the candidate regions under selection included in PopHumanScan Regions with selection signatures in meta-populations Regions with different types of signatures Regions overlap- Linkage Site ping Number of South- disequilib- Frequency protein- Regions cross- candidate European African Asian East-Asian rium (LD) Spectrum Protein coding referenced with Chromosome regions (EUR) (AFR) (SAS) (EAS) (SFS) Changes genes publications 1 214 57 84 48 66 46 176 2 123 152 2 253 77 95 60 81 56 203 2 131 191 3 201 61 81 49 54 34 173 0 116 153 4 241 62 110 68 59 42 207 1 111 173 Downloaded from https://academic.oup.com/nar/article-abstract/47/D1/D1080/5134333 by guest on 05 June 2020 5 166 46 62 41 52 32 140 1 85 119 6 171 42 83 42 49 29 146 1 81 118 7 144 53 58 49 49 25 125 0 74 115 8 164 41 57 51 48 32 135 0 68 108 9 96 28 33 22 34 12 85 1 51 59 10 141 45 47 45 48 21 126 1 70 103 11 120 35 37 34 44 22 99 2 62 86 12 142 44 52 35 42 24 123 1 87 105 13 82 17 31 27 21 9 75 0 23 51 14 74 19 29 25 22 15 61 0 39 57 15 79 34 29 17 25 12 70 1 47 57 16 88 20 20 41 30 21 68 1 49 69 17 86 30 32 34 24 18 71 1 64 69 18 56 15 23 14 18 7 49 0 24 43 19 67 21 29 21 27 5 62 4 42 38 20 59 19 22 14 23 10 52 1 29 49 21 34 13 11 10 16 7 28 0 11 26 22 39 8 16 8 15 8 35 1 27 32 X 142 44 49 36 37 ND 142 0 39 13 TOTAL 2859 831 1090 791 884 487 2451 21 1453 1986 29.1% 38.1% 27.7% 30.9% 17.0% 85.7% 0.7% 50.8% 69.4% ND = Not Determined overlap protein-coding genes; and 884 (30.9%) in EAS, of show signatures in either LD and/or SFS, as only genes which 424 (48.0%) overlap protein-coding genes (Table 1). overlapping candidate regions under selection detected by Most of the regions (82.5%) show signatures that are unique LD and/or SFS were tested for α (MKT). These results to one single meta-population (Figure 3A): 492 (17.2%) would indicate that the statistics we used in our genome- show signatures that are unique in EUR, 835 (29.2%) are wide scan of selection look at different characteristics of the unique in AFR, 433 (15.1%) are unique in SAS and 603 genetic variability of the region, and that they are largely (21.1%) are unique in EAS. Of the 1090 regions showing complementary. signatures in AFR, 76.6% are unique to AFR; while a lesser percentage––59.2, 54.7 and 68.2%––of the regions showing signatures in EUR, SAS and EAS, respectively, are unique Structural description of the regions to their meta-population. About one third (29.0%) of the candidate regions under selection are shared across popula- Region length. Most of the candidate regions under selec- tions within the same meta-population. This percentage is tion (63.6%) span one single 10-kb window, and the variable higher for candidate regions showing both LD and SFS sig- lengths of candidate regions follows a reversed J-shaped dis- natures (52.7%), it is 33.6% for candidate regions showing tribution (Supplementary Figure S4A). LD signatures only and 27.1% for candidate regions show- ing SFS signatures only. Distance between consecutive regions. The average dis- tance between consecutive candidate regions is ∼1 Mb, and the distribution of distances is also reversed J-shaped (Sup- Types of selection signature plementary Figure S4B). The total number of candidate regions showing distinct types of signatures of selection is: 487 (17.0%) for LD; 2451 Recombination. The average recombination rate of the (85.7%) for SFS; and 21 (0.7%) for Protein Changes (i.e. re- candidate regions is 0.71 cM/Mb, and the distribution of current selection since the split between humans and chim- recombination rates is again reversed J-shaped (Supplemen- panzees) (Table 1). Most of the regions (96.6%) show one tary Figure S4C). There is a strong, negative, non-linear single signature of selection (Figure 3B): 403 (14.1%) show association between recombination rate and both region LD signatures only; and 2358 (82.5%) show SFS signatures length (Supplementary Figure S5A) and distance between only. All genes showing evidence of recurrent selection also consecutive candidate regions (Supplementary Figure S5B). Nucleic Acids Research, 2019, Vol. 47, Database issue D1087 Functional description of the regions POPHUMANSCAN WITH AN EXAMPLE: SELECTION AT THE LACTASE LOCUS Regulation. Most of the candidate regions (90.5%) con- tain at least one regulatory element annotated in the ORe- The introduction of agriculture and cattle domestication in gAnno database, and 80.6% contain TFBSs (Figure 3C). On the Middle East and North Africa ∼10 000 years ago lead the contrary, VISTA enhancers are much less abundant in to strong selection pressure for the ability to digest milk as the genome and they are only found in 23 of the 2859 can- adults. This is accomplished if the enzyme lactase that me- didate regions (0.8%). CpG Islands are also in shortage and tabolizes lactose, encoded by the LCT gene, maintains high they are present in 9.3% of the regions. levels into adulthood, a characteristic that is called lactase persistence. Several variants near the LCT locus show some of the strongest signals of selection in the human genome for Comparative genomics. Nearly all (96.8%) candidate re- those populations that have traditionally practiced dairy- gions overlap phastCons conserved elements. In the case of Downloaded from https://academic.oup.com/nar/article-abstract/47/D1/D1080/5134333 by guest on 05 June 2020 ing, including a genetic variant in an intron of the gene GERP and PhyloP, 88.2 and 74.7% of the regions, respec- MCM6, upstream of LCT (3). tively, overlap constrained bases with score >2. The LCT locus is found inside the longest candidate re- gion under selection reported in PopHumanScan (∼1 Mb). Structural variation. The Database of Genomic Variants The region is located in the long arm of chromosome 2 (DGV) (40) is a very exhaustive database of structural vari- and contains 8 GENCODE protein-coding genes, includ- ants annotated in the human genome. One or more elements ing LCT and MCM6 (Figure 4). Our genome-wide scan of annotated in this database are present in 92.1% of the can- selection has detected signatures at four different statistics didate regions under selection. On the contrary, only 104 re- that span the three types of signatures: LD (iHS and XP- gions (3.6%) overlap validated polymorphic inversions from EHH), SFS (Fu and Li’s D) and Protein Changes (α). LD the manually curated InvFEST database (34). signatures involve basically EUR and AFR populations, while the signature of recurrent selection is more general to the four meta-populations. The region contains thousands Archaic introgression. A total of 1526 of the candidate re- of TFBSs and hundreds of ORegAnno regulatory elements, gions (53.4%) overlap haplotypes introgressed from either it overlaps evolutionary constrained elements, and >95% neanderthals or denisovans. This percentage is expected, as of the region overlaps haplotypes introgressed from Nean- introgressed haplotypes persisting in different present-day derthals. It has been reported in 24 published articles (of the human individuals cover 46.7% of the reference genome set of 268 that we considered). (33). CONCLUSION Cross-references with publications. A percentage of 69.5% of the candidate regions are cross-referenced with at least In summary, our exhaustive approach combining eight dif- one publication, and 36.0% are cross-referenced more than ferent statistics to detect candidate regions under selection once (Figure 3D). in 22 non-admixed human populations has been able to locate distinct signatures in 2859 regions that stand out from the background genomic variability, including abnor- mally long haplotypes, shifts in the SFS or excess of non- Gene ontology analysis synonymous substitutions between our species and chim- Our candidate regions overlap a total of 1447 unique GEN- panzees. Many of these regions probably manifest the foot- CODE protein-coding genes. These were functionally clas- prints of selective sweeps that occurred at different histor- sified into Gene Ontology (GO) terms (46) according to the ical ages, or recurrent selection that has been taking place PANTHER GO-Slim annotation dataset using the PAN- during the last millions of years. The PopHumanScan on- THER Classification System (54) (Supplementary Figure line database is going to facilitate the thorough analysis of S6). In addition, statistically over- and under-represented candidate regions under selection in the human genome by functions were analyzed using the complete GO annota- putting together all these evidences of selection with struc- tion dataset (46) using the same tool (Supplementary Ta- tural and functional annotations of the regions and cross- bles S3–5). Interestingly, among all Biological Process cate- references to previously published articles. Furthermore, gories, regulation of neuron projection development is over- the database can incorporate new data from the scientific represented (fold enrichment 1.88, False Discovery Rate community through specific build-in utilities. All in all, Po- (FDR) 1.23E-02), in addition to cellular component organi- pHumanScan aims to become a central repository to share zation (fold enrichment 1.24, FDR 1.59E-03) (Supplemen- information, guide future studies and contribute to the re- tary Table S4). Finally, several Cellular Component cate- search on human genome adaptation. gories are statistically over-represented, including presynap- tic membrane (fold enrichment 2.72, FDR 4.32E-02) (Sup- DATA AVAILABILITY plementary Table S5). In spite of finding some statistically over-represented GO categories in our genes list, selection Scripts for the PopHumanScan analysis pipeline are avail- signatures seem to be heterogeneous and a detailed analysis able as Jupyter Notebooks at https://github.com/BGD- of each candidate region is required to understand the real UAB/PopHumanScan. All data, tools and support re- story under each selective event. sources provided by the PopHumanScan database are freely D1088 Nucleic Acids Research, 2019, Vol. 47, Database issue Downloaded from https://academic.oup.com/nar/article-abstract/47/D1/D1080/5134333 by guest on 05 June 2020 Figure 4. Signatures of selection detected at the lactase locus. The distribution of iHS values for the CEU population in 10-kb windows along chromosome 2 are displayed to the left; windows with a P-value < 0.0005 or P-value < 0.005 in the empirical distribution are highlighted. The candidate region under selection including the LCT gene is zoomed-in to the right, where all significant signatures at four different statistics spanning three different signature types are represented: Protein Changes, Site Frequency Spectrum (SFS) and Linkage Disequilibrium (LD). Signatures in each population are colored according to its meta-population: Europe (EUR), Africa (AFR), South-Asia (SAS) and East-Asia (EAS). available at https://pophumanscan.uab.cat. Log-in informa- REFERENCES tion to contribute data to PopHumanScan is available upon 1. Nielsen,R., Akey,J.M., Jakobsson,M., Pritchard,J.K., Tishkoff,S. and request. Willerslev,E. (2017) Tracing the peopling of the world through genomics. Nature, 541, 302–310. 2. Racimo,F., Sankararaman,S., Nielsen,R. and Huerta-Sánchez,E. (2015) Evidence for archaic adaptive introgression in humans. Nat. SUPPLEMENTARY DATA Rev. Genet., 16, 359–371. 3. Fan,S., Hansen,M.E.B., Lo,Y. and Tishkoff,S.A. (2016) Going global Supplementary Data are available at NAR Online. by adapting local: a review of recent human adaptation. Science, 354, 54–59. 4. Sabeti,P.C., Schaffner,S.F., Fry,B., Lohmueller,J., Varilly,P., Shamovsky,O., Palma,A., Mikkelsen,T.S., Altshuler,D. and ACKNOWLEDGEMENTS Lander,E.S. (2006) Positive natural selection in the human lineage. Science, 312, 1614–1620. We thank Carla Giner for helpful comments on the PopHu- 5. Hinds,D.A., Stuve,L.L., Nilsen,G.B., Halperin,E., Eskin,E., manScan data and implementation, and Esteve Sanz for Ballinger,D.G., Frazer,K.A. and Cox,D.R. (2005) Whole-genome patterns of common DNA variation in three human populations. help with the informatics infrastructure in which PopHu- Science, 307, 1072–1079. manScan is implemented. We also thank two anonymous 6. International HapMap Consortium. (2005) A haplotype map of the referees for very helpful comments on the PopHumanScan human genome. Nature, 437, 1299–1320. implementation and manuscript. 7. International HapMap Consortium, Frazer,K, Ballinger,D.G., Cox,D.R., Hinds,D.A., Stuve,L.L., Gibbs,R.A., Belmont,J.W., Boudreau,A., Hardenbol,P. et al. (2007) A second generation human haplotype map of over 3.1 million SNPs. Nature, 449, 851–861. FUNDING 8. International HapMap 3 Consortium, Altshuler,D.M., Gibbs,R.A., Peltonen,L., Altshuler,D.M., Gibbs,R.A., Peltonen,L., Ministerio de Economı́a y Competitividad (Spain) Dermitzakis,E., Schaffner,S.F., Yu,F. et al. (2010) Integrating [CGL2017-89160P to Mauro Santos, A.B.]; AGAUR common and rare genetic variation in diverse human populations. Nature, 467, 52–58. (Generalitat de Catalunya) [2017SGR-1379 to Alfredo 9. Sabeti,P.C., Varilly,P., Fry,B., Lohmueller,J., Hostetter,E., Ruiz]; Secretaria d’Universitats i Recerca del Departament Cotsapas,C., Xie,X., Byrne,E.H., McCarroll,S.A., Gaudet,R. et al. d’Empresa i Coneixement (Generalitat de Catalunya) (2007) Genome-wide detection and characterization of positive [FI-DGR2015 to M.C-.Z.]; Departament de Genètica selection in human populations. Nature, 449, 913–918. 10. Akey,J.M. (2009) Constructing genomic maps of positive selection in i de Microbiologia (UAB) [PIF to J.M-.M.]; Servei de humans: where do we go from here? Genome Res., 19, 711–722. Genòmica i Bioinformàtica de la UAB. Funding for open 11. Kelley,J.L., Madeoy,J., Calhoun,J.C., Swanson,W. and Akey,J.M. access charge: Ministerio de Economı́a y Competitividad (2006) Genomic signatures of positive selection in humans and the (Spain). limits of outlier approaches. Genome Res., 16, 980–989. Conflict of interest statement. None declared. Nucleic Acids Research, 2019, Vol. 47, Database issue D1089 12. The 1000 Genomes Project Consortium, Abecasis,G.R., Auton,A., 35. Derrien,T., Estellé,J., Marco Sola,S., Knowles,D.G., Raineri,E., Brooks,L.D., DePristo,M.A., Durbin,R.M., Handsaker,R.E., Guigó,R. and Ribeca,P. (2012) Fast computation and applications of Kang,H.M., Marth,G.T. and McVean,G.A. (2012) An integrated map genome mappability. PLoS One, 7, e30377. of genetic variation from 1,092 human genomes. Nature, 491, 56–65. 36. Gardiner-Garden,M. and Frommer,M. (1987) CpG islands in 13. The 1000 Genomes Project Consortium. (2015) A global reference for vertebrate genomes. J. Mol. Biol., 196, 261–282. human genetic variation. Nature, 526, 68–74. 37. Pennacchio,L.A., Ahituv,N., Moses,A.M., Prabhakar,S., 14. Johnson,K.E. and Voight,B.F. (2018) Patterns of shared signatures of Nobrega,M.A., Shoukry,M., Minovitsky,S., Dubchak,I., Holt,A., recent positive selection across human populations. Nat. Ecol. Evol., Lewis,K.D. et al. (2006) In vivo enhancer analysis of human 2, 713–720. conserved non-coding sequences. Nature, 444, 499–502. 15. Li,M.J., Wang,L.Y., Xia,Z., Wong,M.P., Sham,P.C. and Wang,J. 38. Lesurf,R., Cotto,K.C., Wang,G., Griffith,M., Kasaian,K., (2014) dbPSHP: a database of recent positive selection across human Jones,S.J.M., Montgomery,S.B., Griffith,O.L. and Open Regulatory populations. Nucleic Acids Res., 42, D910–D916. Annotation Consortium. (2016) ORegAnno 3.0: a community-driven 16. Pybus,M., Dall’Olio,G.M., Luisi,P., Uzkudun,M., resource for curated regulatory annotation. Nucleic Acids Res., 44, Carreño-Torres,A., Pavlidis,P., Laayouni,H., Bertranpetit,J. and D126–D132. Downloaded from https://academic.oup.com/nar/article-abstract/47/D1/D1080/5134333 by guest on 05 June 2020 Engelken,J. (2014) 1000 Genomes Selection Browser 1.0: a genome 39. Pollard,K.S., Hubisz,M.J., Rosenbloom,K.R. and Siepel,A. (2010) browser dedicated to signatures of natural selection in modern Detection of nonneutral substitution rates on mammalian humans. Nucleic Acids Res., 42, D903–D909. phylogenies. Genome Res., 20, 110–121. 17. Pybus,M., Luisi,P., Dall’Olio,G.M., Uzkudun,M., Laayouni,H., 40. MacDonald,J.R., Ziman,R., Yuen,R.K.C., Feuk,L. and Scherer,S.W. Bertranpetit,J. and Engelken,J. (2015) Hierarchical boosting: a (2014) The database of genomic variants: a curated collection of machine-learning framework to detect and classify hard selective structural variation in the human genome. Nucleic Acids Res., 42, sweeps in human populations. Bioinformatics, 31, 3946–3952. D986–D992. 18. Casillas,S., Mulet,R., Villegas-Mirón,P., Hervas,S., Sanz,E., 41. Bao,W., Kojima,K.K. and Kohany,O. (2015) Repbase Update, a Velasco,D., Bertranpetit,J., Laayouni,H. and Barbadilla,A. (2018) database of repetitive elements in eukaryotic genomes. Mob. DNA, 6, PopHuman: the human population genomics browser. Nucleic Acids 11. Res., 46, D1003–D1010. 42. Bailey,J.A., Gu,Z., Clark,R.A., Reinert,K., Samonte,R.V., 19. McDonald,J.H. and Kreitman,M. (1991) Adaptive protein evolution Schwartz,S., Adams,M.D., Myers,E.W., Li,P.W. and Eichler,E.E. at the Adh locus in Drosophila. Nature, 351, 652–654. (2002) Recent segmental duplications in the human genome. Science, 20. Charlesworth,B. (1994) The effect of background selection against 297, 1003–1007. deleterious mutations on weakly selected, linked variants. Genet. Res., 43. Benson,G. (1999) Tandem repeats finder: a program to analyze DNA 63, 213–227. sequences. Nucleic Acids Res., 27, 573–580. 21. Smith,N.G. and Eyre-Walker,A. (2002) Adaptive protein evolution in 44. Ward,L.D. and Kellis,M. (2016) HaploReg v4: systematic mining of Drosophila. Nature, 415, 1022–1024. putative causal variants, cell types, regulators and target genes for 22. Harrow,J., Frankish,A., Gonzalez,J.M., Tapanari,E., Diekhans,M., human complex traits and disease. Nucleic Acids Res., 44, Kokocinski,F., Aken,B.L., Barrell,D., Zadissa,A., Searle,S. et al. D877–D881. (2012) GENCODE: the reference human genome annotation for The 45. Zerbino,D.R., Achuthan,P., Akanni,W., Amode,M.R., Barrell,D., ENCODE Project. Genome Res., 22, 1760–1774. Bhai,J., Billis,K., Cummins,C., Gall,A., Girón,C.G. et al. (2018) 23. Casillas,S. and Barbadilla,A. (2017) Molecular Population Genetics. Ensembl 2018. Nucleic Acids Res., 46, D754–D761. Genetics, 205, 1003–1035. 46. The Gene Ontology Consortium. (2017) Expansion of the gene 24. Voight,B.F., Kudaravalli,S., Wen,X. and Pritchard,J.K. (2006) A map ontology knowledgebase and resources. Nucleic Acids Res., 45, of recent positive selection in the human genome. PLOS Biol., 4, e72. D331–D338. 25. Shriver,M.D., Kennedy,G.C., Parra,E.J., Lawson,H.A., Sonpar,V., 47. Brown,G.R., Hem,V., Katz,K.S., Ovetsky,M., Wallin,C., Huang,J., Akey,J.M. and Jones,K.W. (2004) The genomic distribution Ermolaeva,O., Tolstoy,I., Tatusova,T., Pruitt,K.D., Maglott,D.R. of population substructure in four populations using 8,525 autosomal et al. (2015) Gene: a gene-centered information resource at NCBI. SNPs. Hum. Genomics, 1, 274–286. Nucleic Acids Res., 43, D36–D42. 26. Fay,J.C. and Wu,C.I. (2000) Hitchhiking under positive Darwinian 48. UniProt Consortium,T. (2018) UniProt: the universal protein selection. Genetics, 155, 1405–1413. knowledgebase. Nucleic Acids Res., 46, D158–D169. 27. Fu,Y.X. and Li,W.H. (1993) Statistical tests of neutrality of 49. Papatheodorou,I., Fonseca,N.A., Keays,M., Tang,Y.A., Barrera,E., mutations. Genetics, 133, 693–709. Bazant,W., Burke,M., Füllgrabe,A., Fuentes,A.M.-P., George,N. 28. Tajima,F. (1989) Statistical method for testing the neutral mutation et al. (2018) Expression Atlas: gene and protein expression across hypothesis by DNA polymorphism. Genetics, 123, 585–595. multiple studies and organisms. Nucleic Acids Res., 46, D246–D251. 29. Wright,S. (1950) Genetical structure of populations. Nature, 166, 50. Amberger,J.S., Bocchini,C.A., Schiettecatte,F., Scott,A.F. and 247–249. Hamosh,A. (2015) OMIM.org: Online Mendelian Inheritance in 30. Hudson,R.R., Slatkin,M. and Maddison,W.P. (1992) Estimation of Man (OMIM®), an online catalog of human genes and genetic levels of gene flow from DNA sequence data. Genetics, 132, 583–589. disorders. Nucleic Acids Res., 43, D789–D798. 31. Mackay,T.F.C., Richards,S., Stone,E.A., Barbadilla,A., Ayroles,J.F., 51. Koscielny,G., An,P., Carvalho-Silva,D., Cham,J.A., Fumis,L., Zhu,D., Casillas,S., Han,Y., Magwire,M.M., Cridland,J.M. et al. Gasparyan,R., Hasan,S., Karamanis,N., Maguire,M., Papa,E. et al. (2012) The Drosophila melanogaster genetic reference panel. Nature, (2017) Open Targets: a platform for therapeutic target identification 482, 173–178. and validation. Nucleic Acids Res., 45, D985–D994. 32. Casper,J., Zweig,A.S., Villarreal,C., Tyner,C., Speir,M.L., 52. Lyne,R., Sullivan,J., Butano,D., Contrino,S., Heimbach,J., Hu,F., Rosenbloom,K.R., Raney,B.J., Lee,C.M., Lee,B.T., Karolchik,D. Kalderimis,A., Lyne,M., Smith,R.N., Štěpán,R. et al. (2015) et al. (2018) The UCSC genome browser database: 2018 update. Cross-organism analysis using InterMine. Genes, 53, 547–560. Nucleic Acids Res., 46, D762–D769. 53. Bhérer,C., Campbell,C.L. and Auton,A. (2017) Refined genetic maps 33. Vernot,B., Tucci,S., Kelso,J., Schraiber,J.G., Wolf,A.B., reveal sexual dimorphism in human meiotic recombination at Gittelman,R.M., Dannemann,M., Grote,S., McCoy,R.C., Norton,H. multiple scales. Nat. Commun., 8, 14994. et al. (2016) Excavating neandertal and denisovan DNA from the 54. Mi,H., Huang,X., Muruganujan,A., Tang,H., Mills,C., Kang,D. and genomes of melanesian individuals. Science, 352, 235–239. Thomas,P.D. (2017) PANTHER version 11: expanded annotation 34. Martı́nez-Fundichely,A., Casillas,S., Egea,R., Ràmia,M., data from Gene Ontology and Reactome pathways, and data analysis Barbadilla,A., Pantano,L., Puig,M. and Cáceres,M. (2014) InvFEST, tool enhancements. Nucleic Acids Res., 45, D183–D189. a database integrating information of polymorphic inversions in the human genome. Nucleic Acids Res., 42, D1027–D1032.