The genomic distribution of transposable elements is driven by spatially variable purifying selection

Abstract It is widely accepted that the genomic distribution of transposable elements (TEs) mainly reflects the outcome of purifying selection and insertion bias (1). Nevertheless, the relative importance of these two evolutionary forces could not be tested thoroughly. Here, we introduce an experimental system, which allows separating purifying selection from TE insertion bias. We used experimental evolution to study the TE insertion patterns in Drosophila simulans founder populations harboring 1040 insertions of an active P-element. After 10 generations at a large population size, we detected strong selection against P-element insertions. The exception were P-element insertions in genomic regions for which a strong insertion bias has been proposed (2–4). Because recurrent P-element insertions cannot explain this pattern, we conclude that purifying selection, with variable strength along the chromosomes, is the major determinant of the genomic distribution of P-elements. Genomic regions with relaxed purifying selection against P-element insertions exhibit normal levels of purifying selection against base substitutions. This suggests that different types of purifying selection operate on base substitutions and P-element insertions. Our results highlight the power of experimental evolution to understand basic evolutionary processes, which are difficult to infer from patterns of natural variation alone.


INTRODUCTION
Transposable elements (TEs) are DNA sequences that move and amplify within the host genome ( 5 ).Their diverse structures and proliferation mechanisms ( 1 ) are welldocumented across the entire tree of life ( 1 , 6 ).The dynamics of TEs are not only dri v en by their selfish amplification mechanisms, but also by the fitness consequences for the host mediated by TE insertions.This interesting interplay has attracted the attention of evolutionary biologists since the discovery of TEs and continues to do so.
Purifying selection is one of the key evolutionary forces shaping the distribution of TEs in na tural popula tions ( 1 , 6 ) and three different mechanisms of TE-related fitness costs have been discussed ( 7 ).(i) RNAs or proteins r equir ed for transposition are deleterious to the host ( 6 , 8 ), (ii) TE copies at non-homologous sites cause elevated ectopic recombination rates (9)(10)(11), (iii) new TE insertions are deleterious for the host ( 7 , 12-14 ).
Empirical evidence for purifying selection removing TE insertions is either based on 'patterns of absence' in functionally important regions of the host genome (e.g.exons) Figure 1.Schematic ov ervie w of the experimental design.Phase 1 (low efficacy of purifying selection due to small population sizes): acquisition of P-element insertions (r ed r ectangles) during 4.5 years in small populations (isofemale lines with ∼40-50 individuals).Phase 2 (high efficacy of purifying selection due to large population sizes): measuring purifying selection for two classes of P-element insertions: persisting insertions from phase 1 (dark red) and insertions acquired during phase 2 (blue).Isofemale lines with and without P-elements were combined to form three large outbred replica te popula tions in phase 2 ( N = 1250 individuals per replica te).( 15 , 16 ) or low population frequencies in natural populations ( 17 ).This inference makes the implicit assumption of a similar insertion rate of TEs across the entire genome.Since insertion pr efer ences of TEs have been reported ( 1 , 3 , 18 ) the distinction between purifying selection purging deleterious TE insertions and insertion bias is difficult, if not impossible, w hen onl y the genomic distribution of TEs of one time point is used to infer evolutionary TE dynamics.Thus, it has been suggested to study TE dynamics under controlled environmental conditions and varying selection pr essur es ( 19 , 20 ), but such studies are rarely carried out.
In this study, we use Evolve & Resequence (E&R) (21)(22)(23) to distinguish the relati v e importance of insertion bias and purifying selection on the genomic distribution of TE insertions.We take advantage of a natural D. simulans population that is inv aded b y the P-element ( 4 ).The P-element --a 2.9 kb long DNA transposon ( 24 , 25 ) --is one of the best studied TEs in eukaryotes ( 26 ) and the successful invasion of dif ferent Dr osophila species is well-documented ( 4 , 27-29 ).To distinguish insertion bias from selection, our study consisted of two phases, which differed by the opportunity for purifying selection: In phase 1, the P-element proliferated in small populations (isofemale lines) with an acti v e P-element ( 29 , 30 ) for about 4.5 years.In phase 2, we established large outbred populations by mixing isofemale lines with and without P-elements and e volv ed them for 10 generations.This experimental design allowed us to disentangle purifying selection from insertion bias (Figure 1 ).While purifying selection is particularly strong in exons, we also found negati v e selection against P-element insertions in introns and intergenic regions.Very weak selection occurred against P-element insertion sites that have been previously reported to be shared between natural Drosophila melanogaster and D. simulans populations ( 4 ).

