Young inversion with multiple linked QTLs under selection in a hybrid zone

Fixed chromosomal inversions can reduce gene flow and promote speciation in two ways: by suppressing recombination and by carrying locally favoured alleles at multiple loci. However, it is unknown whether favoured mutations slowly accumulate on older inversions or if young inversions spread because they capture pre-existing adaptive quantitative trait loci (QTLs). By genetic mapping, chromosome painting and genome sequencing, we have identified a major inversion controlling ecologically important traits in Boechera stricta. The inversion arose since the last glaciation and subsequently reached local high frequency in a hybrid speciation zone. Furthermore, the inversion shows signs of positive directional selection. To test whether the inversion could have captured existing, linked QTLs, we crossed standard, collinear haplotypes from the hybrid zone and found multiple linked phenology QTLs within the inversion region. These findings provide the first direct evidence that linked, locally adapted QTLs may be captured by young inversions during incipient speciation.

C hromosome inversions play an important role in local adaptation and speciation 1,2 , and selectively important inversions have been identified in many species 3,4 . Selection due to different environmental factors or stages in the life cycle 1 may favour inversions carrying locally adapted alleles at several loci. In addition, established inversions are predicted to accumulate selectively important genetic differences, which may contribute to reproductive isolation during speciation 1 .
Although few studies have identified the actual loci that influence selection on inversions 2,4,5 , rearrangements may be favoured due to gene alterations near breakpoints 6 , chromatin changes 7 or combinations of advantageous, co-adapted alleles 8 . Inversions suppress recombination, so locally advantageous alleles may segregate together, causing higher fitness than recombinant haplotypes 9 . Most evolutionary studies have focused on widespread, older inversions, so we have little knowledge of the evolutionary processes that guide their initial increase in frequency. For example, do inversions drift to higher frequency, and then acquire new advantageous mutations after they are common? Or are multiple linked, advantageous alleles captured in a new inversion, allowing them to spread together? Analysis of younger inversions may elucidate the evolutionary forces controlling the initial spread of chromosome inversions, which therefore influence their role in adaptation and speciation 4,8 .
Related species often differ for chromosome inversions that carry locally favoured alleles at multiple loci 10,11 . A key distinction among models for the evolution of inversions is whether early frequency increase is due to genetic drift or natural selection. Genetic drift might predominate initially, with subsequent accumulation of advantageous variants 12 . Alternatively, the Kirkpatrick-Barton model 9 argues that linked, locally adapted alleles exist first, and subsequently are captured within a new, selectively favoured inversion 13 . In this 'inversion-late' evolutionary sequence 1,5 , linked quantitative trait loci (QTLs), similar to the ancestral haplotype that gave rise to the inversion, may still exist in non-inverted genotypes 9 . Here, we test these predictions of the Kirkpatrick-Barton model. First, we introduce ecologically diverged subspecies of Boechera stricta. Next, we examine a young inversion to infer the selective forces controlling its early increase in frequency. Finally, we cross collinear, standard genotypes from the hybrid zone to ask whether old, linked QTLs can be found within the inversion region.
the most similar Bsi1-inv and Bsi1-std haplotypes ( Supplementary  Fig. 5) are located only 1.7 km apart (Fig. 3b, left). Using measured mutation rates 24 and assuming two years per generation, the inversion dates to ~2.1-8.8 ka (Methods). Because these widely distributed subspecies are ecologically diverged 16 , while the Bsi1 inversion originated very recently, this polymorphism is compatible with the inversion-late Kirkpatrick-Barton model.
Previous analysis of a West-inv × East-std mapping population grown in the inversion zone found many QTLs controlling components of fitness 25 and strong selection on flowering time 17 , influenced by recent climate warming 26 . Native West Bsi1-inv alleles had higher fitness than East Bsi1-std alleles 25 . The inversion (spanning ~4% of the genome) explains 7.5% of heritable variation for lifetime fecundity, and inversion homozygotes differ by 0.56 genetic standard deviations for this trait (t = 3.64, d.f. = 165, P = 0.0004). Here, we focused on the effects of the inversion by examining flowering time in two F7 near-isogenic line (NIL) families from this cross. We found no segregation distortion or deviation from Mendelian ratios (Supplementary Table 2) in either family, and the inversion had significant effects on flowering time in these NIL families in the greenhouse (F-ratio, F 2, 282 = 21.81; P < 0.001), where the Bsi1-std haplotype flowered ~2.0 days faster than the Bsi1-inv haplotype ( Supplementary  Fig. 6). Thus, the inversion controls 40% of the average 5.0 day difference in flowering time between subspecies (Supplementary Table 2d). While this pedigree cannot resolve QTLs within the inversion, this cross shows that the inversion has significant effects on an ecologically important trait that experiences natural selection 17 . Gene(s) within the Bsi1 inversion control flowering in the field, and thus may contribute to reproductive isolation during speciation 27 .
To test for phenotypic effects of the inversion or transcription changes near the breakpoints, we compared closely related inverted and standard haplotypes from the inversion zone, using a sympatric West-inv × West-std F2 cross (Fig. 3b, left). We found high seed germination (> 99%), and no Bsi1 segregation distortion or deviation from Mendelian ratios (Supplementary Table 3). Multivariate analysis of variance (MANOVA) showed a significant effect of Bsi1 on a suite of phenological and morphological traits divergence between these subspecies 15 , with West genotypes growing in sites with more constant and abundant water supply. In and near the hybrid zone, the East and West subspecies show significant ecological differentiation across local environmental gradients, and the geographic distribution of the inversion falls within the typical range of hybrid zone habitats. The West subspecies has a significantly faster growth rate, larger leaf area, less succulent leaves, delayed reproductive time and longer flowering duration 16 , and West × East crosses found QTLs for flowering traits, leaf number, defensive chemistry, herbivore resistance, cold tolerance, overwinter survival, fecundity and lifetime fitness [17][18][19][20] . In addition, Q ST -F ST analysis showed that phenology and some morphological traits have experienced divergent selection between subspecies 16 . In our initial cross 17 , genetic mapping in West × East recombinant inbred lines (RILs) identified a region of suppressed recombination (Supplementary Fig. 1) on linkage group 1 (LG1). Here, we combined short-read sequencing, end-sequencing of bacterial artificial chromosomes (BACs), linkage mapping and physical mapping by Whole Genome Profiling (WGP) to assemble a high-quality genomic sequence ( Supplementary Fig. 2) that enabled evolutionary analysis of the inversion. Chromosome painting ( Fig. 1 and Supplementary  Fig. 3) cytogenetically verified the presence of an 8.4 Mb paracentric inversion in the West genotype (Bsi1; Fig. 2a and Supplementary  Fig. 4), spanning ~31 cM (~4% of the genome, the 'inversion region'; Supplementary Table 1). We developed primers to score the inversion using PCR (Fig. 2b), showing that the common, non-inverted allele (Bsi1-std) has the ancestral orientation found in closely related Capsella and Arabidopsis lyrata 21 , and is found in > 200 East and West populations across the species range ( Fig. 3a and Dryad Data Archive-see Data availability section). The derived, inverted allele (Bsi1-inv) was identified in a cluster of West and hybrid populations in the contact zone (Fig. 3b), where it has risen to high frequency. Genotypes carrying the inversion span a narrow geographic range (~14 km, the 'inversion zone'), where pollen and geological analyses 22 Fig. 7). Using replicated F3 homozygotes from this cross, we found no significant expression differences between inversion and standard haplotypes for any of the seven genes flanking the inversion breakpoints (Supplementary Table 5), although power to detect subtle, quantitative differences is limited. Thus, this comparison of West-inv versus West-std Bsi1 haplotypes found genetic differences for ecologically important traits, but little evidence for functional changes at the inversion breakpoints.
Evidence for natural selection. To test whether selection had favoured this inversion, we analysed patterns of population genetic variation ( Fig. 4 and Supplementary Fig. 8). Bsi1-inv genotypes show lower LG1 polymorphism, lower Tajima's D and more linkage disequilibrium (LD) than Bsi1-std individuals. These patterns are compatible with neutral drift in this partially inbreeding species, or with a selective sweep on the inversion 28 . In contrast, young, derived mutations at high frequency suggest positive directional selection 29 , so we asked whether the frequency of Bsi1-inv is typical of the frequencies of all private derived SNP alleles in the same population (the 54 similar genotypes in the left portion of Supplementary  Fig. 9c and the left portion of Fig. 3b). The inversion is at relatively high frequency (0.63) in this group, but it is not fixed. We found 2,416 SNPs that, like the inversion, are confined to this population.
Among these, only 63 SNPs (2.6%) have derived allele frequencies greater than the frequency of the Bsi1-inv inversion. Hence, the frequency of the Bsi1-inv inversion allele is higher than 97.4% of comparable derived alleles in this population-the inversion is a high frequency outlier, supporting the hypothesis of positive directional selection.
QTLs within the inversion. During speciation, co-adapted gene complexes within inversions might reduce gene flow if they preserve favourable combinations of alleles at multiple loci, reducing the frequency of disadvantageous allelic combinations. Although inversions can be engineered for proof of function in some organisms 30 , this is infeasible in Boechera. However, the Kirkpatrick-Barton model predicts that relatives of the ancestral haplotype that gave rise to the inversion might still exist as non-inverted genotypes nearby. To test for such linked QTLs, we crossed collinear West-std × East-std genotypes from the contact zone and tested for QTLs within the inversion region, using freely recombining F3:F4 families. We found several multivariate QTLs altering ecologically important phenology and development traits (Fig. 5a). To clarify differences among these QTLs, we examined differences in their pleiotropic effects (Fig. 5b), using discriminant function analysis (DFA) to find the axis of greatest divergence between the multivariate trait means at each QTL peak. Each DFA axis quantifies the direction of pleiotropy that controls the effects of a locus, and QTLs influencing these composite traits were mapped across the inversion region. We found several distinct QTL peaks, which show divergent patterns of pleiotropy among these linked loci. Finally, we estimated the time of divergence between these parental genotypes (Fig. 5c). These results show that the QTL haplotypes diverged at least 50-100 ka, and are therefore much older than the inversion, which arose less than 10 ka.  DNA extraction for genomic sequencing. Seed sterilization. Seeds of B. stricta ecotypes (LTM and SAD12) were initially surface sterilized with ethanol (95%) for 2 min followed by treatment with 15% sodium hypochlorite solution (with a few drops of Tween-20) for 25 min on a rotatory shaker (Spex CertiPrep; 125 rpm) at room temperature. Surface sterilized seeds were thoroughly washed multiple times with sterile water and were suspended in 0.1% agar. Seeds were kept for stratification in a cold room at 4 °C for 4 d under dark conditions.
Raising of aseptic seedling cultures. To raise aseptic seedlings of LTM and SAD12 ecotypes, surface sterilized and stratified seedlings were inoculated in 250 ml flasks containing 50 ml of sterile half-strength MS liquid medium (pH 5.7) with 0.5 g l −1 MES and 2% sucrose. About 70-80 seeds were inoculated per flask. To obtain etiolated seedlings, the flasks were covered with aluminium foil and kept on a bench top rotatory shaker at 110 rpm for 20 d at room temperature.
Extraction of genomic DNA from etiolated seedlings. Freshly grown 20-day-old etiolated seedlings of both ecotypes were used for extraction of nuclear DNA. Isolation of nuclei was performed according to ref. 18 with slight modifications, as it allowed isolation of clean high-molecular-size nuclear DNA with minimal contamination from organelle DNA. Briefly, about ten flasks of etiolated seedling cultures were removed from the growing medium and thoroughly washed with ice-cold water and placed in ice-cold ethyl ether for 3 min. Subsequently the seedlings were washed five times with ice-cold TE buffer (pH 7.0). The plant material was quickly blotted on sterile filter paper and homogenized in MEB buffer using a commercial blender. The homogenate was filtered through four layers of cheesecloth followed by two layers of mira cloth. Triton X-100 was added to the filtrate at 0.5% concentration and centrifuged to collect the pellet. The pellet was suspended in MPDB solution and gently layered on 37.5% percoll made with MPDB. Nuclei were collected as a pellet after centrifugation, and the pellet was washed by resuspending it in 20 ml of MPBD followed by centrifugation. The nuclei pellet was again suspended in a MPDB solution and high-molecular-weight nuclear DNA was extracted using the Qiagen genomic tip 100/G protocol as per the manufacturer's instructions, with some modifications. The pellet consisting of nuclei was resuspended in the lysis buffer and incubated at 40 °C for 15 min. RNase A and T1 were added to the suspension and further incubated at 37 °C for an additional 30 min. Subsequently, Proteinase K was added to the suspension at a final concentration of 150 μ g ml −1 and incubated for 3 h at 45 °C with gentle shaking. The suspension was cleared by centrifugation at 8,000 rpm for 15 min. The cleared suspension was added to the pre-equilibrated G-100 column and further purification was performed as per the protocol provided by the manufacturer. The eluate consisting of genomic DNA was further precipitated using isopropanol, washed with 70% ethanol and subsequently suspended in TE (pH 8.0) buffer.
Total RNA isolation. To collect RNA from roots, the LTM-genotype seeds were surface sterilized with 12% bleach and grown under aseptic conditions on Whatman filter paper bridges in glass tissue culture tubes containing half-strength MS media with 1% sucrose. Each tube was inoculated with two to three seeds and allowed to grow. The roots obtained from the multiple seedlings were pooled and RNA was isolated using a plant RNeasy kit from Qiagen. For other tissues, we extracted RNA from soil-grown plants at several stages: 15-day-old seedlings, rosette leaves from juvenile plants, cauline leaves from plants aged 2-3 months. Finally, from 7-month-old plants, we extracted RNA from inflorescences, stems, flowers and siliques. For isolation of total RNA from various tissues at different developmental stages, the tissues were snap frozen in liquid nitrogen and stored at − 80 °C. The frozen tissue was homogenized with a pestle and mortar using liquid nitrogen. About 150 mg of homogenized tissue was used for isolation. Total RNA was isolated using a plant RNeasy kit. Each of the independent total RNA preparations was subjected to DNase treatment to remove traces of contaminating genomic DNA using RNase-free DNase (Promega). For each tissue type, three independent isolations were pooled together and used for further analysis.
Sequencing, assembly and annotation. Eighteen genomic DNA libraries and one transcriptome complementary DNA (cDNA) library were prepared (Dryad Data Archive) for sequencing using Illumina and Roche 454. The genome assembly of B. stricta employed Meraculous (May 2013 build) 31 , as described on Phytozome (https://phytozome.jgi.doe.gov/pz/portal.html#!info?alias= Org_Bstricta): contigs were assembled using a size of 51 oligomers and a minimum depth of 8 from 2 × 150 Illumina HiSeq reads from an unamplified 250 bp whole genome shotgun fragment library. Several discrete rounds of scaffolding were performed with libraries containing inserts ranging in size from 250 bp to 40 kb (sequenced with various Illumina technologies). At each round of scaffolding the minimum pairing threshold was chosen by exhaustive optimization of N50 scaffold length. Assembly gaps were closed with 2 × 150 Illumina HiSeq reads. The resulting assembly comprised 171.9 Mb of contigs in 196.5 Mb of scaffold.
The genome assembly was soft-masked to highlight consensus repeat family sequences predicted de novo by RepeatModeler v.1.0.7. A total of 37,384 RNA-seq transcript assemblies were constructed from 118.7 million 2 × 150 paired-end In summary, population genetic evidence for a selective sweep corresponds with our ecological findings: the inversion affects multiple ecologically important traits, including flowering differences that are expected to increase reproductive isolation between subspecies. The inversion occurs in a hybrid zone, and this genomic region contains multiple ecologically important QTLs, as predicted by the Kirkpatrick-Barton model. However, beyond the potential advantage conferred by recombination suppression 9 , closely related West-inv and West-std haplotypes show divergent phenotypic effects, which might contribute to selection in the hybrid zone. Our analysis of a young inversion is compatible with evolutionary predictions that linked, locally adapted QTLs may be captured by new inversions and contribute to local adaptation or incipient speciation 9 .

Methods
Study area. The inversion was originally detected 21 at the Lost Trail Meadow site (~2,500 m elevation) about 4 km northwest from Lost Trail Pass on the Montana-Idaho border. Pollen and geological analyses 22,23 show that climate and vegetation in this area was strongly influenced by Pleistocene climatic changes, although the site apparently was unglaciated during the Last Glacial Maximum about 21 ka. Inferred plant communities 22,23 suggest that conditions suitable for B. stricta may have existed at or near our study site during the Last Glacial Maximum. We searched extensively in nearby B. stricta habitats; the Dryad Data Archive contains information on 122 genotypes from the inversion zone ( Fig. 3) and 83 comparison genotypes from the across-the-species range.

NATURE ECOLOGY & EVOLUTION
Illumina RNA-seq reads using PERTRAN (Shengqiang Shu, Joint Genome Institute in-house pipeline, http://jgi.doe.gov/wp-content/uploads/2013/11/ CSHL-PERTRAN-Shengqiang-Shu-FINAL.pdf). Loci were determined by BLAT transcript assembly alignments, BLASTX alignments of proteins from arabi (Arabidopsis thaliana), rice, soybean, grape, maize or Chlamydomonas reinhardtii genomes, and/or BLASTX alignments of UniProtKB/Swiss-Prot to the B. stricta genome. Gene models were predicated by homology-based predictors, mainly FGENESH+ and GenomeScan. The best scored predictions for each locus were selected using multiple positive factors, including expressed sequence tag and protein support, and one negative factor: overlap with repeats. The selected gene predictions were improved using PASA r2011_05_20. Improvement includes adding untranslated regions, splicing correction, and adding alternative transcripts. PASA-improved gene model proteins were subject to analysis of their protein homology to the abovementioned proteomes, in order to obtain a Cscore and protein coverage. Cscore is a protein BLASTP score ratio to mutual best hit BLASTP score, and protein coverage is the highest percentage of protein aligned to the best of the homologues. PASA-improved transcripts were selected on the basis of Cscore, protein coverage, expressed-sequence-tag coverage, and the percentage of their coding regions that overlapped with repeats. A transcript was selected if both its Cscore and protein coverage were ≥ 0.5, or if it exhibited expressed-sequence-tag coverage, but < 20% of its coding region overlapped with repeats. For gene models where > 20% of the coding region overlapped with repeats, to be selected, its Cscore and homology coverage had to be at least 0.9 and 70%, respectively. The selected gene models were subject to Pfam (protein family database) analysis and gene models whose protein included > 30% transposable-element-associated Pfam domains were removed.
There are 1,591 annotated genes in the inversion region (below), of which 408 are found in the QTL regions within the inversion.
Whole genome profiling. A Hind III BAC library of B. stricta (Bs_LBa) was built from fresh leaf tissues from the LTM reference genotype 18 . The library contains 18,432 clones, arrayed into 48-384 well plates, with an average insert size of 150 Kb and an estimated genome depth of 11× . The library or clones of interest from it are publicly available at the Arizona Genomics Institute (AGI) web page (http://www.genome.arizona.edu/orders/).
A WGP physical map of B. stricta was built following a previously published protocol 32 . Briefly, BAC library clone plates were arrayed into two-and threedimensional pools, and BAC DNA was extracted from the pooled plates. BAC DNA was digested with EcoRI/MseI and ligated to their respective adapters (P5-EcoRI-tag barcode and P7-MseI), followed by PCR amplification and sequencing with an Illumina HiSeq 2500.
Resulting sequences were deconvoluted with KeyGene proprietary scripts (licensed to AGI), to generate the tag file (Dryad Data Archive), which were assembled with FPC 33 v.9.4. Four test assemblies were performed with a fixed tolerance value of 0 and cutoff values of 1 × 10 −25 , 1 × 10 −20 , 1 × 10 −15 and 1 × 10 −10 , to choose the one producing the best results; based on the number of contigs and number of clones included in those contigs.
Questionable clones were eliminated from the best project (1 × 10 −15 ) and the contigs were merged at a cutoff value of 1 × 10 −9 . Singletons (clones that do not assemble with these parameters), were added at reduced stringency (1 × 10 −9 ). This new project was subject to manual editing, after linking the tag, BAC-end and genome-scaffold (v 1.0) sequences to the physical map 34

NATURE ECOLOGY & EVOLUTION
package 35 , the final edited map was aligned to pseudomolecules of B. stricta and Capsella rubella 36 . Synteny with A. lyrata was visualized using MUMmer 37 . The integration of the WGP physical map with the genetic map and sequence scaffolds (via the WGP tags and BAC-end sequencing) showed high concordance among these three genomic resources. WGP enabled regions with low recombination in the linkage map to be ordered and oriented on the basis of their correspondence with the physical map contigs. In addition, we were able to integrate some scaffolds that were not linked in the genetic map, providing more robust pseudomolecules.
Variant discovery. We aligned each genotype read to the B. stricta LTM hardmasked reference (B. strict_278_hardmasked) with BWA 38 . We used GATK 39 for base quality recalibration, indel realignment, simultaneous SNP and indel discovery via HaplotypeCaller and joint genotyping using the default hard filtering parameters as prescribed by GATK best practice 39 .
Ordering and orienting scaffolds into linkage groups. We created a genetic map based on 159 sixth-generation RILs bred from two parent individuals: LTM (the reference individual) and SAD12. We (re-)sequenced LTM to a genomic depth of ~400× using paired-end Illumina sequencing. These reads were aligned to the reference sequence using bwa mem 40 and samtools mpileup 41 to visualize the alignment. We conservatively verified 105,074,494 positions as homozygous by requiring unambiguous alignment of a single base type with a depth ranging from 160-600. Next, we performed a similar analysis on a set of paired-end Illumina reads from SAD12 (depth ~170× ) and compiled a catalogue of 442,637 discriminatory positions throughout the genome. These markers were homozygous for different bases in LTM and SAD12, allowing identification of ancestry within the RILs. These results imply a nucleotide divergence of 0.000425 between these two accessions.
We initially aligned sequences from 199 barcoded RILs, sequenced to modest depth (a few fold), to the reference sequence and genotyped these at all sites

NATURE ECOLOGY & EVOLUTION
possible within the discriminatory catalogue. During this process, we eliminated a number of contaminated RILs, RILs with too little sequence for reliable genotyping, and RILs that appeared to be heterozygous over most of the genome, suggesting that the source DNA had not originated from a single RIL. The scaffolds in the assembly with sizes ≥ 20 kb were binned into 20-kb blocks and each block with at least three genotypeable sites was genotyped as either LTM, SAD12 or heterozygous on the basis of the consensus of genotypeable sites. A catalogue of breakpoints was constructed, as a list of scaffolds and bins where one or more RILs changed genotype. From this data, we identified 9 redundant RILs, and after removal of these, we were left with 159 RILs to use in the analyses. In addition, we found two misjoins in the original assembly, which identified themselves as being apparently unlinked neighbouring bins. These were Scaffold25219, between bin 134 and 137 (removed bins 135,136 and renamed bins 137-139 to Scaffold100000 bins 0-2) and Scaffold19424 between bins 344 and 345 (renamed Scaffold19424 bins 345-391 to Scaffold200000 bins 0-46). In the final data set, 2,664 crossovers were located (16.7 per RIL), and 2.4% of all blocks were genotyped as heterozygous (an average of 3.12% is expected for F6 RILs, with a large variance). Next, we calculated a matrix of recombination fractions, r, between each pair of bins, defined as the fraction of RILs for which the bins differ in genotype. Bins in each pair with r < 0.15 were clustered by a single-linkage approach.
(Briefly, all pairs of bins were compared, and pairs with r < 0.15 were assigned to the same cluster. Two clusters were merged if any pair of members had r < 0.15.) This algorithm resulted in seven distinct linkage groups. The bins in each linkage group were used as markers for input into MSTmap 42 (v.1). Finally, individual bins within assembled scaffolds were rearranged into the order in which they occur in the scaffolds, and in the orientation, if known, consistent with the genetic map. This step is necessary as the sizes of the 20-kb bins are often smaller than the map resolution, so that the bins appear in random order on the map. During this step, the cumulative physical length of the map (in bases) was also inserted into the map. In addition, occasional unmapped bins, the position of which could be inferred from the context of adjacent mapped bins, were inserted at the appropriate locations (Dryad Data Archive). Finally, from the matrix of markers and genotypes, we inferred the locations of recombination events near and within the inversion (Supplementary Fig. 1).

Identification of centromere positions in the B. stricta reference genome.
To identify the centromere position for each B. stricta linkage group, we downloaded sequences of 16 A. thaliana centromeric BAC clones and looked for their homologues in the B. stricta reference genome using BLAST. These 16 clones were mapped on pericentromeric regions of B. stricta by chromosome painting 43 , thus the positions of their homologues on B. stricta linkage groups indicate the margins of centromeres. The lengths of predicted centromere regions were 2-3 Mb for linkage groups 1, 2, 5, 6 and 7, and 6-7 Mb for linkage groups 3 and 4. As expected, all of these predicted centromere regions showed low recombination rates, except for the parts of these regions on linkage groups 3 and 4 (Supplementary Fig. 4 and Supplementary Table 6).

Identification of a paracentric inversion by comparative chromosome painting.
Painting probes. Previous F2 genetic mapping in the LTM × SAD12 F2 population found an extensive block of recombination suppression on LG1 (Bs1). On the basis of the published karyotype structure of B. stricta SAD12 43, we designed BAC-based painting probes to verify the position and orientation of ancestral genomic blocks on the seven chromosomes of B. stricta LTM by comparative chromosome painting. The A. thaliana BAC contigs used as painting probes are listed in ref. 44 . Chromosome-specific BAC contigs were arranged and differentially labelled following the organization of genomic blocks in the reconstructed karyotype of B. stricta 43 .
Chromosome preparations. Chromosome spreads were prepared according to a previously published protocol 45 with minor modifications. Entire inflorescences were fixed overnight in ethanol:acetic acid (3:1) and then stored in 70% ethanol at 20 °C until further use. Before chromosome spreading, closed flower buds with white (young) anthers were rinsed with distilled water and then with citrate buffer (4 mM citric acid and 6 mM trisodium citrate, pH 4.8) and digested using a pectolytic enzyme mixture (0.3% cellulase, cytohelicase, and pectolyase; all Sigma) in citrate buffer at 37 °C for 3-6 h, and then kept in citrate buffer until used. An individual flower bud was put on a microscope slide and dissected using needles in a drop of citrate buffer to form a fine suspension. Then 15 to 30 μ l of 60% acetic acid was pipetted into the cell suspension, which was then spread over the slide. The slide was placed on a hot plate at 50 °C for ~0.5-2.0 min. The spread chromosomes and nuclei were then fixed by pipetting 100 μ l of ethanol:acetic acid (3:1) fixative around the suspension drop. The slide was tilted to remove the fixative, dried using a hairdryer and postfixed with 4% formaldehyde in distilled water for 10 min and left to air dry.

Fluorescence in situ hybridization and comparative chromosome painting.
Before pipetting the probe onto selected slides, the slides were treated with pepsin (0.1 mg ml −1 ; Sigma) in 0.01 M HCl for 3-6 min, postfixed in 4% formaldehyde in 2× saline sodium citrate buffer (300 mM sodium chloride and 30 mM trisodium citrate, pH 7.0) for 10 min, dehydrated in an ethanol series (70, 80 and 96%; 3 min each), and air dried. The labelled BAC clones (and rDNA probes) were pipetted together and the DNA was precipitated by adding a 0.1 volume of 3 M sodium acetate (pH 5.2) and a 2.5 volume of ice-cold 96% ethanol, kept at − 20 °C for at least 30 min, and centrifuged at 13,000 g at 4 °C for 30 min. The pellet was dried using a desiccator and resuspended in 20 μ l of hybridization buffer (50% formamide and 10% dextran sulfate in 2 × SSC) per slide. The probe and chromosomes were denatured together on a hot plate at 80 °C for 2 min and hybridization was carried out by placing the slides in a moist chamber at 37 °C for 40-63 h. Post-hybridization washing was performed in 20% formamide in 2× SSC at 42 °C. Signal detection and amplification were achieved as follows: biotin-dUTP was detected by avidin-Texas Red (1:1,000; Vector Laboratories) and amplified using goat anti-avidin-biotin (1:200, Vector Laboratories) and avidin-Texas Red; digoxigenin-dUTP was detected by mouse anti-digoxigenin (1:250; Jackson Immuno-Research Laboratories) and goat anti-mouse Alexa Fluor 488 (1:200; Molecular Probes); Cy3-dUTP was observed directly. DNA was counterstained with 4′ ,6-diamidino-2-phenylindole (2 μ g ml −1 ) in Vectashield (Vector Laboratories). The hybridization signals were analysed using an Olympus BX-61 epifluorescence microscope equipped with fluorochrome-specific excitation and emission filters (AHF Analysentechnik), and a Zeiss Axio-Cam CCD camera. The monochromatic images were pseudo-coloured and processed using Photoshop CS5 (Adobe).
Mapping SNPs to pseudomolecules. In the genetic map, 208 scaffolds were ordered and oriented into 7 LGs. Misjoins were found in the original assembly of Scaffold25219 and Scaffold19424, which were split into four scaffolds, two with the original names (Scaffold25219 and Scaffold19424) and two with new names Scafffold100000 and Scaffold200000. In Scafffold25219, the first 135 bins (bin: 0-134; position: 1-2,700,000) were found to be unlinked with the last three bins (bin: 137-139; position: 2,740,001-2,797,113 bp), while the position of bins 135-136 (position: 2,700,0001-2,740,000) between them could not be determined. Thus, the first 135 bins were retained in Scaffold25219 (a shorter one; position: 1-2,720,000), the last three bins were moved to a new Scaffold100000 (position 2,740,001-2,797,113 of Scaffold25219 was changed to position 1-57,113 of Scaffold100000), and bins 135-136 were excluded and not used in map. In Scaffold19424, the first 345 bins (bin: 0-344; position: 1-6,900,000) were found to be unlinked to the remaining 47 bins (bin: 345-391; position: 6,900,001-7,828,947), thus bins 0-134 were kept in Scaffold19424 (position: 1-6,900,000) and the rest were moved into a new scaffold200000 (position 6,900,001-7,828,947 of Scaffold19424 was changed to position 1-928,947 on Scaffold200000).
Accordingly, SNP coordinates were extracted from the vcf file (scaffold name and position on scaffold) and examined as to whether they fell in a bin of the linkage map on the basis of its scaffold position. If so, its position on the pseudomolecules (g) was calculated as g = A − |g − c|, where A is the end position of

NATURE ECOLOGY & EVOLUTION
Gene prediction. We obtained the 10 kb region on each side of each breakpoint in the LTM reference genome (inversion) and created pseudomolecules corresponding to the corresponding regions around each breakpoint in the standard haplotype. All four sequences were submitted to the Augustus gene prediction algorithm (http://bioinf.uni-greifswald.de/augustus/), and we found no open reading frames spanning the breakpoints in either the inversion or standard haplotypes. Thus, the chromosomal inversion does not disrupt existing open reading frames or create new ones (Supplementary Fig. 7). Genes predicted by Augustus largely correspond to the transcriptome-assisted annotation in the B. Stricta genome v1.2.
Gene expression. To investigate whether the inversion event changed the expression of flanking genes, we compared two nearby populations using F2 progeny from the LTM × SDM cross (representing inversion and standard haplotypes, respectively). On the basis of the inversion genotypes, we chose 11 F2 plants that were homozygous for the inversion and 12 that were homozygous for standard haplotypes, providing sample sizes well above those recommended by power analysis 51 . All 23 plants were from cross CL9.1 where LTM was the mother. These F2 plants were self-fertilized for one generation, and two-month-old rosette leaves from Bsi1 homozygous F3 plants were used to compare the expression of seven genes flanking the two breakpoints and syntenic to the region in A. thaliana. We used a Sigma Spectrum Plant Total RNA kit to extract the RNA, and a Thermo Scientific DyNAmo cDNA Synthesis kit to synthesize the cDNA. We used a Thermo Scientific DyNAmo SYBR Green qPCR kits for quantitative PCR (primers in Dryad Data Archive). As in previous gene expression studies of B. stricta 18 , we used ACTIN2 (ACT2) as a reference and calculated the expression as Δ C t = C t ACT2 − C t gene , where C t is the cycle threshold, C t ACT2 is the C t value for ACT2 and C t gene is the C t value for each gene. Fold gene expression relative to ACT2 was calculated as Δ C 2 t . Since Δ C 2 t has a skewed distribution, we analysed results from both Δ C 2 t and Δ C t . For each gene separately, ANOVA was performed with inversion status as a fixed effect (Supplementary Table 5). In addition, we examined the statistical power using JMP Pro v.12.0.1 (SAS) and via simulations in R.

Dissecting the inversion region in a collinear cross. The East × West
'Parker × Ruby' cross used for QTL mapping was developed from two parents in the East-West contact zone: one in Parker Meadow (RP105, Parker, East subspecies; 44° 37′ N, 114° 31′ W) and one at Ruby Creek (RP109, Ruby, West subspecies; 45° 33′ N, 113° 46′ W). On the basis of nucleotide similarity ( Supplementary  Fig. 8a, left), the West-std parent (RP109) used in this cross was very similar to the most recent common ancestor that gave rise to the inversion and its closest West-std relatives (IN086 and IN087). The F1 hybrid was self-fertilized to produce F2 plants, and subsequent generations were propagated by self-fertilization and single-seed descent to create 153 independent genetic lines (families). In each line, multiple F4 progeny from the same F3 plant were used in a randomized complete block design totalling 1,714 individuals. The phenotypic least-squares means were calculated in JMP 8 to represent the genotypic value for their F3 parent.
Genotyping by sequencing. We used genotyping by sequencing (GBS) 52 to identify SNPs in this cross. In each family, DNA from each accession was extracted from 0.1 g of young leaf tissue following a dark treatment of about three days. Tissue was stored at − 80 °C, flash-frozen in liquid nitrogen and homogenized prior to genomic DNA extraction using Qiagen DNeasy Plant Mini kits. DNA concentration was measured using a Qubit fluorometer (Turner Biosystems, Invitrogen) from at least ten pooled F4 individuals to represent the genotype of their F3 parent. This protocol uses an adaptor design that is compatible with TruSeq adaptors and indexes while allowing paired-end sequencing. The combination of 48 unique barcodes with four different TruSeq indexes allowed multiplexing of 192 samples (153 F3:F4 families, 19 inbred individuals of the Parker parent, and 20 inbred individuals of the Ruby parent). For sequencing, we used an Illumina HiSeq-2000 or HiSeq-2500 at the Duke Genome Sequencing and Analysis Core Resource. From this population we obtained ~249 million reads with unambiguous barcodes (Sequence Read Archive accession, SRP075905). Read pairs were assigned to genotypes and parents by custom Perl code, and low-quality bases in the end of reads were trimmed. SNPs were called using GATK Best Practices (above). In addition, these protocols were used to assay GBS SNPs (Sequence Read Archive, accession SRP075997) in 122 genotypes across the inversion area.
Because the LTM × SAD12 cross does not recombine in the inversion on chromosome 1, we used the Parker × Ruby cross to infer the linkage map for bins within the inversion. The mean sequencing depth of GBS SNPs in the Parker and Ruby parents was ~20× , with a standard deviation of ~40× . We focused on 52,827 bi-allelic SNPs where the two parents had sequencing depths ≥ 6× and ≤ 100× (the mean depth plus two standard deviations), and were homozygous for different alleles. We binned the genome into 100 kb windows, and for each individual we counted the number of reads from the two parents, and calculated the proportion of Parker-derived reads. For an individual, if the cumulative depth of all SNPs in a window was less than 20× , it was classified as missing. We removed windows with ≥ 30% missing data. We used a hard-genotype cutoff based on the proportion of Parker reads: a Parker proportion ≤ 0.2 was defined as homozygous Ruby, a Parker proportion ≥ 0.7 was defined as homozygous Parker, and the remaining windows the bin that the SNP fell in, g is the SNP's position on the scaffold, and c is the end coordinate of the bin on the scaffold. For SNPs that did not fall in any bins but were located on regions 2,740,001-2,797,113 of Scaffold25219 or 6,900,001-7,828,947 of Scaffold19424, their positions on Scaffold100000 and Scaffold200000 were calculated as g = g 1 − 2,740,000 or g = g 1 − 6,900,000, respectively (g 1 is the position on the original scaffold and g is the position on the new scaffold). After that, their positions on pseudomolecules (g′ ) were calculated as described above.
Delimiting and genotyping the inversion. After chromosome painting narrowed down possible locations of the inversion breakpoints in A. thaliana, we blasted these Arabidopsis genes onto the B. stricta LTM genome and identified putative scaffolds containing the two breakpoints. To identify the exact breakpoints, we mapped paired-end Illumina reads from SAD12 onto the LTM scaffolds with BWA 38 , followed by de novo assembly with Velvet 46 for read-pairs having one or both ends near the putative breakpoint region in the LTM genome.
The primers Inv-A, Inv-B and Inv-E were designed to amplify the inversion region using Primer 3 (ref. 47 ). Primer Inv-D was used for Sanger sequencing the breakpoints in the two karyotypes for the confirmation of breakpoints. Locations of the breakpoints and related primers are shown in the Dryad Data Archive.
Thermocycling consisted of heating to 94 °C for 3 min, then 33 cycles of 94 °C for 20 s, 55 °C for 1 min, 72 °C for 30 s, followed by 72 °C for 6 min. Amplified fragments were separated by electrophoresis on a 1% agarose gel stained with SYBR Safe (Thermo Fisher). Fragments were ~300 bp in accessions containing the inversion, and ~600 bp in accessions without the inversion.
Differences between subspecies. From our earlier study of ecological divergence between East and West subspecies 16 , we used least-squares means from 24 genotypes (Supplementary Table 1) to quantify differences in flowering time. In the greenhouse, East genotypes flower an average of 5.0 days before West genotypes.
Here we report on three East × West crosses, which all gave healthy advanced generation progeny in F2, F4 or F6 generations. More generally, intrinsic reproductive isolation might be present in some instances, especially if complicated by occasional polyploidy or apomixis 48 .

Effect of the inversion on NIL flowering time.
In the LTM × SAD12 (MT West-inv × CO East-std ) cross 17 , we generated two F7 NIL families from the F6 population to investigate the inversion effect. Two F6 heterogeneous inbred families (HIFs: 116A and 120A) 49 heterozygous for the inversion and mostly homozygous elsewhere in the genome were self-fertilized, and 144 F7 plants from each parent were grown in a completely randomized design in the greenhouse. After three months, plants were vernalized at 4 °C for six weeks, and flowering time was recorded. We also genotyped the inversion status in all individuals using the inversion primers Inv-A, Inv-B and Inv-C. We used fixed effects ANOVA to test the equation, Flowering Time = Family + Inversion + Family × Inversion, and verified compatibility with parametric statistical assumptions.
For families from this MT West-inv × CO East-std cross, we also analysed lifetime fecundity in the field, using Dryad-archived data (http://dx.doi.org/10.5061/dryad. rp3pc) from the Lost Trail field site in the inversion zone 20 . We analysed lifetime fecundity in an East × West recombinant inbred F6 population grown in nature for three years.

Effects of the inversion in a sympatric cross.
The inversion might be favoured if this rearrangement alters genes around the break points and gives it higher fitness over its ancestral standard haplotypes from the same population. To test this, we generated crosses between the reference accession, LTM (inversion), and the accession, SDM, which has the standard haplotype. These genotypes are of the same subspecies (West) and were collected from similar habitats ~1.7 km apart. Two crosses were generated, one with LTM as mother (CL9.1) and the other with SDM as mother (CL10.1). The F1 hybrids were self-fertilized, and 1,000 F2 plants were grown in the greenhouse. The inversion was genotyped in all individuals, and we measured 11 life history traits (Survival after vernalization (binary), flowering (binary), width at 4 weeks, width at 10 weeks, flowering time (days after vernalization), flowering width, flowering height, flowering rosette number, flowering leaf number, fruit number, lifetime fitness and inversion genotype) on 994 individuals.
To test the phenotypic effects of the inversion, we permuted the inversion genotypes in the crosses 1,000 times, and performed a MANOVA on the vector of traits using R 53 (v.3.1.2). Because R only reports a type-I sum of squares for MANOVA, we obtained a type III sum of squares by performing several MANOVAs, and reported the effect of each factor when added to the model last. We modelled Multivariate Traits = Cross + Inversion + Cross × Inversion, with cross type (either CL9.1 or CL10.1 as maternal parent), three inversion genotypes (inversion, standard, or heterozygote), and their interaction as fixed-effect predictor variables. P-values were obtained by comparing the true Wilks' Lambda statistic to those from 1,000 permutations (Supplementary Table 4). Because MANOVA showed that multivariate means differed between Bsi1 genotypes, the univariate traits were analysed by ANOVA with a P = 0.05 significance threshold 50 .

9
© 2017 Macmillan Publishers Limited, part of Springer Nature. All rights reserved.

NATURE ECOLOGY & EVOLUTION
were classified as heterozygous. These cutoffs gave the correct hard-genotype proportions for an F3 population (about 37.5% of either homozygote and 25% heterozygous genotypes). After determining the 'hard genotype calls' , we calculated the proportion of the three genotypes in each window. To remove windows with excessive proportions of heterozygotes or either homozygote, we excluded windows where the proportion of the three genotypes was beyond the upper or lower 5% of their respective whole-genome distribution. Given this was a cross in the F3 generation, which had only experienced two generations of recombination, we regarded it as extremely improbable that two recombination events would occur within 1 Mb of each other. Therefore, in each individual, if two recombination events were identified within 1 Mb, we changed the genotype calls of markers between two breakpoints to missing. The linkage map was constructed with MSTmap. With only three windows (markers) dropped out, MSTmap constructed seven linkage groups, with scaffold-chromosome assignment that was consistent with the LTM × SAD12 cross. The final linkage map was refined in multiple steps: we first removed markers that were > 5 cM away from both flanking markers and re-constructed the map. In the new map, we imputed a missing marker genotype if the upstream and downstream markers with data had the same genotype call and were not both > 10 cM away from the missing position. Another linkage map was constructed, and the imputation process was repeated since the new map had slightly changed the order of some markers, making some missing genotypes now imputable. The final linkage map was then constructed with 1,010 markers. Thus, we used the Parker × Ruby cross to infer the linkage map for bins within the inversion and we inferred the position of polymorphisms within the inversion from their positions in each ordered and oriented contig.
This collinear East × West 'Parker × Ruby' greenhouse experiment consisted of 12 blocks, each with one F4 plant from each of the 153 lines and multiple Parker and Ruby individuals (N = 1,714). Seeds were stratified in 4 °C for four weeks and planted in 'Cone-tainers' (Ray Leach SC10; Stuewe and Sons), with soil composition and greenhouse conditions as previously described 16 . When rosettes were 11 weeks old, all leaves from three blocks of plants were harvested for rosette-and leaf-morphology measurements as described 16 . At 12 weeks of age, the remaining nine blocks were vernalized at 4 °C for 6 weeks, then returned to the same greenhouse conditions for phenology measurements. All traits were measured as previous described 16 , except that (1) no physiological traits were measured and (2) leaf width/length ratio was used instead of leaf shape morphometrics because the leaf-shape landscape points were highly correlated 16 (Supplementary Table 7).
For QTL mapping, all individual-level measurements were transformed to family-level least-squares means in JMP (N = 153 families). Due to the skewed distribution of most phenotypes, all traits (except binomial traits or plant stages) were log-transformed at the individual level. For greenhouse measurements, Block and Genotype were considered as random effects. We used the R libraries qvalue and MASS for MANOVA, false discovery rate and discriminant function analysis. For each marker across the inversion interval and adjacent contigs, we fitted a MANOVA model: Multivariate Phenology Traits = Block + Single Marker Genotype. At each marker we computed the Wilks' Lambda, permuted the vector of phenotypes with respect to marker genotypes 1,000 times, and finally corrected for false discovery rates across the inversion region 54 .
DFA was performed in R for a marker at the peak of each multivariate QTL (bins: Scaffold-25219_50000, Scaffold13671_1150000, and Scaffold26675_450000), defining the direction of trait variation along the axis of greatest differentiation among the multivariate means of these marker genotypes. Subsequently, at these three markers we used the eigenvector for the first DFA axis to compute the phenotypic projection of each genotype on the trait axis identified by DFA. These new trait values were used for univariate QTL mapping across the inversion region, with correction for false discovery rates across the inversion region.
Finally, we annotated SNPs using snpEff 55 v.4.2 and found 17 genes that are orthologues of the flowering-time genes of A. thaliana 56  Population genetics. SNP filtering. SNPs called by GATK were filtered in each of three groups: sets consisting of 35 inversion genotypes (G INV ) and 87 standard genotypes (G STD ) from the inversion zone, and a 'Reference Population' set of 83 genotypes from across the species range, outside the inversion zone (G RP ). Next, data from these single groups were pooled together into two datasets, the complete dataset (G INV + G STD + G RP ) and the inversion zone genotypes (G INV + G STD ). In each dataset, genotypes supported by less than two reads were assigned as missing, and SNPs were discarded if they met the following criteria in any of the groups: (1) detected in fewer than 50% of individuals; (2) mean depth > 20; (3) more than one variant allele was observed; (4) sites where the proportion of heterozygous genotypes was > 15% (B. stricta is predominantly inbred 17 , hence high heterozygosity may indicate paralogous loci); (5) reference or variant alleles were indels. After filtering, 75,737 and 43,722 SNPs were retained for the G INV + G STD + G RP and G INV + G STD sets, respectively. The following analyses used these SNPs.
Ancestral allelic state. Several population genetic statistics require information on the ancestral state of segregating variants. Therefore, we sequenced a Boechera holboellii (B. holboellii; reference genotype 'Panther'; location: 45° 18.198′ N, 114° 22.599′ W) individual with a mean depth of ~400× (Joint Genome Institute project ID: 1051698). Short reads were mapped to the B. stricta reference genome using BWA, genotypes were called using GATK with default settings, and sites were filtered out if they had a depth of < 20 or > 800, or if they were heterozygous. A Python (v.2.7.9) script assigned ancestral alleles as the shared allele between B. holboellii and B. stricta. In total, we obtained information on the ancestral state for 63,182 (83.4%) of the 75,737 G INV + G STD + G RP SNPs and for 35,100 (80.3%) of the 43,722 G INV + G STD SNPs.
Windows. We partitioned the genome into 300 kb windows. The sequence coverage of each window was calculated by counting the number of available sites (both variant and non-variant sites obtained from GATK using the '-allSites' argument) per window passing the following quality filters: detected in at least 50% of individuals, mean depth ≤ 20, and proportion of heterozygous genotypes ≤ 15% in each of the three groups (G INV , G STD , and G RP ). Windows were distributed on a per-scaffold basis, beginning at position 1 of a scaffold, and were then oriented along pseudomolecules according to the linkage map. After excluding windows with sequence coverage < 5 kb, 2,801 windows were retained for downstream analyses.
Phylogenetic relationship and population structure analyses. Python scripts were used to generate alignments from genotypes in the vcf file, with missing and heterozygous loci coded as 'N' . With genome-wide SNPs, we verified the East-West population structure using principal component analysis in EIGENSOFT 57 v.6.0.
To investigate the origin of the Bsi1-inv haplotype we used data from 901 SNPs in the inversion region of the genome to examine relationships among 35 Bsi1-inv and 83 Bsi1-std haplotypes. (Four admixed Bsi1-std genotypes were excluded.) First, we constructed neighbour-joining trees using MEGA 58 v.6.06 with 1,000 bootstrap steps. Second, maximum-likelihood (ML) trees were constructed using RAxML 59 v.8.0.0. One thousand ML trees were generated to find the best-scoring ML tree, and topological robustness was investigated using 1,000 non-parametric bootstrap replicates. Neighbour-joining and ML trees were displayed using FigTree 60 .
Site frequency spectrum and summary statistics. To control for error rates and variable coverage of short-read sequencing, we used a probabilistic method implemented in ANGSD 61 v.0.911 to estimate the site frequency spectrum (SFS) and related population genetic statistics. A number of filtering steps were performed: (1) removed reads with a minimal mapping quality of 30 and bases with a minimal quality score of 30 (-minMapQ and -minQ); (2) removed sites with information from less than 50% of individuals (-minInd); (3) removed sites with a P-value higher than 1 × 10 −6 (-snp_pval); (4) assigned genotypes as missing if the depth was less than two for an individual; (5) only used genotypes with a posterior probability higher than 0.95; (6) removed sites that did not pass previous filtering criteria (above). Because B. stricta is predominantly inbreeding 17 , we estimated inbreeding coefficients for each individual in ngsTools 62 v.1.0 and incorporated them into the calculation of SFS in ANGSD. Using genotype likelihoods based on the GATK genotyping model 63 , we estimated folded and unfolded SFS and derived a set of population genetic summary statistics in 300 kb windows. For the G INV and G STD groups, we estimated Tajima's D 64 on the basis of both folded and unfolded SFS, and calculated Fay and Wu's H 29 , and Fu and Li's D and F 65 using the unfolded SFS. The inferred ancestral allelic state (based on short reads from the B. holboellii genotype 'Panther') aligned to the B. stricta reference genome was used to estimate the unfolded SFS. We compared population genetic statistics in the inversion region, on other chromosomes (LG2 to LG7), and in Block D 21 , a chromosome region with unusually high polymorphism in related species, perhaps due to clusters of genes for nucleotide-binding site-leucine-rich repeat proteins 66 .
F ST , D xy , D A and π. Genetic differentiation (F ST ) 67 between G INV and G STD groups was estimated using VCFtools 68 v.0.1.12 and we used Python scripts to estimate nucleotide diversity (π ) 69 within each group, pairwise nucleotide divergence (D xy ) 70 and net pairwise nucleotide divergence (D A ) 70 between groups. All parameters were calculated on a per-site basis, and then averaged to obtain window-based estimates (300 kb windows). The window-based F ST was calculated by averaging per site F ST across all variable loci, and the window-based π , D xy and D A were estimated by averaging all sites (both variable and monomorphic) passing the initial quality filters for each window (above).
Linkage disequilibrium. The level of linkage disequilibrium (LD) between the inversion and all SNPs on chromosome 1 (including those in the inverted region) was estimated as the mean-squared correlation (r 2 ) for 118 samples from inversion zone accessions (35 Bsi1-inv and 83 Bsi1-std) using plink 71 v.1.90 . (Four admixed Bsi1-std samples were excluded; see above.) In addition, to compare LD between the inversion and flanking SNPs, we identified ten randomly chosen comparator SNPs from LG2 to LG7, having similar derived allele frequencies (25-35%; close to the 29.7% frequency of Bsi1-inv in the inversion zone), and then calculated LD

NATURE ECOLOGY & EVOLUTION
between them and flanking SNPs along the chromosomes. For LG1 (Fig. 5), strong LD (r 2 > 0.4) between the inversion and SNPs within and outside the inverted region extended to the end of the chromosome (more than 10 Mb). In contrast, LD with the ten comparator SNPs from LG2 to LG7 declined quickly to background levels ( Supplementary Fig. 8b).
Age of the inversion. To estimate the age of the inversion, we phased the 35 Bsi1-inv genotypes using Beagle 72 v.4.1 with default settings, and estimated pairwise genetic distances among haplotypes on the inverted genome region. After that, the divergence time between haplotypes was calculated as T = d/2μ, where d is pairwise genetic distance and μ is the mutation rate. We used a mutation rate of 7 × 10 −9 per site per generation in A. thaliana 24 ) with a generation time of two years in B. stricta. To avoid biases due to missing data, we only considered 1,675 haplotype pairs with less than 50% missing sites in the inversion region. We approximated the lower bound for the age of the inversion on the basis of the mean pairwise distances 73 , and the upper bound from the maximum of all pairwise distances among inversion genotypes.
Divergence time between QTL alleles in the collinear cross. Using the same mutation rate and generation time, we estimated D A = D xy = 2μT, and solved for T between these two alleles in 200-kb non-overlapping windows 71 . This window size was chosen to ensure sufficient polymorphic sites within each interval.
Frequency distribution of derived alleles. We tested whether Bsi1-inv had unusually high frequency among derived SNP alleles in the inversion zone population. For this analysis, we excluded one admixed, genetically divergent inversion individual (IN019; Supplementary Fig. 9c, right), and analysed Bsi1-inv and closely related Bsi1-std genotypes from the inversion zone ( Supplementary Fig. 9c, left; first principal component < 0.0). These genotypes comprise 54 samples (34 Bsi1-inv and 20 Bsi1-std), with 2,416 SNPs that, like the inversion, segregate among these genotypes and are monomorphic in the rest of our collection.
Code availability. All code is available in the Dryad Data Archive at http://dx.doi. org/10.5061/dryad.tt743.
Data availability. The assembly and annotation of the Boechera stricta genome are available at https://phytozome.jgi.doe.gov/pz/portal.html#!info?alias= Org_Bstricta and have been deposited under GenBank accession number MLHT00000000. The BAC end sequences have been deposited under GenBank dbGSS accession numbers KS412618 to KS448200. The short reads of B. stricta (LTM, SAD12 and recombinant inbred lines) have been deposited under GenBank accession numbers SRP078672, SRP048481 and SRP079414, respectively. The short reads from the RNA-seq transcriptome have been deposited under GenBank accession number SRX1971488. The short reads of outgroup B. holboellii have been deposited under GenBank accession number SRP078889, and the short reads of the GBS data of B. stricta have been deposited under GenBank accession numbers SRP075905 and SRP075997. All SNPs used in population genetic analyses have been deposited into dbSNP under accession numbers ss2136554982 to ss2136641742. All other data are available in the Dryad Data Archive at http://dx.doi.org/10.5061/ dryad.tt743.