ARTICLE https://doi.org/10.1038/s41467-019-12173-x OPEN Evolutionary and functional impact of common polymorphic inversions in the human genome Carla Giner-Delgado 1,2,13, Sergi Villatoro 1,13, Jon Lerga-Jaso1,13, Magdalena Gayà-Vidal 1,3, Meritxell Oliva1, David Castellano 1, Lorena Pantano 1, Bárbara D. Bitarello 4, David Izquierdo 1, Isaac Noguera1, Iñigo Olalde5, Alejandra Delprat1, Antoine Blancher6,7, Carles Lalueza-Fox 5, Tõnu Esko8, Paul F. O’Reilly 9, Aida M. Andrés4,10, Luca Ferretti 11, Marta Puig 1 & Mario Cáceres 1,12 1234567890():,; Inversions are one type of structural variants linked to phenotypic differences and adaptation in multiple organisms. However, there is still very little information about polymorphic inversions in the human genome due to the difficulty of their detection. Here, we develop a new high-throughput genotyping method based on probe hybridization and amplification, and we perform a complete study of 45 common human inversions of 0.1–415 kb. Most inversions promoted by homologous recombination occur recurrently in humans and great apes and they are not tagged by SNPs. Furthermore, there is an enrichment of inversions showing signatures of positive or balancing selection, diverse functional effects, such as gene dis- ruption and gene-expression changes, or association with phenotypic traits. Therefore, our results indicate that the genome is more dynamic than previously thought and that human inversions have important functional and evolutionary consequences, making possible to determine for the first time their contribution to complex traits. 1 Institut de Biotecnologia i de Biomedicina, Universitat Autònoma de Barcelona, Bellaterra, Barcelona 08193, Spain. 2 Departament de Genètica i de Microbiologia, Universitat Autònoma de Barcelona, Bellaterra, Barcelona 08193, Spain. 3 CIBIO/InBIO Research Center in Biodiversity and Genetic Resources, Universidade do Porto, Vairão, Distrito do Porto 4485-661, Portugal. 4 Department of Evolutionary Genetics, Max Planck Institute for Evolutionary Anthropology, Leipzig, Saxony 04103, Germany. 5 Institute of Evolutionary Biology, CSIC-Universitat Pompeu Fabra, Barcelona 08003, Spain. 6 Laboratoire d’immunologie, CHU de Toulouse, IFB Hôpital Purpan, Toulouse 31059, France. 7 Centre de Physiopathologie Toulouse-Purpan (CPTP), Université de Toulouse, Centre National de la Recherche Scientifique (CNRS), Institut National de la Santé et de la Recherche Médicale (Inserm), Université Paul Sabatier (UPS), Toulouse 31024, France. 8 Estonian Genome Center, Institute of Genomics, University of Tartu, Tartu 51010, Estonia. 9 Social, Genetic, and Developmental Psychiatry Centre, Institute of Psychiatry, Psychology and Neuroscience, King’s College London, London SE5 8AF, UK. 10 UCL Genetics Institute, Department of Genetics, Evolution and Environment, University College London, London WC1E 6BT, UK. 11 Big Data Institute, Li Ka Shing Centre for Health Information and Discovery, University of Oxford, Oxford OX3 7LF, UK. 12 ICREA, Barcelona 08010, Spain. 13These authors contributed equally: Carla Giner-Delgado, Sergi Villatoro, Jon Lerga-Jaso. Correspondence and requests for materials should be addressed to M.C. (email: mcaceres@icrea.cat) NATURE COMMUNICATIONS | (2019)10:4222 | https://doi.org/10.1038/s41467-019-12173-x | www.nature.com/naturecommunications 1 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-019-12173-x I n the last decade a great effort has been devoted to char- diverse range of inversions in multiple individuals are necessary acterizing all the variation in the human genome1–5, which to determine their functional and evolutionary impact. opens the door to determining the genetic basis of phenotypic Here, we have developed a new high-throughput genotyping traits and disease susceptibility. Nevertheless, despite the initial method and we have characterized in detail 45 common poly- expectations, a significant fraction of the genetic risk for common morphic inversions. By combining accurate inversion genotypes and complex diseases is still unexplained6,7. Furthermore, not all in 551 individuals of different populations and the available variants have been studied at the same level of detail. In parti- genomic information, we show that a large fraction of inversions cular, inversions are a type of structural variant that changes the are not linked to other variants and have occurred recurrently. In orientation of a genomic segment, usually without gain or loss of addition, several of them have signatures of selection and/or DNA, and they often have highly-identical inverted repeats (IRs) functional effects, emphasizing the role of inversions in the at their breakpoints. These characteristics make inversion detec- human genome. tion very challenging, even with next-generation sequencing methods, and they have been largely overlooked in humans8,9. Genome-wide inversion discovery has been typically based Results on genome sequence comparison10,11 or paired-end mapping High-throughput genotyping of inversions. We focused on a (PEM)4,12, although recent studies have exploited newer techni- representative set of 45 paracentric inversions from the InvFEST ques that could be especially useful for inversion detection, such database19, which comprised most of those experimentally vali- as long-read sequencing13–15, Strand-seq16, BioNano optical dated by PCR-based techniques when the project started19 and maps17, or a combination of them18. In most cases around corresponds approximately to half of the estimated number of 100–200 inversions have been predicted, with a maximum of 786 real variants with >5% frequency in human populations4,14 predictions in the 1000 Genomes Project (1000GP)4,19. However, (Supplementary Fig. 1, Supplementary Data 1). These inversions these methods are not suitable for high-throughput genotyping, were originally detected in the comparison of the hg18 and HuRef and with few exceptions4,16,20, just a reduced number of indivi- genome assemblies10 or a fosmid PEM survey in nine indivi- duals (1–15) have been analyzed. Moreover, the presence of duals12, and between 1 and 36 of them have been identified in repetitive sequences at the breakpoints influences the inversions different recent studies (Supplementary Data 1). The main lim- that can be detected by each technique and results in high error itation for inversion genotyping was due to breakpoint IRs, that rates for inversion validation compared to other variants4,18,20,21. had to be of less than 25–30 kb and with target sites for certain Apart from the intensely-studied 17q21.31 and 8p23.1 restriction enzymes at both sides but not within the IRs26, which inversions22,23, genotyping efforts have been restricted to a small excluded previously-known large inversions mediated by big number of inversions and samples. For example, four other large repeat blocks19. Overall, the studied inversions are located inversions have been genotyped by FISH24 and five smaller throughout the genome (37 in the autosomes, 7 in chr. X and 1 in inversions by PCR25 in 27 and 42 individuals of four populations, chr. Y), with sizes ranging from 83 bp to 415 kb. Also, 24 (53%) respectively. In addition, PCR and inverse PCR (iPCR) have been have been generated by non-allelic homologous recombination used for targeted studies of 34 inversions in 70–90 Europeans21,26 (NAHR) between >90% identical IRs (from 654 bp to 24.2 kb). and a more worldwide characterization of three inversions25,27,28. The rest (47%) were probably generated by non-homologous Also, although inversion genotypes might be predicted based on mechanisms (NH), including 18 with small deletions or insertions SNP data, these methods can only detect inversions above a in the derived allele that may have been created in a single certain size or associated with specific SNP combinations and the complex FoSTeS/MMBIR event4,21 (Supplementary Fig. 1, Sup- error rate can be high21,23,25,29. Therefore, it is not yet clear how plementary Data 1). many polymorphic inversions really exist in humans and very Of those, 41 inversions were genotyped simultaneously using little is known about their global frequency and distribution19. high-throughput assays derived from the multiplex ligation- Actually, inversions have been a model in evolutionary biology dependent probe amplification (MLPA) technique42. For inver- for almost 90 years30,31 and there are numerous examples of their sions with simple breakpoint sequences (17), we carried out phenotypic consequences and adaptive significance in diverse directly custom MLPA assays with minor modifications. How- organisms, from plants to birds32. One of their main effects is ever, for inversions with repetitive sequences at the breakpoints related to recombination, since single crossovers within the (24), which are difficult to detect by most techniques, we inverted region in heterozygotes generate unbalanced gametes developed a new method combining the principles of iPCR26 and and, at the same time, the resulting inhibition of recombination MLPA42 named iMLPA. In both cases, two pairs of oligonucleo- could protect favorable allele combinations30,31. In addition, tide probes were used to interrogate the two alternative inversion breakpoints can directly alter the expression patterns of orientations for each inversion, orientation 1 (O1) and orienta- adjacent genes9,33. tion 2 (O2) (Fig. 1). Four additional inversions not initially From the little information available, it is clear that inversions included in the MLPA-like assays were tested independently by can have important consequences in humans9. Inversions are PCR or iPCR (Supplementary Data 1). The 45 inversions were associated with haemophilia A34, increased risk of neurodegen- genotyped in 551 individuals from seven populations studied in erative diseases35–37, autoimmune diseases23,29 or mental dis- HapMap and 1000GP1,3 with African (AFR) (YRI, LWK), orders38. They could also predispose to other genomic European (EUR) (CEU, TSI), South-Asian (SAS) (GIH) or rearrangements with negative phenotypic consequences in the East-Asian (EAS) (CHB, JPT) ancestry, here referred as offspring9. Moreover, there is evidence that the 17q21.31 inver- population groups (Supplementary Table 1). sion increases the fertility of the carriers and has been positively MLPA and iMLPA inversion genotypes were carefully selected in Europeans22. Finally, inversions have been shown to validated through several analyses and quality controls (Fig. 1): affect gene expression23,28,29,39. However, most of these effects (1) comparison with 3377 available genotypes19 (see Supplemen- are associated with just the two best-known inversions. Attempts tary Data 1 for data source); (2) identification of potential iPCR to associate inversions with gene-expression and phenotypic or iMLPA problems caused by restriction site polymorphisms; variation in large datasets have been limited to those with simple and (3) association between inversions and other variants (see breakpoints, and only a couple of additional candidates have been below). As part of the validation, we repeated by PCR or iPCR identified so far4,40,41. Thus, specific genotyping studies of a 2160 extra genotypes with discrepancies or possible errors plus 2 NATURE COMMUNICATIONS | (2019)10:4222 | https://doi.org/10.1038/s41467-019-12173-x | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-019-12173-x ARTICLE A B C D Orientation 1 A C B D Orientation 2 45 inversions PCR 1 17 24 3 iPCR MLPA iMLPA A B A B A B Restriction enzyme site A C A C A C Genotyping in 551 individuals CEU CHB TSI JPT YRI GIH LWK Quality controls Comparsion with Restriction site Haplotype and tag SNPs published genotypes polymorphism analysis PCR validation Final genotypes (24,355) Fig. 1 Schematic representation of inversion genotyping strategy. High-throughput genotyping of the 45 inversions in 551 individuals from different populations was done by MLPA (17), iMLPA (24), regular PCR (1) and iPCR (3). To avoid confusion about the ancestral status, the two inversion orientations have been named as 1 (AB and CD breakpoints) and 2 (AC and BD breakpoints), using the hg18 genome assembly as reference (Supplementary Data 1). In MLPA and iMLPA, two pairs of oligonucleotide probes (represented in top of the genome sequence) that are able to hybridize contiguously to the target region through a specific sequence complementary to the genome (light grey) were used to interrogate the two alternative orientations of each inversion. These probes, which include a stuffer sequence of variable length (black), are ligated together in a subsequent step and the resulting products are amplified for all the analyzed inversions at the same time with common primers (in green). IRs or other repetitive sequences at the breakpoints are represented as a dark pointed rectangle NATURE COMMUNICATIONS | (2019)10:4222 | https://doi.org/10.1038/s41467-019-12173-x | www.nature.com/naturecommunications 3 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-019-12173-x NH NAHR 1.00 Frequency DAF 0.75 0.50 MAF 0.25 0.00 Present Age 1 Mya 2 Mya FST Selection LSFS Pos. LSFS Bal. NCD1 NCD2 Gene mut. Function Geuv. eQTLs 1 1 18 2 3 1 5 4 GTEx eQTLs 2 2 21 1 1 3 1 4 7 4 1 4 9 16 GWAS hits GWAS enrich. HsInv0097 HsInv0379 HsInv0284 HsInv0068 HsInv0092 HsInv0098 HsInv0260 HsInv0102 HsInv0095 HsInv0059 HsInv0058 HsInv0031 HsInv0105 HsInv0041 HsInv0045 HsInv0409 HsInv1053 HsInv0201 HsInv0006 HsInv1116 HsInv0004 HsInv0061 HsInv1051 HsInv0072 HsInv0030 HsInv0341 HsInv0209 HsInv0340 HsInv0396 HsInv0347 HsInv0040 HsInv0124 HsInv0266 HsInv0055 HsInv0832 HsInv0241 HsInv0397 HsInv0344 HsInv0114 HsInv0393 HsInv0403 HsInv0374 HsInv1124 HsInv0389 HsInv0278 LWK CHB CEU GIH Strong None Population Evidence YRI JPT TSI Weak ND Fig. 2 Evolutionary and functional information for human polymorphic inversions. Frequency: inversion frequencies in the 480 unrelated individuals from seven populations, showing either the derived allele frequency (DAF) if the ancestral orientation is known or the minor allele frequency (MAF) according to the global frequency of the inversion otherwise (which enables the MAF to be higher than 0.5 in specific populations). Age: average age for 22 inversions in which it can be calculated from the divergence between orientations using three different substitution rate estimates. Selection: summary of inversion selection signatures from FST, LSFS (positive or balancing selection) and NCD1 and NCD2 tests. Populations where the signal was detected are indicated by different colors in the corners of each cell, with alternating vertical and diagonal crosses to avoid visual overlap. Criteria for classification of strong and weak selection evidence are explained in Supplementary Data 7. Function: functional effects of inversions summarizing direct gene mutations, which include gene or transcript disruption (strong) or exchange of genic sequences (weak) (Supplementary Table 3), eQTLs in the GEUVADIS or GTEx datasets (showing the number of affected genes and labeled as strong if inversion is lead eQTL for at least one gene) (Supplementary Data 8 and Supplementary Data 9), and association with GWAS hits (strong if for one of the associations P < 1 × 10−6) (Supplementary Table 4) or GWAS signal enrichment (strong if enrichment empirical test P < 0.01 in both GWAS databases). Source data are provided as a Source Data file others randomly selected, including the whole set of 551 African and non-African individual carry the O2 allele, individuals for the three inversions that had the highest error respectively, for 28 and 24 inversions. rate (HsInv0045, HsInv0055 and HsInv0340) (Supplementary Data 1). This showed that the new inversion genotyping technique is very robust, with missing data (0.7%) and genotype Nucleotide variation and haplotype distribution. Thanks to the errors for MLPA (0.1%) and iMLPA (0.9%) accumulating mainly accurate genotypes, we were able to explore the linkage dis- in specific problematic inversions or DNA samples (Supplemen- equilibrium (LD) between inversions and neighboring variants tary Fig. 2A). When compared to the global inversion data from (SNPs and small indels) from HapMap and 1000GP1,3. While the 1000GP4, just 14 of our inversions (31%) were detected, with most NH inversions (20/21) have variants in complete LD (r2 = nine of them having 2.5–71.5% incorrect genotypes and 1) either inside or up to 100 kb from the breakpoint, among the extremely low genotype agreement in the only two inversions 24 NAHR inversions only HsInv0040 and HsInv1051 have at mediated by large IRs in common (Supplementary Fig. 2B). least one such variant (Fig. 3a; Supplementary Table 2). Max- Therefore, we have generated the largest and most accurate imum r2 values between the remaining 22 NAHR inversions and dataset of different types of inversions in humans (Supplementary 1000GP variants range from 0.14 to 0.91 (Fig. 3a). Data 2). Also, we checked the presence of shared SNPs (not including As expected, the 45 inversions show correct genetic transmis- indels) in both orientations. Consistent with recombination sion in the 30 CEU and 30 YRI father-mother-child trios, and inhibition in heterozygotes, most NH inversions do not have allele frequencies do not deviate from Hardy-Weinberg equili- shared SNPs within the inverted region in accessible 1000GP brium in any population (P > 0.01). Minor allele frequencies positions or HapMap data, with the exception of a few (MAF) range globally from 0.5% to 49.7%, with 41 inversions individual SNPs that might be genotype errors or gene spread through several population groups and only the four with conversion events (Fig. 3a; Supplementary Table 2). Outside the lowest frequencies being present in a single population group of the inversion, the average proportion of shared 1000GP SNPs (HsInv0097, HsInv0284 and HsInv1051 in Africa; HsInv0379 in increases progressively after the last fixed variant, up to ~20% East Asia) (Fig. 2; Supplementary Table 1). On average, an (Fig. 3b). This allowed us to define an extended area on each 4 NATURE COMMUNICATIONS | (2019)10:4222 | https://doi.org/10.1038/s41467-019-12173-x | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-019-12173-x ARTICLE a NH inversions NAHR inversions c Out of africa R2a (5) 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.21 0.12 0.39 0.23 0.09 1.43 0.30 0.58 0.34 0.29 0.38 0.24 0.62 0.50 0.13 1.55 1.33 2.45 0.69 0.98 0.96 1.24 1.42 5.86 IR/Inv length R1b (39) P-M45 R1a (16) 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 0.60 1.00 1.00 1.00 1.00 1.00 0.63 0.88 0.75 0.81 0.90 0.76 0.71 0.91 0.74 0.66 0.62 0.54 0.29 0.29 0.55 0.36 0.77 0.27 0.50 0.19 0.64 0.14 Max LD (r 2) Q 0.00 O3 (13) 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.01 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.07 0.08 0.12 0.16 0.07 0.17 0.30 0.05 0.26 0.14 0.15 0.19 0.00 0.39 0.18 0.41 0.24 0.33 0.40 0.50 Frac. shared O2 (13) (1000GP) O1 (1) 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.43 0.17 0.00 0.00 0.00 0.00 0.00 0.33 0.55 0.75 1.00 0.00 1.00 0.50 0.83 1.00 1.00 1.00 1.00 1.00 1.00 1.00 Frac. shared N (2) IJK-M523 (HapMap) K Recurrence T (3) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 2 4 1 3 2 3 5 0 1 7 1 0 4 4 3 events L (2) J2 (16) 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.07 0.29 0.05 0.05 0.00 0.41 0.03 0.33 0.25 0.00 0.43 0.11 0.30 0.00 0.46 Freq. chimp J1 (1) I2 (1) 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.23 0.00 0.00 0.31 0.22 0.00 0.08 0.11 0.00 0.15 0.33 0.15 Freq. gorilla I1 (6) H (9) CF-P143 G (2) HsInv0004 HsInv0031 HsInv0092 HsInv0097 HsInv0260 HsInv0284 HsInv0379 HsInv0095 HsInv1053 HsInv0041 HsInv0058 HsInv0045 HsInv0098 HsInv0068 HsInv0105 HsInv1116 HsInv0006 HsInv0102 HsInv0059 HsInv0201 HsInv0409 HsInv0040 HsInv1051 HsInv0061 HsInv0114 HsInv0030 HsInv0209 HsInv0389 HsInv0072 HsInv0341 HsInv0347 HsInv0266 HsInv0340 HsInv0124 HsInv0393 HsInv0396 HsInv0374 HsInv0241 HsInv0397 HsInv0832 HsInv0344 HsInv0055 HsInv0278 HsInv0403 HsInv1124 F C3 (2) C5 (9) b Fixed Shared C1 (2) BT-M42 D2 (10/1) of variants 0.4 Fraction E2 (4) DE-M145 0.2 E1b1a1 (61) E1b1b1 (8) 0.0 E1a (5) –20 0 20 –20 0 20 B (1) A 80 70 60 50 40 30 20 10 0 Relative position (kb) Time (kya) Fig. 3 Evidence of unique or recurrent origin of NH and NAHR generated inversions. a Summary of the ratio between breakpoint IR and inversion length, maximum LD with neighboring variants in 1000GP and HapMap, fraction of shared SNPs inside the inverted region from the total number of SNPs analyzed in 1000GP or HapMap (which includes less SNPs and results in higher shared fractions), estimated number of recurrence events from haplotype analysis, and DAF or MAF for the inversions in non-human primates. Gray cells indicate values which could not be calculated. NAHR inversions show clear differences in LD with nearby variants, proportion of shared polymorphisms between orientations, recurrence events, frequency in ape species and other measures, with values that could be associated with higher levels of recurrence represented in stronger shades of green. b Distribution of fixed (brown) and shared (green) variants between the two orientations estimated from inversion genotypes and 1000GP data. The shared/fixed fraction with respect to the total number of polymorphic variants was computed for the whole inverted region plus 10-kb overlapping windows with a 5-kb step size in the flanking regions, and the horizontal axis represents the distance of the window central position to the inversion breakpoints (indicated by dashed lines). Thin lines represent individual inversions and thick ones represent the average for each inversion class. c Overview of the human chr. Y phylogenetic tree showing five different possible inversion events in the HsInv0832 region (red arrows). Divergence dates (bottom scale) and tree topology are based on Poznik et al77. For each chr. Y haplogroup, the number of males genotyped for HsInv0832 with each orientation is indicated in parenthesis (O1 in dark grey and O2 in red). Only the main branches and those including genotyped individuals are shown, with some characteristic mutations indicated in the tree. For haplogroup E, one of the two possible scenarios is represented, with the alternative being two inversion events in E1a and E1b1a1 haplogroups. Source data are provided as a Source Data file side of the inversion of no or little recombination between Based on the results of the different analyses, we inferred the orientations, which ranges from 0 to more than 20 kb minimum number of recurrent inversion events and their (Supplementary Table 2). In contrast, 20 of the 24 NAHR distribution in human populations (Fig. 3a; Supplementary inversions have a considerable number of shared SNPs Data 3). However, this relies on having differentiated haplotype scattered throughout the inverted and flanking regions (Fig. 3a, clusters in which the existence of the alternative orientation b; Supplementary Table 2). cannot be easily explained by other factors (such as gene Next, we visualized the haplotype diversity and distribution conversion or genotype/phasing errors of a few variants). Another across orientations and populations using haplotype networks problem is that recombination generates mixed sequences and a new representation integrating a hierarchical clustering of between haplotype groups and makes it difficult to accurately haplotypes and the differences between them (Supplementary quantify recurrence. Thus, a nice example is chr. Y inversion Fig. 3). This analysis was focused on 1000GP data, although HsInv0832, in which there is no recombination. We used consistent results were obtained with HapMap SNPs. After taking available haplogroup information of 232 males from the known into account possible phasing errors, two clearly differentiated chr. Y genealogical tree (Supplementary Data 4) to identify five patterns were observed again. In the 20 NH inversions with independent inversion events in the last ~60,000 years (Fig. 3c). sequence variation information, the haplotypes of one of the This results in an inversion rate of 5.35 × 10–5 per generation (see orientations tend to cluster together, supporting a unique origin Methods), which is ~1000 times higher than that of single bases. of the inversion (Supplementary Fig. 3). In NAHR inversions, this For the inversions in which it is possible to quantify recurrence, is true only for HsInv0040, HsInv0061 and HsInv1051, with the we estimated a total of 40 additional inversion and re-inversion other 21 having O1 and O2 haplotypes mixed throughout the events (ranging for each inversion from 0 to 7 with an average of network and hierarchical cluster, including in many cases 2.2) (Supplementary Data 3). Of those, 12 are distributed globally identical haplotypes with both orientations (Supplementary and could precede the out-of-Africa migration, 17 are restricted Fig. 3). Such pattern is consistent with a multiple origin of these to African individuals, and 11 probably appeared more recently in inversions and explains the absence of fixed SNPs and the high non-African populations. The fact that many of the recurrence number of shared SNPs between orientations. events are shared by several individuals indicates that they are not NATURE COMMUNICATIONS | (2019)10:4222 | https://doi.org/10.1038/s41467-019-12173-x | www.nature.com/naturecommunications 5 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-019-12173-x a All inversions (MAF) b c AFR *** NH inversions NAHR inversions EAS  = 0.18 * 1.00 1.00 EUR  = 0.32 *** SAS *  = – 0.17 ***  = 0.15 Derived allele frequency (DAF) All ** Minor allele frequency (MAF)  = – 0.13 **  = 0.14 * NH inversions (DAF) 0.75 0.75  = – 0.12 *** Population group AFR  = – 0.12 *** EAS EUR 0.50 c = – 0.13 *** 0.50 SAS c = – 0.03 All c = – 0.08 NAHR inversions (MAF) c = – 0.07 ** AFR *** 0.25 0.25 EAS EUR SAS * All ** 0.00 0.00 0.1 0.2 0.3 0.4 0.5 –2 –1 0 1 2 (0,1/5] (1/5,2/5] (2/5,1] (1,5] Mean frequency log10 (Genetic length (4Ner)) IR size/inversion size Fig. 4 Determinants of inversion frequency in human populations. a Observed mean inversion frequency per population group and mechanism of generation (black dots) compared with that expected from the detection method simulations using a null distribution of 10,000 SNP samples. Graphs represent DAF only if the ancestral orientation is known for all the inversions included and MAF otherwise. b Logarithmic robust regression between DAF of NH inversions and genetic length (measured in 4Ner units) in the different population groups, showing a significant negative trend for all inversions (solid line; β regression coefficient) and for inversions larger than 2 kb to correct in part the detection bias against low-frequency small inversions (dashed line; βc regression coefficient). c Boxplots showing a positive relationship between the frequency of the minor allele in all populations together and IR/ inversion length ratio for NAHR inversions (centre, median; bounds of box, first and third quartiles; and whiskers, smallest/largest value within 1.5 interquartile range). β indicates the robust regression coefficient in every population group. In b and c dots indicate the frequency of each inversion with population groups represented with the same colors as in a. Two-tailed t-test (a) and robust regression t-test (b and c) P values: *P < 0.05; **P < 0.01; and ***P < 0.001. Source data are provided as a Source Data file artifacts from lymphoblastoid cell line (LCL) culture. In addition, Moreover, we checked the presence of the breakpoint as part of the validation of inversion genotypes, all the recurrence sequences of 19 NH inversions without IRs in available events were confirmed by checking at least one of the supporting Neanderthal, Denisovan and two ancient modern human individuals by PCR/iPCR. We have therefore extended consider- genomes (Supplementary Data 6). Five inversions (HsInv0004, ably the previous recurrence analysis of some of the inversions in HsInv0006, HsInv0201, HsInv0409, and HsInv1116) showed the just the CEU population21,26. derived orientation in the Neanderthal or Denisovan genomes (including two with the derived orientation in both). These inversions are distributed through African and non-African Ancestral orientation and inversion age. The published data on populations, which suggests that they are not the result of the ancestral orientation of 32 of the inversion regions21,26,28 was introgression and they appeared before the divergence of the most complemented and expanded by experimentally testing 42 recent Homo groups, around 550,000–750,000 years ago (ya)43. inversions for which the human or modified assays generated Finally, we dated more precisely 22 unique inversions from the reliable results in a panel of 23 chimpanzees (40 inversions) and sequence divergence between orientations (Fig. 2; Supplementary seven gorillas (41 inversions) (Supplementary Data 5). Inversion Data 6). Six inversions were estimated to have appeared more orientation was also assessed in available genome assemblies of than 1 Mya, including four with the derived orientation in both species plus orangutan and rhesus macaque (Supplementary Neanderthal or Denisovan, plus HsInv0031 and HsInv0058. For Data 6). In total, we could infer the ancestral allele for 29 HsInv0006, the estimated age (407,795–495,470 ya) is slightly inversions, with 15 showing the ancestral and 14 the derived allele more recent than the Neanderthal-Denisova and modern human in the human reference genome. For the 21 NH inversions, divergence, but it is very close to its lower bound. Ages of the rest orientation was consistent in all the non-human samples and of inversions are consistent with their geographical distribution, genomes analyzed, and with existing deletions and insertions with inversions restricted to either African (HsInv0097 and occurring in the derived allele (Supplementary Data 1 and Sup- HsInv0284) or non-African populations (HsInv0379) having a plementary Data 6). In contrast, 14 of the 22 NAHR inversions relatively recent origin. Only in two cosmopolitan inversions age experimentally genotyped were polymorphic in at least one of estimates are lower than population split times (HsInv1053 with a the apes (six in chimpanzees, two in gorillas and six in both negative age and HsInv0095 with 22,582–41,258 ya), probably species) and several had opposite orientation in different primate due to an underestimation of the divergence between orientations assemblies (Supplementary Data 6). This agrees with previous caused by the limited sequence information available. analyses in a much smaller set of samples21,26, but we found five additional polymorphic inversions in each species. In fact, the lower number of polymorphic inversions in gorilla indicates that Analysis of inversion frequencies. To evaluate whether there are more inversion regions might be identified as polymorphic in selective pressures acting globally on inversions, we compared non-human primates by analyzing more individuals. Given the the inversion frequency spectrum with that of SNPs sampled species divergence times, the most likely scenario is that shared according to neutral expectations (see Methods). Although the inversions have appeared independently in chimpanzee and ascertainment bias associated with inversion detection predicts gorilla lineages, providing additional support for inversion an enrichment of high-frequency inversions (~0.20 expected recurrence (Fig. 3a). MAF), the observed frequencies tend to be higher than expected 6 NATURE COMMUNICATIONS | (2019)10:4222 | https://doi.org/10.1038/s41467-019-12173-x | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-019-12173-x ARTICLE a CCDC144B-002 b r 2 with 3 30 15 Expression 2 HsInv1051 HsInv0030 1 INDEL 0 SNP 0.75 –log10(P value) –1 10 –2 SV 0.4 O2/O2 O1/O2 O1/O1 20 –log10(P value) 0 5 0 –1 Mb –500 kb TSS 500 kb 1 Mb 10 Genome position (chr17) 12 3 RP11-326C3.7 2 Expression 1 HsInv0124 0 INDEL CTRB2 CTRB1 –log10(P value) –1 SNP 8 –2 SV 75.21 Mb 75.24 Mb 75.27 Mb 75.30 Mb –3 O2/O2 O1/O2 O1/O1 Genome position (chr16) 4 e 0 Alzheimer′s disease –50 kb –25 kb TSS 25 kb 50 kb 1.0 –log10(P value) Genome position (chr11) 5.5 5.0 r 2 with HsInv0031 4.5 c d 4.0 3.5 O2/O2 3.0 HsInv0058 ** ** 0.5 (n = 1) Expression O1/O2 (n = 16) HsInv0347 ** O1/O1 (n = 70) HsInv0030 * F. trans. 0.0 18.73 Mb 18.75 Mb 0 10 45 85.18 Mb 85.19 Mb 85.20 Mb Genome position (chr17) Deviation in GWAS signals Genome position (chr16) Fig. 5 Examples of inversion functional effects. a Manhattan plots of logarithm-transformed linear regression t-test P values for cis-eQTL associations in LCLs of transcript CCDC144B-002 and gene RP11-326C3.7, showing inversions HsInv1051 (top) and HsInv0124 (bottom) as lead variants, together with boxplots of rank normalized expression and inversion genotype (centre, median; bounds of box, first and third quartiles; whiskers, smallest/largest value within 1.5 interquartile range; and notches, 95% confidence interval of the median). b Manhattan plot of logarithm-transformed GTEx P values of pancreas eQTLs for CTRB1 (red) and CTRB2 (blue) mapping around the HsInv0030 region (orange bar). Top lead eQTLs for each gene were the variants in highest LD with the inversion (rs9928842, r2 = 0.75; rs8057145, r2 = 0.73; black arrows). c Schematic diagram of average RNA-Seq expression profile in HsInv1051 genotypes of GEUVADIS LCL reads mapped to the inverted allele, showing the generation of a fusion transcript (F. trans.) with new 3′ sequences that loses the last two thirds of the gene exons. Inversion breakpoint is indicated in pink. d GWAS Catalog (orange) and GWASdb (blue) signals enrichment within individual inversions and flanking regions (±20 kb). Dot represents the mean and error bars the 0–0.95 confidence interval of the difference between the observed number of GWAS signals and a null distribution from 1000 random genomic regions for each inversion. One-tailed empirical test P value: **P < 0.01; *P < 0.05. e LD (r2) between HsInv0031 (orange bar) and 1000GP variants (black points) or Alzheimer’s disease GWAS signals (diamonds)78. Source data are provided as a Source Data file in all population groups (Fig. 4a). When inversions are sepa- (FST), which can identify positive selection leading to a rapid rated by the generation mechanism, the increase in frequency is increase in allele frequency in some areas. The global FST value was significant only in NAHR inversions (Fig. 4a). We also investi- 0.11 for autosomal inversions, 0.21 for chr. X inversions, and 0.73 gated how diverse genomic variables affect autosomal and chr. X for the one in chr. Y, with the largest frequency variation between inversion frequencies. The most significant predictor is inver- continents. Three inversions were within the top 1% of the FST sion genetic length, which is negatively correlated with DAF of distributions derived from SNPs with the same frequency (Fig. 2; NH inversions and explains 23–55% of the frequency variance in Supplementary Data 7): HsInv0040, with frequency differences the different population groups (Fig. 4b), followed by physical between European populations, and HsInv0389 and HsInv0059, length (19–54% variance explained). For NAHR inversions, only with high frequency in Africa or East Asia, respectively. Eleven 5–14% of MAF variance is accounted by their genetic length, more inversions fell within the top 5% of the empirical FST dis- whereas the ratio between IR and inversion length is positively tribution, and were considered to have weak evidence of selection correlated with inversion frequency in all population groups and (Fig. 2; Supplementary Data 7). explains 13–45% of MAF variance (Fig. 4c). This suggests that Second, we applied a novel test based on the frequency the higher frequency of NAHR inversions might be related to spectrum of linked sites (LSFS), which is well suited to detect their repeated generation by recurrence. deviations from neutrality in low-recombination regions, such as inversions. We used optimized tests to identify positive and long-term balancing selection maintaining a polymorphism Selection on human inversions. To investigate signals of natural in several populations, and significance was assessed empiri- selection acting on specific inversions, we first measured inversion cally using the null LSFS distribution from autosomes. Only 18 frequency differences between populations using the fixation index unique autosomal inversions with nearby perfect tag variants NATURE COMMUNICATIONS | (2019)10:4222 | https://doi.org/10.1038/s41467-019-12173-x | www.nature.com/naturecommunications 7 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-019-12173-x that could be reliably phased were analyzed, including most NH Supplementary Data 9). By searching for eQTL signals in inversions and HsInv0040 (Fig. 2; Supplementary Data 7). The moderate LD (r2 ≥ 0.6) with some of the recurrent inversions, strongest signals (empirical test P < 0.01) were in HsInv0201 for we found additional potential expression differences associated balancing selection and HsInv0006 and HsInv0031 for positive with HsInv0124 and in the genes affected by HsInv0030, which selection. In addition, four other inversions showed weaker exchanges the first exon and promoter of chymotrypsinogen evidence of balancing or positive selection. Consistent with the precursor genes CTRB1 and CTRB2 expressed only in FST results, in HsInv0006 and HsInv0059 positive selection was pancreas21,25, and HsInv0340, which disrupts the long non- detected in those populations with increased DAF (Fig. 2; coding gene LINC00395 expressed in testis (Supplementary Fig. 5; Supplementary Data 7). Supplementary Data 9). In total, 17/27 of inversion-gene As independent confirmation of balancing selection signatures, associations in LCLs were also identified in the smaller GTEx we also used the recently developed non-central deviation sample (Supplementary Fig. 5). statistics, NCD1 and NCD2, which detect site frequency spectrum To assess if inversions were the main cause of the observed shifts towards an equilibrium frequency and an excess of expression changes, we performed a joint eQTL analysis in LCLs polymorphic sites44. However, the results of these tests summar- including our inversions together with SNPs, indels and structural ize the data of all the SNPs in a region and are not necessarily variants from the 1000GP3,4. Two inversions, HsInv0124 and linked to the inversion, as before. Focusing on signals detected in HsInv1051, were the most likely causal variant for two genes and at least three populations, we found respectively four and six three transcripts in LCLs (Fig. 5a). Six other inversions show the inversions with strong and weak signatures of balancing selection highest LD (r2 ≥ 0.9) with variants reported as first or second lead for NCD1 or NCD2 (Fig. 2; Supplementary Data 7). Many of eQTL in a given tissue by the GTEx project (Supplementary these candidates could not be analyzed with the LSFS method Fig. 6). Similarly, for recurrent inversions HsInv0124 and because of the lack of tag SNPs or correspond to low-frequency HsInv0030, eQTL significance in GTEx data increases with LD inversions, such as HsInv1051 and HsInv0209, that are unlikely with the inversions, supporting their causal role (Fig. 5b). In the targets of selection, but consistent results were found for general, some of the strongest effects are related to inversions HsInv0201. affecting exonic sequences, although the consequences can be complex and need to be investigated in detail. For example, HsInv1051 breaks the CCDC144B gene and the apparent Effect of inversions on genes and gene expression. As previously upregulation of specific isoforms (Fig. 5a) is actually due to the described, some of the analyzed inversions can have important creation of a fusion transcript with new sequences at 3′ (Fig. 5c; effects on genes4,9,21,25–28. Although half of our inversions (21/ Supplementary Fig. 7). HsInv0124 is the lead variant for the 45) are located in intergenic regions, three of them invert genes, antisense RNAs RP11-326C3.7 and RP11-326C3.11, which over- eight are located within introns, seven might exchange gene lap respectively the IFITM2 and IFITM3 genes located at the sequences overlapping the IRs at the breakpoints, and six affect breakpoints, and it has opposite effects in the two pairs of genes more directly through the inversion or deletion of an overlapping transcripts (Supplementary Data 8; Supplementary internal exon (HsInv0102, HsInv0201) and the disruption of the Data 9). Also, HsInv0102 removes the RHOH isoforms with the whole gene (HsInv0340, HsInv0379, HsInv1051) or an alternative alternative non-coding exon that gets inverted, but its effect is transcript isoform (HsInv0124) (Fig. 2 and Supplementary masked by a more frequent lead eQTL (SNP rs7699141) acting in Table 3). the same direction. On the other hand, HsInv0058 is associated We measured the effect of the 42 autosomal and chr. X with chr. 6 MHC haplotypes APD, COX, DBB, QBL and SSTO, inversions with MAF > 0.01 on expression of nearby genes (<1 which extend ~4 Mb and harbor important functional differ- Mb away) by linear regression between inversion genotypes and ences47, suggesting that other variants in these haplotypes are LCL transcriptome data from the GEUVADIS consortium45. To responsible for the observed effects. increase statistical power and reliability, the analysis was replicated in two datasets: (1) 173 CEU, TSI and YRI individuals with inversion genotypes; and (2) the complete GEUVADIS set of Inversions and phenotypic traits. The role of inversions in 445 European (358) and African (87) individuals in which the phenotypic variation was investigated using available genome- genotypes of 33 inversions could be imputed accurately wide association studies (GWAS) data. We found a 1.26- and (Supplementary Fig. 4A). Considering the largest sample size 1.95-fold increase in GWAS Catalog48 and GWASdb49 variants in for each inversion, we uncovered eight inversions significantly the inversion and flanking regions (±20 kb). The top inversion associated with LCL expression of 27 genes and 44 transcripts driving this result was the MHC-inversion HsInv0058, but (Supplementary Fig. 4B; Supplementary Data 8), with highly HsInv0030 and HsInv0347 showed similar enrichment of GWAS concordant results for those analyzed in both datasets (7/7 genes hits in both datasets (Fig. 5d; Supplementary Fig. 8A-B). GWAS and 11/12 transcript effects were replicated) (Supplementary signals close to the latter two inversions are consistent with their Fig. 4C). As negative control, no associations were observed by effect on genes, involving type 1 and 2 diabetes, pancreatic cancer, permuting inversion genotypes relative to expression levels insulin secretion, and cholesterol and triglyceride levels for (Supplementary Fig. 4D). Moreover, significant expression effects HsInv0030, and glaucoma and optic disc and nerve characteristics were robust when applying different analysis approaches (see for HsInv0347, which is associated to the expression of c14orf39 Supplementary Fig. 4E-F and Supplementary Methods), and (SIX6OS1) and SIX6, related to eye development. inversions acting as expression quantitative trait loci (eQTLs) We also explored whether inversions were in strong LD (r2 ≥ located significantly closer to the transcription start site (TSS) of 0.8) with known GWAS hits in the population where the the differentially expressed genes (<100 kb) (Supplementary association was reported. That is the case of HsInv0004, which is Fig. 4D). in complete LD in Europeans with a nearly genome-wide Next, we examined inversion expression effects in other tissues significant GWAS SNP related to asthma susceptibility in through variants already reported as eQTLs in the GTEx children and another one associated with body mass index in project46. We found 62 genes with eQTLs in different tissues in asthmatic children (Supplementary Table 4). Moreover, high LD (r2 ≥ 0.8) with 11 of the 26 analyzed inversions, HsInv0006 is linked to schizophrenia in Ashkenazi Jews and including seven not detected in LCL data (Supplementary Fig. 5; glaucoma in Europeans (Supplementary Table 4). Remarkably, 8 NATURE COMMUNICATIONS | (2019)10:4222 | https://doi.org/10.1038/s41467-019-12173-x | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-019-12173-x ARTICLE a Hslnv0201 b 1.00 Hslnv0201 LSFS bal. P < 0.01 0.75 NCD1 P < 0.01 15 r 2 with inversion 0.50 0.25 0.50 0.75 Derived allele frequency 0.25 1.00 10 –log10(P value) Immune response 0.00 Non-infected Salmonella Hslnv0006 1.00 FST AFR P < 0.05 5 0.75 LSFS pos. P < 0.01 0.50 0 0.25 RefSeq SPINK14 SPINK6 MiTranscriptome T291020 0.00 147,520,000 147,540,000 147,560,000 147,580,000 147,600,000 R S R S R EA SA AF EU AM Genome position (chr5) c chr1 205.025 Mb 205.050 Mb 205.075 Mb 205.100 Mb 205.125 Mb 205.150 Mb 205.175 Mb 205.200 Mb 205.225 Mb 205.250 Mb Coding HsInv0006 genes CNTN2 TMEM81 (first intron) TMCC2 RBBP5 DSTYK Non-coding RP11-536L3.4 RP11-383G10.5 genes rs12142514 Decreased expression Glaucoma-GWAS hit 2 Increased expression (inversion r = 1) Fig. 6 Integrative evolutionary and functional analysis of inversion candidates. a Frequency of HsInv0201 and HsInv0006 across worldwide human populations from 1000GP Ph3 (colored dots) and in each population group (colored bars) estimated from global inversion tag SNPs (rs200056603 and rs79619752, respectively), showing a summary of test results for balancing (bal.) and positive (pos.) selection. b Manhattan plot of logarithm-transformed linear regression t-test P values for cis-eQTL associations with expression variation of gene SPINK6 in infected and non-infected primary macrophages52, showing SNPs in complete LD with HsInv0201 (orange bar) as lead eQTLs in Salmonella infection. Variants are represented as circles with varying size depending on the LD (r2) with the inversion. The location of the genes in the region is shown below, including a new SPINK6 isoform discovered by MiTranscriptome assembly79, whose first exon is removed by a deletion associated to the inversion breakpoints (light blue bar). c Diagram of HsInv0006 (orange bar) genomic region showing the effect of the inverted allele on the expression of neighboring genes in different tissues according to the GTEx data and the inversion tag SNP in Europeans associated to increased risk of Glaucoma (Supplementary Table 4). Arrowheads indicate the direction of transcription of the genes. Source data are provided as a Source Data file several of these inversions affect gene expression as well (Fig. 2). Despite the effort to include as many inversions as possible, A good example is HsInv0031, which is associated with lower only those validated and with well-characterized breakpoints levels of FAM92B in cerebellum and is in almost perfect LD in could be analyzed19. In addition, important limitations exist for Europeans (r2 = 0.98) with SNP rs2937145 associated with the study of those mediated by IRs. Although it is not clear how Alzheimer’s disease risk (Fig. 5e; Supplementary Table 4). many of them are real, in the InvFEST database there are ~100 Nevertheless, for many inversions the low LD with SNPs means inversion predictions with 1–25 kb inverted segmental duplica- that any effects would be missed in typical array-based GWAS tions (SDs) at the breakpoints, which could be potentially inter- (Supplementary Fig. 8C). rogated using our method, and ~250 with larger SDs, which are currently out of reach of the methodology19. Therefore, further improvements and complementary strategies, such as the possi- Discussion bility of making directed cuts in specific positions or estimating This work represents the most exhaustive and accurate study of the distance between regions separated by larger repeats, are human inversions so far, including a significant fraction of necessary to expand the range of analyzed inversions. Never- common inversions19, and it is focused on inversions mediated by theless, the present work already offers a good picture of the IRs, which generally escape detection. New genome-wide tech- contribution of common inversions to genetic diversity, adapta- niques are providing a more complete inversion catalogue13–18, tion, and phenotypic traits in humans. but they tend to be laborious and expensive, and population data In particular, we found clearly contrasting patterns for inver- in well-characterized individuals are still scarce. In addition, in sions generated by homologous and non-homologous mechan- many cases the high genotype error rate has precluded identifying isms, supporting a high-degree of recurrence of all inversions tag variants that could be used to infer inversion effects4,40. Thus, mediated by highly-identical IRs except three (two of which have the new genotyping method developed here and the reliable very low frequency). Recurrent inversions are characterized by information from multiple individuals of different populations low LD with other variants, a large amount of shared SNPs and generated are crucial to fill the void in the knowledge of this type shared haplotypes between O1 and O2 chromosomes (Fig. 3). of variation. Also, recurrence is not limited to humans but extends to other NATURE COMMUNICATIONS | (2019)10:4222 | https://doi.org/10.1038/s41467-019-12173-x | www.nature.com/naturecommunications 9 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-019-12173-x great ape species. Similar results had been found in previous frequency. Therefore, this may explain the observed enrichment small-scale studies21,24,26, but thanks to the analysis of many of inversions with functional and selection signals. Future in- more individuals and populations, multiple controls (such as depth studies of the identified candidates and other inversions genotype confirmation by PCR of both breakpoints and different will help uncover their real role in human evolution and the sources of nucleotide variation data) and a detailed analysis of unexplained part of the genetic basis of complex traits. haplotype relationships, we have obtained a better estimate of the independent inversion events. Specifically, we found nine more Methods recurrence events in the five inversions originally predicted as Human and ape DNA samples. We used 550 human samples included in the last recurrent in humans and that seven inversions considered to be phase of the HapMap Project and many of them in 1000GP phase 1 (Ph1) and 3 unique or lacking information in CEU are recurrent as well21,26. (Ph3)1–3, which belong to seven populations of four main population groups (AFR, EUR, SAS, and EAS), plus individual NA15510 of unknown origin (see Supple- This suggests that, like other repeats, IRs are rearrangement mentary Table 1 and Supplementary Data 2 for details). Most individuals were hotspots and the genome is more dynamic than previously unrelated, but 70 are either children of mother-father-child trios (30 in YRI and 30 thought. However, recurrence has been detected only through the in CEU) or cryptic first and second-degree relatives (9 in LWK and 1 in GIH)2,3. sequences associated with the inversion, and we need more direct Genomic DNAs of 70 CEU and 10 YRI samples and NA15510 were extracted from ways to quantify inversion generation rates precisely. LCLs commercialized by the Coriell Cell Repository (Camden, NJ, USA), while the rest of DNAs were acquired from Coriell21,26,28. Chimpanzee and gorilla DNAs As for other types of recurrent changes50,51, the lack of asso- include six already used to genotype most of the inversions from frontal cortex ciation between SNPs and many inversions and poor coverage by tissue samples of the Banc de Teixits Animals de Catalunya (N457/03, Z01/03 and common arrays (Supplementary Fig. 8C) means that their phe- Z02/03) or LCLs from Barcelona Zoo individuals (PTR1211, PTR1213, and notypic effects have been largely missed in GWAS. For example, PTR1215)21,26, plus 19 chimpanzee and five gorilla DNAs extracted from a previously-existing collection of primate LCLs of one of the authors (Supple- half of the NAHR inversions cannot be accurately predicted with mentary Data 5). Ape samples comprise four mother-father-child trios and one typical imputation algorithms (Supplementary Fig. 4A). We have father-son pair in chimpanzees and one father-son pair in gorillas. As for humans, found 23 and 22 inversions with different selection and functional cells were grown in 75-ml flasks to nearly confluency and high-molecular weight signatures, respectively (Fig. 2). More importantly, although not genomic DNA was obtained using a standard phenol-chloroform extraction26. All procedures that involved the use of human and non-human primate samples were all analyses could be applied to every inversion, there is a sig- approved by the Animal and Human Experimentation Ethics Committee (CEEAH) nificant enrichment of inversions with both effects on genes or of the Universitat Autònoma de Barcelona. gene expression and selection signals directly linked to them (Fisher’s exact test P = 0.0320) (see Supplementary Methods). Experimental genotyping of inversions. Initial genotypes for 41 inversions were This combination of the two independent types of evidence obtained by newly developed assays based on the MLPA technique (Supplementary strongly indicates that inversions can have important con- Data 1), which has been widely used for genome copy-number analysis and consists sequences in humans. on the multiplex amplification of fragments of different sizes with common primers that are fluorescently labeled, and their detection by capillary electrophoresis42. One particularly interesting example is HsInv0201, an old Most NH inversions (17) were genotyped simultaneously from 100–150 ng of DNA inversion (>1.5 Mya) with intermediate frequency around the with a slightly modified MLPA protocol using two pairs of probes that bind spe- globe and clear signals of balancing selection (Fig. 6a), which cifically at the breakpoint sequences of each orientation (AB and CD or AC and deletes an exon of SPINK14 and is lead eQTL for two nearby BD), with one of the probes that could be common to both pairs (Fig. 1). For 24 genes (Supplementary Data 9). Moreover, the inversion haplotype inversions with IRs or other repetitive sequences at the breakpoints, we developed a new iMLPA method that uses a combination of iPCR and MLPA. iMLPA requires is the main responsible for the lower SPINK6 expression during some extra processing of the DNA with a restriction enzyme that cuts at each side of immune response to Salmonella infection52,53 (Fig. 6b). In fact, the breakpoint repetitive sequences, followed by self-circularization of all the HsInv0201 eliminates the promoter and first exon of a putative digested DNA molecules together by ligation in diluted conditions with T4 DNA novel SPINK6 transcript (Fig. 6b) and it is in high LD (r2 = 0.971 ligase and DNA purification. Then, MLPA was carried out as usual with a pair of probes that bind specifically at the self-ligation site of the circular molecules from in EUR) with SNPs accounting for plasma levels of SPINK6 each orientation (see Supplementary Methods for MLPA and iMLPA details). protein54. Together with the role of several of the affected genes Supplementary Data 10 lists the sequences and concentrations of the probes used in lung and extracellular mucosas, this suggests that the inversion for MLPA (68) and iMLPA (87). could be related to immune response. In HsInv0006, its particular PCRs and iPCRs were carried out to genotype seven inversions in the 551 human samples (including four not analyzed in the MLPA/iMLPA assays) and to distribution pattern and the selection test results point to positive validate many of the MLPA/iMLPA genotypes (Fig. 1; Supplementary Data 1). selection of the derived allele in Africa (Fig. 6a). Furthermore, the Multiplex PCRs and iPCRs of each inversion were done with primers flanking inversion is located within DSTYK first intron21 and is associated either the breakpoint (PCR) or the self-ligation site of circularized molecules with expression changes in the proximal genes, including DSTYK (iPCR) from the two orientations21,26 (Supplementary Data 11). For the restriction site polymorphism analysis, we downloaded dbSNP (version 142) SNPs and indels upregulation in different tissues (Fig. 6c; Supplementary Fig. 5). around the inversion region (±50 kb) and determined all possible restriction site DSTYK deletion causes pigmentation problems and elevated cell gains or losses affecting the iPCR/iMLPA experiments. To ensure that the death after ultraviolet irradiation55. Thus, positive selection on genotypes were completely accurate, in most potential discrepancies both these traits could explain the inversion increase in Africa. Inci- breakpoints of each orientation were tested. dentally, the inverted orientation has been linked to higher risk of For inversion genotyping in chimpanzees and gorillas, first we ran the same MLPA and iMLPA assays as described above, and those inversions that did not glaucoma in Europeans (Supplementary Table 4) and glaucoma is work were genotyped by PCR or iPCR. In some cases, new chimpanzee or gorilla more common and severe in individuals from African ancestry56. specific primers and restriction enzymes for iPCR were used to overcome human Other interesting candidates include HsInv0031, HsInv0059, assay problems26 (Supplementary Data 5; Supplementary Data 11). However, this HsInv0124, HsInv0340, and HsInv0389 (Fig. 2). was not always possible due mainly to deletions or genome gaps, and a few inversions in one or both species could not be tested. All polymorphic inversions in Inversions differ from other genetic variants because of their chimpanzees or gorillas were validated by PCR or iPCR of at least one breakpoint expected negative consequences in fertility resulting from the to make sure that there were no errors in the iMLPA results. generation of unbalanced gametes by recombination30,31, which is exemplified by the reduction in frequency with genetic length Analysis of nucleotide variants associated with inversions. We measured (Fig. 4b) and the small number of inversions described compared pairwise LD (r2) between inversions and overlapping and neighboring biallelic to CNVs19. According to this, there could be a maximum length variants up to 200 kb at each side of the inversion from 1000GP Ph3 (including for an inversion to behave neutrally in terms of its fertility effects. SNPs and indels for 434 unrelated individuals) and HapMap release 27 (including fewer SNPs but all the 480 unrelated individuals) using either plink v1.9057 or Above that size, some type of compensatory selection, perhaps Haploview v4.158. Variants located within the breakpoint interval and associated related to advantageous regulatory changes on nearby genes, deletions or IRs were excluded for this and subsequent analysis to avoid possible would be necessary for the inversions to reach a certain genotyping errors. Supplementary Data 12 lists the inversion tag variants with r2 ≥ 10 NATURE COMMUNICATIONS | (2019)10:4222 | https://doi.org/10.1038/s41467-019-12173-x | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-019-12173-x ARTICLE 0.8 from 1000GP and HapMap data considering all the individuals together, as well non-recombining flanking region (up to 20 kb), excluding breakpoint intervals and as each population and population group. This analysis allowed us to detect a few associated IRs or indels to avoid sequence errors. Confidence intervals of age inversion genotypes that did not match those expected from the tag variants and estimates were calculated by bootstrap sampling the same number of total most of them were confirmed to be genotyping errors by independent PCR or iPCR individuals with replacement 1000 times and using both a constant substitution validation (Supplementary Fig. 2). SNPs and indels around inversion regions were rate and two local substitution rates from the divergence with chimpanzee and further classified directly from the genotype data according to its distribution gorilla (see Supplementary Methods). across orientations in fixed (r2 = 1), shared (unambiguously polymorphic in both orientations) and polymorphic in O1 or O2 chromosomes using in-house perl and bash scripts21,26. To minimize possible genotyping errors, only the most reliable Inversion frequency analyses. To control the effect of the study design ascer- variants according to the 1000GP strict accessibility mask were included. Non- tainment bias in the observed frequency of inversions, we simulated the detection recombining flanking regions were defined from the 1000GP data (which has more and genotyping process in biallelic SNPs from 1000GP Ph3. The process was resolution than HapMap) as the longest sequence outside the breakpoints up to a simulated in two steps: (1) selection of those SNPs for which the alternative allele is maximum distance of 20 kb that: (1) does not contain reliable shared variants present in 1000GP individuals matching the demographic and gender composition compatible with a crossing-over event between orientations; and (2) includes most of the detection panel (nine individuals for 38 autosomal or chr. X inversions of the fixed variants (Supplementary Table 2). detected from the fosmid PEM data12 and one individual for six inversions detected exclusively from the genome assembly comparison10); and (2) for each of the PEM inversions, generate a random sample of 10,000 SNPs from the total pool Phasing and visualization of haplotypes. Each orientation haplotypes were of polymorphic SNPs in the PEM panel according to the simulated detection determined following two complementary strategies to minimize errors and obtain probability calculated from the SNP frequency and the inversion characteristics more robust results. First, after testing several commonly-used phasing programs, (see Supplementary Methods). Mean and median frequencies of inversions and we selected PHASE 2.159 because it avoids switch errors in inversion heterozygotes SNPs in the 434 1000GP Ph3 individuals were compared by sampling 10,000 sets of by fixing the phase of the two orientation alleles added at the breakpoint SNPs without replacement, with one matched-SNP per inversion at a time, and positions21,26. Phasing was done independently for the YRI, LWK, EUR, SAS, and empirical P-values were estimated as twice the fraction of samples with values more EAS populations or population groups, using the available trio information when extreme or equal than the observed. possible and five iterations (−x5) for HapMap data and two iterations (−x2) with Different genomic variables were tested to determine their effect on inversion the hybrid model (−MQ) for 1000GP data. Only variants within the inverted frequency: (1) physical length of the inverted region; (2) inversion genetic length; region plus 20 kb of flanking sequence for 1000GP Ph1 data (which was the only (3) number of genes within the inversion or breakpoint regions; (4) distance to the one available at that time) or 200 kb for HapMap data were selected. Second, we closest coding gene; (5) number of mammalian constrained sites inside the took advantage of the 1000GP Ph3 phased haplotypes3 to impute directly the inversion;65 (6) direct functional effect of inversions on genes (Supplementary inversion orientation based on the presence of perfect tag variants (r2 = 1). For Table 3); and (7) size of breakpoint IRs. Inversion genetic length was estimated as inversions without perfect tag variants, only homozygotes or hemizygotes for each the cumulative 4Ner for all the genotyped chromosomes using the 1000GP Ph3 orientation, in which the inversion status of the haplotypes can be assigned SNP data and the LDhat v2.2 rhomap function66. Due to high correlation between unambiguously, were analyzed. several variables, to keep only those with significant regression coefficients and The relationships between the different haplotypes of each dataset were reduce the effect of potential outliers in reduced samples, we built robust regression visualized by combining several data sources and representation methods. Phased models of the autosomal and chr. X inversion frequency in each population group 1000GP Ph1 and HapMap haplotypes were used to build Median-Joining (MJ) by stepwise forward selection of predictors with the lmrob function from networks60 with the NETWORK v.4.6.1.3 software (www.fluxus-engineering.com). robustbase R package67. In addition, we devised our own representation of the haplotype relationships and the distribution of nucleotide changes along the sequence, named integrated haplotype plot (iHPlot), which combines a hierarchical clustering, distance matrix Inversion selection tests. Frequency differences between populations were cal- and visual alignment of the alleles in each polymorphic position, plus the haplotype culated with vcftools (v0.1.15) using the FST statistic68. FST values were compared orientation and the populations in which it is present (see Supplementary Methods with empirical null distributions in the same individuals obtained from bialellic for details). SNPs accessible according to the strict mask, with a defined ancestral allele, and Inversion origin was estimated from the information of the different phasing that have similar frequency and chromosome type (autosome or chr. X) as the and haplotype visualization strategies, which overall showed consistent results inversion (see Supplementary Methods). (see Supplementary Methods). For inversions with perfect tag variants, the analysis LSFS tests are a newly developed family of neutrality tests especially appropriate was based mainly on 1000GP Ph3 haplotypes (including when possible the flanking for inversions, which are a direct application of nearly optimal linear tests for non-recombining region), which allowed a better discrimination of haplotypes, neutrality69. The summary statistic was the frequency spectrum of variants closely filtering of accessible SNPs according to the pilot accessibility mask and had less linked to the inversion, including their linkage pattern (nested or disjoint) with the phasing errors due to the use of sequences from more individuals3. For inversions inverted allele70, and we tested strong positive or balancing selection coefficients. without tag variants, we relied mainly in the phased 1000GP Ph1 haplotypes LSFS was calculated from biallelic 1000GP Ph3 SNPs in the genotyped individuals, iHPlots, since all the genotyped individuals in common could be analyzed. after removing those with a GERP score65 higher than 2 and within 0.5 Mb of any Inversion recurrence events were conservatively estimated by identifying clusters of of the inversions in our dataset. Only 18 autosomal inversions unambiguously haplotypes with both orientations that differ significantly from all others with the phased into the 1000GP Ph3 haplotypes with perfect tag SNPs within 20 kb from same orientation as the potential recurrence event, after eliminating possible the breakpoints were analyzed using 3 kb non-overlapping windows localized phasing errors in inversion heterozygotes (see Supplementary Methods). Around within the inverted or non-recombining flanking regions (skipping the 1–25 individuals supporting each of the independent recurrence events were breakpoints, IRs and indels to avoid genotype errors). Inversion windows were validated by PCR or iPCR (in this case testing mostly both breakpoints) to discard compared against the empirical LSFS computed around all autosomal SNPs and possible inversion genotyping errors (Supplementary Data 3). Also, we validated tests were conditioned on the inversion frequency in the different populations. the genotype of many more individuals with unexpected orientation-haplotype Each population and window was tested separately and population P values of the combinations and of other inversions with a high proportion of shared SNPs for same windows were combined via Edgington’s method71, whereas the results which recurrence events could not be clearly identified. HsInv0832 inversion rate across different windows of an inversion were combined using a conservative and was estimated from the publicly available information of 232 of the 282 genotyped an approximate approach (see Supplementary Methods). males (Supplementary Data 4) by calculating the total number of mutations and NCD1 and NCD2 statistics44 to test long-term balancing selection acting on generations along the genealogical tree of the analyzed Y chromosomes using the autosomal and chr. X inversion regions were computed for three target frequencies approach of Repping et al.61 and Hallast et al.62 (see Supplementary Methods). (0.3, 0.4, and 0.5) in overlapping windows of 2 kb (with 1 kb step), defined with the same criteria as in the LSFS test, using 1000GP Ph3 SNPs (accessible according to the pilot mask) from all individuals of the seven studied populations. Only Ancestral orientation and inversion age estimate. Besides the experimental windows with a minimum of eight informative sites (either human polymorphisms analysis in chimpanzees and gorillas, the ancestral state of inversions was com- or fixed differences with chimpanzee in NCD2 only) and at least 16.7% of positions plemented by bioinformatic and manual inspection of four of the best non-human covered by hg19-panTro4 alignments were considered. Finally, an empirical P- primate genome assemblies (see Supplementary Methods). Due to their fragmented value was assigned for each inversion region and population by comparing the tests status, the orientation of several available ancient hominin genomes (Supple- results with a null genome-wide distribution obtained by sampling regions of the mentary Data 6) was determined by identifying the reads that span the O1 or O2 same size as the inversion (see Supplementary Methods). breakpoints of 19 NH inversions without IRs using a library with four 100-bp sequences centered at the two breakpoints of both orientations (Supplementary Data 13) and a slightly modified version of the BreakSeq pipeline21,27,28. Gene-expression analysis in LCLs. We analyzed 42 inversions (excluding two Age of unique inversions was obtained with the usual divergence-based with MAF < 1% and HsInv0832 in chr. Y) in 173 experimentally genotyped indi- approach63,64, using the pairwise differences between orientations from all SNPs viduals (42 CEU, 84 TSI, and 47 YRI) with GEUVADIS45 and 1000GP Ph33 data. (excluding indels) in the available 1000GP Ph3 haplotypes and the largest of the Besides, we imputed 33 inversions in the complete set of 445 GEUVADIS indi- two average pairwise differences within O1 or O2 sequences. To have more viduals in common with 1000GP Ph3 (89 CEU, 91 TSI, 86 GBR, 92 FIN, and 87 information in short inversions, we considered the inverted region and any extra YRI) using the already identified perfect tag SNPs (r2 = 1) (19 inversions) or NATURE COMMUNICATIONS | (2019)10:4222 | https://doi.org/10.1038/s41467-019-12173-x | www.nature.com/naturecommunications 11 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-019-12173-x IMPUTE v2.3.272 (14 inversions with >90% average imputation accuracy) (Sup- References plementary Fig. 4A). LCL raw RNA-Seq reads (ArrayExpress experiment E-GEUV- 1. The International HapMap 3 Consortium. Integrating common and rare 1) were aligned against the GRCh38.p10 human genome with STAR v2.4.2a73 genetic variation in diverse human populations. Nature 467, 52–58 (2010). using GENCODE version 26 annotations74. Gene-expression levels were estimated 2. The 1000 Genomes Project Consortium. An integrated map of genetic as reads per kilobase per million mapped reads (RPKM) and transcript expression variation from 1092 human genomes. Nature 491, 56–65 (2012). levels were quantified with RSEM v1.2.3175. cis-eQTL analysis was done through 3. The 1000 Genomes Project Consortium. A global reference for human genetic linear regressions implemented in QTLtools v1.176, considering 850 genes and variation. Nature 526, 68–74 (2015). 3318 transcripts expressed in at least 20% of the samples and with TSS within 1 Mb 4. Sudmant, P. H. et al. An integrated map of structural variation in 2504 human from inversions. First, we carried out a targeted study to test only the association genomes. Nature 526, 75–81 (2015). with the genotypes of each inversion. Second, we performed a joint eQTL analysis 5. Walter, K. et al. The UK10K project identifies rare variants in health and with inversions and neighboring 1000GP variants to estimate their contribution to disease. Nature 526, 82–90 (2015). observed gene-expression changes and identify lead eQTLs. Expression values were 6. Manolio, T. A. et al. Finding the missing heritability of complex diseases. adjusted by gender, the first three principal components from 1000GP data (cor- Nature 461, 747–753 (2009). responding to continent, population and population structure) and a set of 7. Eichler, E. E. et al. Missing heritability and strategies for finding the expression principal components to capture technical confounding factors (for underlying causes of complex disease. Nat. Rev. Genet. 11, 446–450 (2010). genes and transcripts, respectively, 10 and 20 in the experimental dataset and 35 8. Alkan, C., Coe, B. P. & Eichler, E. E. Genome structural variation discovery and 50 in the imputed dataset) (see Supplementary Methods). Next, they were and genotyping. Nat. Rev. Genet. 12, 363–376 (2011). transformed to match normal distributions N(0,1) to avoid false positive associa- 9. Puig, M., Casillas, S., Villatoro, S. & Cáceres, M. Human inversions and their tions due to outliers and significance was established at 5% false discovery functional consequences. Brief. Funct. Genom. 14, 369–379 (2015). rate (FDR). 10. Levy, S. et al. The diploid genome sequence of an individual human. PLoS Biol. 5, e254 (2007). Inversions as eQTLs in other tissues and conditions. Gene-expression effects in 11. Catacchio, C. R. et al. Inversion variants in human and primate genomes. other tissues for 26 inversions were examined by looking whether their highest Genome Res. 28, 910–920 (2018). associated SNPs across all populations (r2 ≥ 0.8) have been identified as cis-eQTLs 12. Kidd, J. M. et al. Mapping and sequencing of structural variation from eight in different human tissues by the GTEx project (GTEx Analysis Release v7)46. Also, human genomes. Nature 453, 56–64 (2008). we extended the analysis to seven recurrent inversions with SNPs in moderate LD 13. Huddleston, J. et al. Discovery and genotyping of structural variation from (r2 ≥ 0.6). To determine the potential causal variants, we checked if those SNPs in long-read haploid genome sequence data. Genome Res. 27, 677–685 (2017). highest LD with the inversions were reported as being the first or second lead eQTL 14. Audano, P. A. et al. Characterizing the major structural variant alleles of the in a specific tissue. The same strategy was applied to link inversions to immune human genome. Cell 176, 663–675 (2019). eQTLs associated to the transcriptional response to Listeria and Salmonella of 15. Shao, H. et al. npInv: accurate detection and genotyping of inversions using primary macrophages from African-American (n = 76) and European-American long read sub-alignment. BMC Bioinforma. 19, 261 (2018). (n = 99) individuals52 and of macrophages differentiated from 123 induced plur- 16. Sanders, A. D. et al. Characterizing polymorphic inversions in human ipotent stem cell (iPSC) lines of European origin to Salmonella plus interferon γ genomes by single cell sequencing. Genome Res. 26, 1575–1587 (2016). stimulation53. 17. Li, L. et al. OMSV enables accurate and comprehensive identification of large structural variations from nanochannel-based single-molecule optical maps. Association of inversions with GWAS SNPs. Association of inversions with Genome Biol. 18, 230 (2017). specific traits or diseases was based on the NHGRI Catalog of published GWAS 18. Chaisson, M. J. P. et al. Multi-platform discovery of haplotype-resolved (release 2017-07-17)48 and the GWASdb (release 2015-08-19)49 databases. To structural variation in human genomes. Nat. Commun. 10, 1784 (2019). remove redundant entries, the strongest signal per locus (±100 kb genomic region) 19. Martínez-Fundichely, A. et al. InvFEST, a database integrating information of was selected. Only inversions in high LD (r2 ≥ 0.8) with a GWAS signal in the polymorphic inversions in the human genome. Nucleic Acids Res. 42, studied population or the closest one available were considered (Supplementary D1027–D1032 (2014). Table 4). To investigate if the number of GWAS signals in the inversion and 20. Hehir-Kwa, J. Y. et al. A high-quality human reference panel reveals the flanking regions (±20 kb) was higher than expected by chance, we crossed GWAS complexity and distribution of genomic structural variants. Nat. Commun. 7, Catalog and GWASdb signals with 1000GP variants and grouped together those in 12989 (2016). high LD (r2 ≥ 0.8) and associated to the same phenotype. Enrichment P values were 21. Vicente-Salvador, D. et al. Detailed analysis of inversions predicted between calculated by comparison with a null distribution from 1000 random genomic two human genomes: errors, real polymorphisms, and their origin and regions as a background model for each inversion, controlling by inversion size and population distribution. Hum. Mol. Genet. 26, 567–581 (2017). SNP frequencies (average SNP frequency ±0.01 and Chi-square test P > 0.05 for the 22. Stefansson, H. et al. A common inversion under selection in Europeans. Nat. number of SNPs with MAF <0.2 and ≥0.2 compared to the inversion region) and Genet. 37, 129–137 (2005). excluding gaps and chr. Y. To test the enrichment in specific inversions, we 23. Salm, M. P. A. et al. The origin, global distribution, and functional impact of repeated the same analysis computing a one-tailed empirical P-value for each inversion. the human 8p23 inversion polymorphism. Genome Res. 22, 1144–1153 (2012). 24. Antonacci, F. et al. Characterization of six human disease-associated inversion polymorphisms. Hum. Mol. Genet. 18, 2555–2566 (2009). Statistical information. Details of the statistical tests are described in the corre- 25. Pang, A. W. C., Migita, O., Macdonald, J. R., Feuk, L. & Scherer, S. W. sponding sections. In many cases P-values were derived from empirical genome- Mechanisms of formation of structural variation in a fully sequenced human wide null distributions from at least 1000 random samples and two-tailed tests genome. Hum. Mutat. 34, 345–354 (2013). were always used, except when specifically mentioned. 26. Aguado, C. et al. Validation and genotyping of multiple human polymorphic inversions mediated by inverted repeats reveals a high degree of recurrence. Reporting summary. Further information on research design is available in PLoS Genet. 10, e1004208 (2014). the Nature Research Reporting Summary linked to this article. 27. Lucas-Lledó, J. I., Vicente-Salvador, D., Aguado, C. & Cáceres, M. Population genetic analysis of bi-allelic structural variants from low-coverage sequence data with an expectation-maximization algorithm. BMC Bioinforma. 15, 163 (2014). Data availability 28. Puig, M. et al. Functional impact and evolution of a novel human All data described in this article are available in the Supplementary Information and in polymorphic inversion that disrupts a gene and creates a fusion transcript. the InvFEST database (http://invfestdb.uab.cat/). In addition, inversion genotypes have PLoS Genet. 11, e1005495 (2015). been deposited in the dbVar database (https://www.ncbi.nlm.nih.gov/dbvar/) under 29. González, J. R. et al. A common 16p11.2 inversion underlies the joint accession number nstd169 [https://www.ncbi.nlm.nih.gov/dbvar/studies/nstd169/]. The susceptibility to asthma and obesity. Am. J. Hum. Genet. 94, 361–372 (2014). source data underlying Figs. 2, 3a-c, 4a-c, 5a, c, e and 6a-b and Supplementary Figs. 1, 2a- 30. Hoffmann, A. A. & Rieseberg, L. H. Revisiting the impact of inversions in b, 4a, f, 5 and 8a-c are provided as a Source Data file. All other relevant data are available evolution: from population genetic markers to drivers of adaptive shifts and upon request. speciation? Annu. Rev. Ecol. Evol. Syst. 39, 21–42 (2008). 31. Kirkpatrick, M. How and why chromosome inversions evolve. PLoS Biol. 8, Code availability e1000501 (2010). Computer code is available upon request. 32. Wellenreuther, M. & Bernatchez, L. Eco-evolutionary genomics of chromosomal inversions. Trends Ecol. Evol. 33, 427–440 (2018). 33. Imsland, F. et al. The Rose-comb mutation in chickens constitutes a structural Received: 21 December 2018 Accepted: 27 August 2019 rearrangement causing both altered comb morphology and defective sperm motility. PLoS Genet. 8, e1002775 (2012). 12 NATURE COMMUNICATIONS | (2019)10:4222 | https://doi.org/10.1038/s41467-019-12173-x | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-019-12173-x ARTICLE 34. Lakich, D., Kazazian, H. H., Antonarakis, S. E. & Gitschier, J. Inversions 67. Maechler, M. et al. robustbase: Basic Robust Statistics R package version 0.93- disrupting the factor VIII gene are a common cause of severe haemophilia A. 2. (2018). Nat. Genet. 5, 236–241 (1993). 68. Weir, B. S. & Cockerham, C. C. Estimating F-statistics for the analysis of 35. Myers, A. J. et al. The H1c haplotype at the MAPT locus is associated with population structure. Evolution 38, 1358–1370 (1984). Alzheimer’s disease. Hum. Mol. Genet. 14, 2399–2404 (2005). 69. Ferretti, L., Perez-Enciso, M. & Ramos-Onsins, S. Optimal neutrality tests 36. Zabetian, C. P. et al. Association analysis of MAPT H1 haplotype and based on the frequency spectrum. Genetics 186, 353–365 (2010). subhaplotypes in Parkinson’s disease. Ann. Neurol. 62, 137–144 (2007). 70. Ferretti, L. et al. The neutral frequency spectrum of linked sites. Theor. Popul. 37. Webb, A. et al. Role of the tau gene region chromosome inversion in Biol. 123, 70–79 (2018). progressive supranuclear palsy, corticobasal degeneration, and related 71. Edgington, E. S. An additive method for combining probability values from disorders. Arch. Neurol. 65, 1473–1478 (2008). independent experiments. J. Psychol. 80, 351–363 (1972). 38. Okbay, A. et al. Genetic variants associated with subjective well-being, 72. Howie, B. N., Donnelly, P. & Marchini, J. A flexible and accurate genotype depressive symptoms and neuroticism identified through genome-wide imputation method for the next generation of genome-wide association analyses. Nat. Genet. 48, 624–633 (2016). studies. PLoS Genet. 5, e1000529 (2009). 39. de Jong, S. et al. Common inversion polymorphism at 17q21.31 affects expression 73. Dobin, A. et al. STAR: ultrafast universal RNA-seq aligner. Bioinformatics 29, of multiple genes in tissue-specific manner. BMC Genom. 13, 458 (2012). 15–21 (2013). 40. Chiang, C. et al. The impact of structural variation on human gene expression. 74. Harrow, J. et al. GENCODE: The reference human genome annotation for Nat. Genet. 49, 692–699 (2017). The ENCODE Project. Genome Res. 22, 1760–1774 (2012). 41. Kehr, B. et al. Diversity in non-repetitive human sequences not found in the 75. Li, B. & Dewey, C. N. RSEM: accurate transcript quantification from RNA-Seq reference genome. Nat. Genet. 49, 588–593 (2017). data with or without a reference genome. BMC Bioinforma. 12, 323 (2011). 42. Schouten, J. P. et al. Relative quantification of 40 nucleic acid sequences by 76. Delaneau, O. et al. A complete tool set for molecular QTL discovery and multiplex ligation-dependent probe amplification. Nucleic Acids Res. 30, e57 analysis. Nat. Commun. 8, 15452 (2017). (2002). 77. Poznik, G. D. et al. Punctuated bursts in human male demography inferred 43. Prüfer, K. et al. The complete genome sequence of a Neanderthal from the from 1,244 worldwide Y-chromosome sequences. Nat. Genet. 48, 593–599 Altai Mountains. Nature 505, 43–49 (2014). (2016). 44. Bitarello, B. D. et al. Signatures of long-term balancing selection in human 78. Pérez-Palma, E. et al. Overrepresentation of glutamate signaling in genomes. Genome Biol. Evol. 10, 939–955 (2018). Alzheimer’s disease: network-based pathway enrichment using meta-analysis 45. Lappalainen, T. et al. Transcriptome and genome sequencing uncovers of genome-wide association studies. PLoS One 9, e95413 (2014). functional variation in humans. Nature 501, 506–511 (2013). 79. Iyer, M. K. et al. The landscape of long noncoding RNAs in the human 46. GTEx Consortium. Genetic effects on gene expression across human tissues. transcriptome. Nat. Genet. 47, 199–208 (2015). Nature 550, 204–213 (2017). 47. Horton, R. et al. Variation analysis and gene annotation of eight MHC haplotypes: the MHC Haplotype Project. Immunogenetics 60, 1–18 (2008). Acknowledgements 48. MacArthur, J. et al. The new NHGRI-EBI Catalog of published genome-wide We thank Cristina Aguado, Olga Dolgova, Teresa Soos, Esteban Urrea, David Vicente- association studies (GWAS Catalog). Nucleic Acids Res. 45, D896–D901 (2017). Salvador, and Roser Zaurín for help with inversion genotyping, Antonio Barbadilla, Ruth 49. Li, M. J. et al. GWASdbv2: an update database for human genetic variants Gómez, José Ignacio Lucas-Lledó, Sebastián Ramos-Onsins, and Alfredo Ruiz for help identified by genome-wide association studies. Nucleic Acids Res. 44, with the evolutionary analysis, Mariona Bellet, Robert Castelo, Diego Garrido, and D869–D876 (2016). Roderic Guigó for help with the gene-expression analysis, Sònia Casillas and Alexander 50. Handsaker, R. E. et al. Large multiallelic copy number variations in humans. Martínez-Fundichely for help with the inversion selection, and Xavier Estivill, Tomàs Nat. Genet. 47, 296–303 (2015). Marquès-Bonet, Aurora Ruiz-Herrera, the Coriell Institute for Medical Research, the 51. Gymrek, M. et al. Abundant contribution of short tandem repeats to gene Barcelona Zoo and the Banc de Teixits Animals de Catalunya (BTAC) for providing the expression variation in humans. Nat. Genet. 48, 22–29 (2016). human and non-human primate samples used in this study. This work was supported by 52. Nédélec, Y. et al. Genetic ancestry and natural selection drive population research grants ERC Starting Grant 243212 (INVFEST) from the European Research differences in immune responses to pathogens. Cell 167, 657–669 (2016). Council under the European Union Seventh Research Framework Programme (FP7), 53. Alasoo, K. et al. Shared genetic effects on chromatin and gene expression BFU2013-42649-P and BFU2016-77244-R funded by the Agencia Estatal de Investiga- indicate a role for enhancer priming in immune response. Nat. Genet. 50, ción (AEI, Spain) and the European Regional Development Fund (FEDER, EU), and 424–431 (2018). 2014-SGR-1346 and 2017-SGR-1379 from the Generalitat de Catalunya (Spain) to M.C., 54. Sun, B. et al. Genomic atlas of the human plasma proteome. Nature 558, a PIF PhD fellowship from the Universitat Autònoma de Barcelona (Spain) to C.G.D., a 73–79 (2018). La Caixa Doctoral fellowship to J.L.J., and a FPI PhD fellowship from the Ministerio de 55. Lee, J. Y. W. et al. Large intragenic deletion in DSTYK underlies autosomal- Economía y Competitividad (Spain) to M.O. and I.N. M.G.V. was supported in part by recessive complicated spastic paraparesis, SPG23. Am. J. Hum. Genet. 100, POCI-01-0145-FEDER-006821 funded through the Operational Programme for Com- 364–370 (2017). petitiveness Factors (COMPETE, EU) and UID/BIA/50027/2013 from the Foundation 56. Tielsch, J. M. et al. Racial variations in the prevalence of primary open-angle for Science and Technology (FCT, Portugal). glaucoma. The Baltimore Eye Survey. JAMA 266, 369–374 (1991). 57. Purcell, S. et al. PLINK: a tool set for whole-genome association and Author contributions population-based linkage analyses. Am. J. Hum. Genet. 81, 559–575 (2007). M.C. conceived the inversion genotyping strategy, devised the study and oversaw all the 58. Barrett, J. C., Fry, B., Maller, J. & Daly, M. J. Haploview: analysis and steps; S.V., M.P., and M.C. designed the genotyping assays; S.V., D.I., A.D., and M.P. visualization of LD and haplotype maps. Bioinformatics 21, 263–265 (2005). carried out experiments; C.G.D., M.G.V., I.O., C.L.F., M.P., and M.C. analyzed evolu- 59. Stephens, M. & Donnelly, P. A comparison of bayesian methods for haplotype tionary data; C.G.D., J.L.J., D.C., B.B., I.N., P.O., A.A., and L.F. performed selection tests; reconstruction from population genotype data. Am. J. Hum. Genet. 73, J.L.J., M.O., L.P., T.E., and M.C. analyzed functional effects; A.B. provided samples; M.C., 1162–1169 (2003). C.G.D., J.L.J., M.G.V., L.F., and M.P. wrote the paper and all the authors contributed 60. Bandelt, H. J., Forster, P. & Röhl, A. Median-joining networks for inferring comments to the final version of the manuscript. intraspecific phylogenies. Mol. Biol. Evol. 16, 37–48 (1999). 61. Repping, S. et al. High mutation rates have driven extensive structural polymorphism among human Y chromosomes. Nat. Genet. 38, 463–467 (2006). Additional information 62. Hallast, P., Balaresque, P., Bowden, G. R., Ballereau, S. & Jobling, M. A. Supplementary Information accompanies this paper at https://doi.org/10.1038/s41467- Recombination dynamics of a human Y-chromosomal palindrome: rapid GC- 019-12173-x. biased gene conversion, multi-kilobase conversion tracts, and rare inversions. PLoS Genet. 9, e1003666 (2013). Competing interests: The authors declare no competing interests. 63. Hasson, E. & Eanes, W. F. Contrasting histories of three gene regions associated with In(3L)Payne of Drosophila melanogaster. Genetics 144, Reprints and permission information is available online at http://npg.nature.com/ 1565–1575 (1996). reprintsandpermissions/ 64. Corbett-Detig, R. B. & Hartl, D. L. Population genomics of inversion polymorphisms in Drosophila melanogaster. PLoS Genet. 8, e1003056 (2012). Peer review information Nature Communications thanks Megan Dennis and the other, 65. Davydov, E. V. et al. Identifying a high fraction of the human genome to be anonymous, reviewer(s) for their contribution to the peer review of this work. Peer under selective constraint using GERP++. PLoS Comput. Biol. 6, e1001025 reviewer reports are available. (2010). 66. Auton, A. & McVean, G. Recombination rate estimation in the presence of Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in hotspots. Genome Res. 17, 1219–1227 (2007). published maps and institutional affiliations. NATURE COMMUNICATIONS | (2019)10:4222 | https://doi.org/10.1038/s41467-019-12173-x | www.nature.com/naturecommunications 13 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-019-12173-x Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/ licenses/by/4.0/. © The Author(s) 2019 14 NATURE COMMUNICATIONS | (2019)10:4222 | https://doi.org/10.1038/s41467-019-12173-x | www.nature.com/naturecommunications