Set up and maintenance of experimental populations
We set up three independently maintained replicate populations in June 2015 using 191 isofemale lines collected from a natural D. simulans population in Tallahassee, Florida ( 31 ).It has been estimated that 25-44% of these isofemale lines contain at least one copy of the P-element ( 4 , 29 ).At the time of set up, the isofemale lines had been maintained for 4.5 years at small population sizes ( ∼40-50 individuals).Each replica te popula tion consisted of fiv e 300 ml plastic bottles with 70 ml standard Drosophila medium, and each of these bottles was set up with 191 mated female flies -one female from each of the isofemale lines (corresponding to the F0 flies).The offspring of all fiv e bottles in a gi v en replicate were mixed before the set-up of each new generation.Replicates were maintained with non-overlapping generations and a census size of 1250 flies in a cycling hot temperatur e r egime (12 h light and 28 • C; 12 h dark and 18 • C).

P-element detection
Isofemale lines.To test whether P-element abundance has increased in isofemale lines over time (Supplementary Figure S1, Supplementary File 1), we estimated P-element copy numbers in three different isofemale lines (I116, I174, I211) after 2 and 6 years of maintenance in the laboratory.Isofemale lines were collected from a na tural Dr osophila population in Tallahassee Florida in November 2010 ( 31 ) and have been kept at small population sizes ( ∼40-50 individuals) e v er since.The details of the genomic sequencing and crossing schemes can be found in ( 4) and ( 32 ).For each time point, one individual male from an isofemale line was crossed with a virgin female from the P-element free M252 reference strain ( 4 , 33 ).From each cross, one single female offspring was sequenced.We demultiplexed barcoded files ( -maximumMismatches 3 -maximumN 2 ) and trimmed raw paired-end reads ( -mottQualityThreshold 15 -disab le5pT rim -minReadLength 25 ) with ReadTools (v1.5.2) ( 34 ).We estimated P-element copy numbers with De viaTE ( 35 ).Because De viaTE does not take paired-end read information into account ( 35 ), we only used the first read from each read pair for the analysis.To avoid systematic biases in P-element copy number estimates due to differ ent r ead lengths ( 35 ) (2 y ears: 100 bp, 6 y ears: 150 bp), we clipped the quality-trimmed reads to a maximum read length of 100 bp using the program BBDuk of BBMap (v38.87;sourceforge.net/ projects / bbma p).We ma pped the r esulting trimmed r eads with b wa b wasw (v0.7.17; -M ) ( 36 ) to a r efer ence containing the P-element consensus sequence of Drosophila ( 37 ) and three single-copy genes ( rpl32 , traffic jam, rhino ).DeviaTE estimates P-element copy numbers per haplotype by comparing the read depth of single copy genes with the read depth of the P-element consensus sequence ( 35 , 38 ).
Experimental populations (phase 2).We sequenced the ancestral populations ( n = 3, generation 0, females only) and the e volv ed populations ( n = 3, generation 10, mixed sexes) using the Pool-Seq ( 39 ) approach.Genomic DNA was extracted with a standard salting out protocol ( 40 ) from the pool of all females for a gi v en replicate in the ancestral population (5 × 191 flies = 955 flies) and for half of the pooled individuals from a given replicate in the F10 generation (approx.600 flies).For the ancestral populations, paired-end libraries were generated using the NEB-Next Ultra II DNA Library Prep Kit, starting from 1000 ng DNA sheared with a Covaris S2 (Covaris, Inc.Woburn, MA, USA).Libraries were subjected to double-sided size selection for an insert size of 300 bp using AMPure XP beads (Beckman Coulter, Carlsbad, CA) and amplified with 4 PCR cycles using dual index primers.For generation 10, pair ed-end libraries wer e generated using the NEBNext Ultra II FS DNA Library Prep Kit (New England Biolabs, Ipswich, MA).The protocol was modified to use only 10% of the provided reagents in each step, a target insert size of 300 bp and amplification with dual-index primers and 5 PCR cy cles.For the e xperimental populations in phase 2, all libraries were sequenced on a HiSeq X Ten with a read length of 2 × 150 bp.
For each replicate population, P-element insertions are classified as 'lost', if they are detected in any replicate population after phase 1 (i.e. generation 0), but are not detected in the focal replica te popula tion after phase 2 (i.e. generation 10), whereas 'persistent' P-element insertions are detectable after phase 2. P-element insertions, which occurred during phase 2 in the focal replica te popula tion, and are not detectable after phase 1 in any replicate population, are classified as 'phase-2 (i.e.new insertions).Gi v en that most Pelement insertions occur at a low frequency, not all insertions can be r e-discover ed, which could lead in a false classification of P-element insertions.Ne v ertheless, gi v en that the frequency of P-element insertions is very similar across genomic annotation features (Supplementary Figure S2A), we expect the same re-discovery and false classification rate for all genomic annotation features.P-element frequency distribution.The joint analysis mode of PopoolationTE2 is well-suited to study orthologous Pelement insertion sites in different samples, but it is not wellsuited for allele frequency estimates across multiple evolutionary replicates: low frequency P-element insertions of the same replicate population sequenced in a time-resolved manner ar e mor e likely to be detected and included in the frequency estimates ( 41 ).To estimate the frequency distribution of P-element insertions, we generated a ppileup file (PopoolationTE2; v1.10.03;ppileup ( 41 )) that contains only the three replica te popula tions sequenced at generation 0. PopoolationTE2 estimates population frequencies of TEs as the proportion of physical coverage supporting a TE insertion to the total physical coverage of the population.To obtain roughly the same sensitivity as in the joint analysis with six samples, we subsampled the resulting ppileup to a uniform physical coverage of 60x per sample, which removed 19.12% of the genomic sites from the analysis.We identified and filtered P-element insertions as described above (with the exception of masking regions with an av erage cov erage < 60 × instead of < 50 ×), before estimating frequencies of individual P-element insertion sites with PopoolationTE2 (v1.10.03 pairupSignatures , Supplementary File 3) ( 41 ).

P-element c lassification (separ ate analysis).
To confirm the robustness of our results, we reanalyzed the P-element insertions for each replicate population independently.For this, we generated for each replicate population ( n = 3) an individual ppileup file with PopoolationTE2 (v1.10.03;ppileup ) ( 41 ) that contained the respecti v e e xperimental popula tion sequenced a t genera tion 0 and 10.To obtain roughly the same sensitivity as in the joint analysis with six samples, all resulting ppileup files were subsampled to a uniform physical coverage of 60x per sample ( 41 ), which removed 15.61%, 19.13% and 16.27% of the genomic sites from the replicate-specific analysis.P-element signatur es wer e identified as described above (with the exception of masking regions with an average coverage < 60 × instead of < 50 ×) before calling final P-element insertions with PopoolationTE2 (v1.10.03 pairupSignatures , Supplementary File 4) ( 41 ).
In a replica te popula tion, a P-element insertion is classified as 'lost', if it is present after phase 1 (i.e. generation 0), but is not detected after phase 2 (i.e. generation 10).'Persistent' P-element insertions are detectable after both phase 1 and phase 2, and 'phase-2' P-element are only detectable after phase 2. In the 'separate analysis', P-element insertions are considered to be independent across replicate popula tions (i.e.no identifica tion of orthologous insertion sites).The only exception is the r eplicate fr equency spectrum: here, P-element insertions with a distance of less than 1 kb across replicate populations are considered to be identical.Identical P-element insertions were determined with bedtools merge (v2.29.2) ( 44 ).
We observed that the joint and separate analyses resulted in qualitati v ely similar r esults.Thus, we pr esent the joint analysis in the main part of this manuscript and show the results of the separate analysis in the Supplement.
Natural population.To compare the genomic distribution of P-element insertions of our experimental populations to a natural population, we annotated previously reported Pelement insertion sites from a wild D. simulans population from South Africa with an established P-element invasion ( 4 ).For this, we obtained P-element insertion sites together with accompanying annotations and r efer ence genomes from the original analysis by Kofler et al. ( 4 , 33 ).

Characterization of P-element insertion sites
Recombination rate.We used the D. simulans recombination map from Howie et al. (file: Dsim recombination map LOESS 100kb 1.txt) ( 32 ) to check the average recombination rate in cM / Mb for each P-element insertion site in experimental populations.We used the 10th percentile ( = 2.85 cM / Mb) of the genome-wide recombina tion ra te as a threshold determining whether P-element insertions are in a region with low recombina tion ra te.

Annotations.
Genomic features.For both, experimental and natural populations, we used bedtools sort, merge, intersect and complement (v2.29.2) ( 44 ) to extract dif ferent fea tures from the D. simulans annotation ( 4 , 33 ) in order to determine the annotation type of each P-element insertion.In case different annotation tracks were ov erlapping, we recor ded the most detailed annotation (e.g.coding sequence over translated region) and the annotation with the presumably higher functional importance (e.g.exon over intron).
For each population, we calculated the expected number of P-element ( P exp ) insertions in a respecti v e genomic feature as follows: With f l being the total genomic feature length (e.g. total length of all exons, Supplementary File 5), P obs-total the total number of observed P-element insertions on the main chromosome arms, and G l the total length of the main chromosome arms.The enrichment of P-element insertions (Figure 2 A, B) is defined as the fraction of observed ( P obs ) to expected number of P-element insertions for the genomic feature of interest: P obs P exp .
We used a Cochran-Mantel-Haenszel test (CMH-test) to test, whether the abundance of lost and persistent P-element counts grouped by annota tion fea ture dif fer across experimental populations.The CMH-test allows to test of independence of stratified categorical data (i.e.P-element type abundances across replicate populations) ( 45 ).
Origin r eco gnition complex (ORC) binding sites.For experimental populations, we detected ORCs in D. simulans scrutinizing ORC sequences extracted from the D. melanogaster r efer ence genome (v5.53) ( 3 ) as described pr eviously for the natural South African population ( 4 ).Following ( 4 , 46 ), we mapped these sequences to our D. simulans r efer ence genome ( 33 ) using b wa b wasw (v0.7.17; -M ) ( 36 ).We sorted and merged overlapping ORC sequences with bedtools merge (v2.29.2) ( 44 ) and retained only ORC regions with a minimum length of 100 bp, resulting in 5788 non-overlapping ORC regions for experimental populations (Supplementary File 6).We used the R-package Ge-nomicRanges (v1.42.0) ( 47 ) to calculate the overlap in base pairs of ORCs with different annotation tracks in D. simulans ( 33 ).We note that this analysis assumes that ORCs are shared between D. melanogaster and D. simulans , but gi v en the high sequence conservation between both species ( 48 ), we do not consider that this assumption affected our analysis.Rather, if D. melanogaster ORCs are no longer functional in D. simulans , this makes our analyses conservati v e, as non-functional ORCs are included into the ORC set.
For enrichment analyses (Figure 2 B), we calculated the number of expected P-element insertions in ORCs for each replica te popula tion using equa tion ( 1 ) with f l being the total genomic length of all ORCs.For phase-2 P-element insertions f l is reduced by the total length of ORCs that are already occupied by phase-1 P-element insertions.Shared sites.We further tested whether P-element insertions are enriched in P-element insertion sites shared between D. melanogaster and D. simulans from South Africa ( 4 ) ('shared sites').For this, we extracted P-element insertions that have been previously classified as 'shared sites' ( ±1000 bp) from the D. melanogaster r efer ence genome (v5.53) ( 4 ) using bedtools getfasta (v2.29.2) ( 44 ) and mapped them with bwa bwasw (v0.7.17 -M ) ( 49 ) to our D. simulans r efer ence genome ( 33 ).We used samtools (v1.11) ( 42 ) to filter for primary alignments with a minimum mapping quality of 15.We used bedtools (v2.29.2) ( 44 ) to sort and merge overlapping r egions, r esulting in 337 regions classified as shared sites for experimental populations (Supplementary File 7).We used the R-package GenomicRanges (v1.42.0) ( 47 ) to calculate the overlap in base pairs of shared sites with different annotation features of the D. simulans annotation ( 33 ).
We used Fisher's Exact Test to assess whether shared sites have a different frequency of persistent P-elements than 3 UTR, CDS, 5 UTR, introns and intergenic regions in experimental populations (Figure 2 C).For each pairwise test (e.g.shared sites against CDS), we excluded P-elements in shared sites that are also positioned in the genomic feature of interest.We used the Benjamini-Hochberg correction to adjust p-values for multiple testing across all replicate populations.
Ag gr egation calculation.Insertion bias into ORCs / shared and selection operating against P-element insertions outside of ORCs / shared sites results in an aggregation of Pelement insertions in ORC / shared sites.To disentangle insertion bias from selection, we contrasted the aggregation of P-elements into these regions between phase 1 and phase 2. If no selection is operating, the same le v el of aggregation is expected in phase 1 and phase 2, but with selection, the observed aggregation should be stronger in phase 2. We determined the aggregation of P-element insertions in a genomic region as the fraction of P-element insertions occurring into the respecti v e region of interest (i.e.either ORCs or shared sites in our analysis).The calculation of the aggregation in phase 1 is rather straightforward: The aggregation is simply the number of P-element insertions into a genomic region divided by the total number of detected phase-1 P-element insertions.The inference of aggregation in phase 2 is a bit more complicated, because we need to account for the possibility that lost phase-1 Pelement insertions are restored by newly acquired P-element insertions during phase 2. For this, we first calculated the fraction of phase-1 P-element insertions outside ORCs and shared sites that are rediscovered after phase 2 (i.e. the fraction of persistent P-elements).This quantity --which we termed the 'rediscovery rate' --provides due to low reinsertion probabilities (caused by a lack of strong insertion bias) a good estimator for the expected fraction of persistent phase-1 P-elements in ORCs / shared sites if P-elements in these genomic regions would e volv e under the same evolutionary forces (i.e.identical selection pressure and no reinsertions due to insertion bias).The difference between the actually observed and -based on the rediscovery rate -expected number of persistent P-element insertions in ORCs / shared sites in phase 2 provides an estimate for the re-insertion rate.This calculation assumes a constant selection pr essur e on P-element insertions r egardless of their genomic position (i.e.inside / outside of ORCs / shared sites).As a next step, we estimated the number of P-element insertions that occurred during phase 2 as the sum of observed phase-2 P-element insertions and the number of expected re-insertions.Finally, the aggregation of P-elements into ORCs / shared sites during phase 2 is estimated as the fraction of P-insertions in these r egions, corr ected for reinsertions.We note that the aggregation analysis is not affected by non-uniform insertion bias among ORCs / shared sites.

Diversity measures in experimental populations
We evaluated whether base substitutions, similar to Pelements, experience less purifying selection in ORCs and shared sites.For this, we investigated the level of intraspecific polymorphisms in the ancestral experimental populations ( θ π ) and the le v el of interspecific sequence conservation ( 50 ) in ORCs and shared sites for three genomic features (5 UTR, introns, and intergenic r egions; Figur e 3; Supplementary Files 8-9).We chose 5 UTR, introns, and intergenic regions for the comparison because more than 90% of the regions annotated as ORCs or shared sites overlap with these three annotation features.

P-element proliferates in D. simulans isofemale lines
Isofemale lines were established from a natural D. simulans population, which is being invaded by the P-element ( 4 ).Lines with an acti v e P-element are expected to accumulate additional P-element copies until the piRNA defense system is established ( 57 ).We confirmed this anticipated increase in P-element copy number by sequencing three isofemale lines with an acti v e P-element ( 4) after 2 and 6 years of maintenance.We estimated P-element copy number per haplotype for each time point with DeviaTE ( 35 ).The Pelement copy number increased in all three lines to at least 11 copies per haplotype (Supplementary Figure S1, Supplementary File 1), which is broadly consistent with the distribution of P-elements in natural and experimental popula tions a t the end of a P-element invasion ( 4 , 30 ).The incr ease in cop y number confirms that the P-element has further spread in isofemale lines carrying the transposon.

Heterogeneous distribution of P-element insertions
We combined isofemale lines with and without P-elements to set up three large experimental populations with a population size of 1250 individuals each.We performed Pool-Seq for the founder populations and three replicate populations, which e volv ed for 10 generations in a hot environment fluctuating between 18 and 28 • C. We used Popoola-tionTE2 ( 41 ) to identify P-element insertion sites from Pool-Seq data of these thr ee r eplica te popula tions after phase 1 (founder popula tion, genera tion 0) and phase 2 (generation 10).We detected 1040 phase-1 P-element insertions on the main chromosome arms: 812 on the autosomes and 228 on the X-chromosome.After phase 1, we found a pronounced underr epr esentation of P-element insertions in coding regions.(joint analysis: χ 2 = 45 .524 , d f = 1 , P < 0 .001 , see Material & Methods: Genomic features , Supplementary File 2-3, see Supplementary Table S1 for replicate-specific (i.e.'separate') analysis).The 5 UTR harbored more P-element insertions than expected under a random insertion process (joint analysis: χ 2 = 129 .29 , d f = 1 , P < 0 .001 , Supplementary Table S1 (separate analysis)).Since this excess of Pelement insertions in the 5 UTR coincides with a previously reported insertion bias ( 2 ), we also measured whether the P-element is overr epr esented in origin recognition complex (ORC) binding sites --genomic regions that are considered to be P-element insertion hotspots ( 3 ).As expected from the pr eviously r eported insertion pr efer ence ( 3 ), we found a pronounced overr epr esentation of P-elements in ORCs (joint analysis: χ 2 = 303 .29 , d f = 1 , P < 0 .001 , Supplementary Table S1 (separate analysis)).

Purifying selection shapes the genome-wide P-element distribution
Phase-1 P-element insertions occurred at small effecti v e population sizes ( N e ) with low selection efficacy and were subjected to efficient selection due to larger N e in phase 2. Thus, the comparison of P-element distributions in phase 1 and phase 2 shows the impact of purifying selection on P-element insertions.We explored the fate of P-element insertions from phase 1, which either persisted from phase 1 to phase 2 ('persistent' P-element insertions) or got lost in phase 2 ('lost' P-element insertions, Figure 1 ).The abundance of persistent and lost P-elements differed significantly across genomic features (joint analysis: Cochran-Mantel-Haenszel Test, M 2 = 69 .249 , d f = 4 , P < 0 .001 , Supplementary Table S1 (separate analysis)) demonstrating a nonrandom loss of P-elements.Our analysis does not distinguish between P-element loss caused by purifying selection, random genetic drift, or excision.While selection is operating locally (i.e.varies across genomic features), the excision rate and random genetic drift should be homogeneous across the genome and thus, different genomic features are affected similarly by these random processes.Furthermore, the heterogeneous distribution of lost and persistent Pelement insertions cannot be explained by a systematic difference in P-element frequencies across genomic features in phase 1 (Kruskal-Wallis rank sum test , χ 2 = 8 .7561 , d f = 4 , P = 0 .06749 , Supplementary Figure S2A).
Gi v en this strong signal for heterogeneous selection across the genome, we determined the enrichment ( = observ ed / e xpected) of persistent P-element insertions for different genomic features (Figure 2 A).The enrichment of persistent P-element insertions in experimental populations is remar kab ly similar to enrichment of P-element insertions in a natural D. simulans population with an established P-element invasion ( 4 ): with the exception of intergenic regions and 5 UTRs (Figure 2 A), persistent P-element insertions displayed signs of purifying selection (i.e.enrichment < 1) across the genome.We observed the lowest aver age fr action of persistent P-elements in coding r egions, r eflecting the expected strong selection against TE insertions in coding sequences ( 7 , 15 ).These results suggest that the majority of P-element insertions are deleterious.Furthermore, the very efficient selection against P-element insertions implies that unlike most de novo mutations, P-element insertions ( 58 , 59 ) are most likely not recessi v e, as the low population frequency of each insertion after phase 1 makes selection against recessi v e P-element insertions in phase 2 less efficient.
The pattern for ORCs was rather unexpected.The enrichment of persistent P-element insertions in ORCs was about 85% higher than the enrichment of phase-1 P-element insertions (Figure 2 B) and resembled the enrichment of P-element insertions in ORCs of a natural D. simulans population ( 4 ).This incr eased over-r epr esentation indicates weaker purifying selection in ORCs than in the remaining genome.Yet, the enrichment analysis assumes that no new P-element insertions occurred during phase 2, which could r estor e a phase-1 P-element insertion after it was lost by drift, excision, or selection during phase 2. Hence, a strong insertion bias may be misinterpreted as a signal of low purifying selection.Because only 5.2% of the ORCs contained P-element insertions at the beginning of phase 2, such 'recurrent insertions' do not occur at a sufficiently high rate to explain the increased enrichment.Nevertheless, if the strength of insertion bias differs among ORCs, re-insertions in the same site are more likely and this analysis is not sufficient to distinguish between purifying selection and insertion bias.
With ORCs, a pr eferr ed P-element insertion target ( 2 , 3 ), experiencing less purifying selection than other genomic regions, we also evaluated another category of pr eferr ed insertion targets: sites that share P-element insertions between natural D. melanogaster and D. simulans populations ('shared sites').Because the split between D. melanogaster and D. simulans predates the P-element invasion in both species ( 4 , 28 , 60 ), the presence of shared P-element insertion sites has been attributed to a strong insertion bias of the P-element ( 4 ).We observed the same pattern for shared sites as for ORCs -the enrichment of persistent Pelements in shared sites overall agrees with the enrichment in a natural D. simulans population.The enrichment of persisting P-element insertions was e v en stronger in shared sites than in ORCs (162%, Figure 2 B).Moreover, all genomic features, e v en introns and intergenic regions, lost more phase-1 P-element insertions than shared sites (joint analysis: P adj < 0.001 for all comparisons after multiple testing correction, pairwise Fisher's exact tests between each annota tion fea tur e and shar ed sites per experimental population; Figure 2C; Supplementary Table S1 (separate analysis)).We interpret this pattern as evidence for reduced selection against P-element insertions in ORCs / shared sites.
Further evidence for low purifying selection in ORCs and shared sites comes from the probability to detect the same persistent P-element insertion in multiple replicate populations after phase 2. P-elements subjected to purifying selection are either segregating at a low frequency in e volv ed replicates or they have already been removed in a subset of the replicates.In both cases, we expect that most of these P-elements will be detected in a single replica te popula tion.Phase-1 P-element insertions exposed to strong purifying selection will be removed in all replicates during phase 2 and are classified as lost.We compared the percentage of four classes of persistent P-element insertions (within / outside ORCs / shared sites) that could be detected in one, two, and all three replicate experimental populations.Consistent with weaker purifying selection in ORCs / shared sites we observed that these genomic regions shared more persistent P-element insertions across at least two replicates than the remaining genome (Figure 2 D).The comparison of ORCs and shared sites suggests weaker purifying selection in shared sites than in ORCs (Figure 2 B-D).
The enrichment of persisting P-elements in ORCs / shared sites pro vides tw o important insights.First, the strength of selection against P-element insertions is heterogeneous across the entire genome.Second, the insertion bias of the Pelement is probably overestimated when relati v e P-element abundances of natural populations are used.
Since insertion bias and low le v els of purifying selection provide the same pattern of enrichment in ORCs and shared sites, we were interested to disentangle these two processes.In absence of reduced purifying selection in OCRs / shared sites, the aggregation of P-element insertions in ORCs / shared sites should be similar in phase 1 (with weak purifying selection due to small N e ) and phase 2 (with stronger purifying selection due to larger N e ).We estimate the aggregation of the P-element in a genomic region (i.e.ORCs / shared sites) as the fraction of P-element insertions occurring into this region.While the calculation of the le v el of aggregation in phase 1 is straightforward, estimating this quantity in phase 2 is complicated by the difficulty to distinguish truly persistent P-element insertions from P-element insertions that were lost in phase 2, but wer e r eplaced by new P-element insertions later on.Assuming that ORCs and shared sites experience the same evolutionary forces (i.e.purifying selection, drift, stochastic sampling during Pool-Seq) as other genomic regions, we calculated the 'rediscovery rate', which is the probability to detect a phase-1 P-element after phase 2. In phase 2, the total number of new P-element insertions into OCRs / shared sites can be calculated by adding the expected number of re-insertions (a quantity that can be calculated based on the rediscovery rate, see M&M section Aggregation calculations for details) to the observed number of P-element insertions.Consistent with weaker purifying selection in ORCs / shared sites, the le v el of aggr egation incr eased by 75% (joint analysis: shared sites, SEM = 29%, Supplementary Table S1 (separate analysis of each replicate)) and 43% (ORCs, SEM = 12%, Supplementary Table S1 (separate analysis of each replicate)) between phase 1 and phase 2. Assuming that the insertion pr efer ence of the P-element does not change between phase 1 and phase 2, the increase in the le v el of aggregation can be attributed to reduced purifying selection in ORCs and shared sites.

Purifying selection operates differently on TE insertions and base substitutions
We evaluated whether base substitutions, similar to P-elements, experience less purifying selection in ORCs / shared sites.We found no evidence for r elax ed purifying selection against base substitutions.Neither were intraspecific polymorphisms ( θ π ) elevated nor was interspecific sequence conservation ( 50 ) lower in ORCs and shared sites (Figure 3 ).This pattern persisted when we only considered ORCs and shared sites with persistent P-element insertions (Supplementary Figure S4 (joint analysis), Supplementary Figure S5 (separate analysis)).Interestingly, our findings suggest an opposite trend, indicating potentially enhanced purifying selection on base substitutions within ORCs / shared sites.It is important to ackno wledge, ho we v er, that our study was not explicitly designed to investigate the underlying mechanisms driving these patterns.Overall, our result suggests that for many genomic regions selection against TE insertions is dri v en by evolutionary forces tha t dif fer from selection against base substitutions.
Ectopic recombination is frequently considered a major selection force against TE insertions (9)(10)(11).To protect TE insertions in ORCs and shared sites from ectopic recombination, a lower local recombination rate would be r equir ed.Because the recombination map in D. simulans is rather uniform across almost the entire chromosome arms ( 32 , 61 ) comparison of the P-element density in genomic regions with high and low recombina tion ra tes is not very powerful.Contrasting P-element insertions in regions with low recombina tion ra tes against remaining Pelement insertions results in no difference between persistent and lost P-element insertions (joint analysis: Cochran-Mantel-Haenszel test, M 2 = 1 .2501 , d f = 1 , P = 0 .2635 , Supplementary Table S1 (separate analysis)).Furthermore, we cannot test the presence of recombination cold spots at ORCs and shared sites, because the currently available recombination map for D. simulans is too coarse ( 32 ).

DISCUSSION
Contrasting the P-element insertions in populations with low and strong selection efficacy suggested that ORCs and shar ed sites ar e not only pr eferr ed insertion targets of the Pelement, but P-element insertions in these genomic regions ar e mor e likely to be tolerated than in other parts of the genome.
One key assumption of our analyses is that properties of P-element insertions such as population frequencies, and insertion probabilities inside and outside of ORCs / shared sites were the same during phase 1.In fact, all P-element insertions occur only at a low frequency during phase 1 (Supplementary Figure S2).Ne v ertheless, P-element insertions in ORCs and shared sites have slightly higher frequencies than those outside (Supplementary Figure S2B).This higher frequency can be explained by the insertion bias, which results in independent P-element insertions at the same site in different isofemale lines.We caution, that a higher population frequency of P-element insertions in ORCs and shared sites results in a lower rate of loss of Pelement insertions due to drift which may lead to a higher re-discov ery rate, e v en without a difference in purifying selection.Howe v er, the accuracy of P-element insertion frequency estimates is constrained by the physical coverage of the data.Since low frequency estimates are particularly sensiti v e, we caution that potential differences in population frequency may not be very robust.
Figur e 3. Intraspecific pol ymorphism ( θ π ) and interspecific sequence conserva tion measurements.( A ) Average estima tes across ancestral experimental populations ( n = 3) for three different annotation features (x-axis), grouped by an overlap with ORCs.We chose 5 UTR, introns, and intergenic regions for the comparison because more than 90% of the regions annotated as ORCs or shared sites overlap with these three annota tion fea tures.( B ) Average estimates across ancestral experimental populations ( n = 3) for three dif ferent annota tion fea tures (x-axis), grouped by an overlap with shared sites ( C ) Interspecific sequence conservation for thr ee differ ent annotation features, grouped by an overlap with ORCs.( D ) Interspecific sequence conservation for thr ee differ ent annota tion fea tur es, grouped by an overlap with shar ed sites Ne v ertheless, the comparison of the enrichment in natural populations to our experiment suggests that the observed differences between ORCs / shared sites and the rest of the genome are not an artifact of the starting conditions in phase 1.The enrichment pattern for different genomic features is very similar between persistent P-elements in our laboratory experiment and P-elements in a natural D. simulans population with a fully established P-element invasion (Figure 2 A).This suggest that our experiment nicely mirrors the selecti v e forces operating on P-elements in natural populations.Importantly, we also observed a similar consistency between natural populations and our experiment for OCRs / shar ed sites (Figur e 2 B).If the small population frequency differences in phase 1 would have biased our analysis, the enrichment for ORCs / shared sites would have been more pronounced in our experiment than in natural populations.As this is not the case, we conclude that our en-richment estimates are robust to the subtle population frequency differences after phase 1.This implies that our conclusion that P-element insertions in ORCs / shared sites are subject to r elax ed purifying selection remains unaffected.
The coincidence of genomic regions with relaxed purifying selection against P-element insertions and sites with pr eviously r eported insertion pr efer ences is r emar kab le.We propose that the benefit of insertions into genomic regions with reduced selection efficacy provides a strong selection force that may have resulted in the evolution of insertion bias.Hence, insertion bias and reduced purifying selection should not be treated as two distinct phenomena, but rather be considered as cause and consequence.This hypothesis could be validated if distantly related P-elements have a different insertion bias than D. melanogaster / D. simulans Pelements.An experimental invasion of such a P-element in D .melanogaster / D .simulans would re v eal the differences

Figure 2 .
Figure 2. Characteristics of P-elements in phase 1 and phase 2. ( A ) Average enrichment ( = observ ed / e xpected) of different P-element classes by annotation feature across three replicate populations.The horizontal dashed black line marks a homogeneous distribution (i.e.no insertion bias or selection, observed = expected).Nevertheless, since we consider the enrichment in the 5 UTR to reflect reduced purifying selection in combination with insertion bias, the lower enrichment values of all other functional classes indicate purifying selection.White diamonds display the enrichment of P-elements across genomic features for a natural D. simulans population from South Africa with an established P-element invasion ( 4 ).The error bars show the standard error of the mean across experimental populations.( B ) Average enrichment of all phase-1 (red), persistent (dark red), and phase-2 (blue) P-element insertions in ORCs and shared sites.( C ) Average proportion of lost P-element insertions out of all phase-1 P-element insertions by annotation feature , ORCs , and shared sites.( D ) Replicate frequency spectrum ( 31 ) of persistent P-element insertions inside and outside of ORCs and shared sites (!ORC & !shared: P-element insertions that are located outside of OCRs and shared sites, ORC & !shared: P-elements insertions that are restricted to ORCs only, etc.).For separate analysis, please see Supplementary FigureS3.