Thomas Wicker
email: wicker@botinst.uzh.ch

Yeisoo Yu

Georg Haberer

Pradeep Klaus F X Mayer

Reddy Marri

Steve Rounsley

Mingsheng Chen

Andrea Zuccolo

Olivier Panaud

Rod A Wing

Pradeep Reddy Marri

Stefan Roffler

DNA transposon activity is associated

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

T he grass (Poaceae) family contains over 10,000 species and includes the most important agricultural crops such as rice, maize, wheat and barley. Grasses evolved from a common ancestor B70 Myr ago [START_REF] Phylogeny | Phylogeny and subfamilial classification of the grasses (Poaceae)[END_REF] . One unique characteristic of grass genomes is that they contain enormous numbers of DNA (class 2) transposons. For example, the superfamilies DTT_Mariner and DTH_Harbinger alone are present in at least 40,000 copies in grass genomes [START_REF] Schnable | The B73 maize genome: complexity, diversity, and dynamics[END_REF][START_REF]Genome sequencing and analysis of the model grass Brachypodium distachyon[END_REF] . Interestingly, the vast majority of DNA transposons in grasses are non-autonomous, meaning that they rely for their transposition on enzymes encoded by a small number of 'mother' elements elsewhere in the genome [START_REF]Genome sequencing and analysis of the model grass Brachypodium distachyon[END_REF][START_REF] Yang | A rice Tc1/mariner-like element transposes in yeast[END_REF] . Furthermore, these small non-autonomous transposons were reported to preferably insert near genes [START_REF]Genome sequencing and analysis of the model grass Brachypodium distachyon[END_REF][START_REF] Bureau | Mobile inverted-repeat elements of the Tourist family are associated with the genes of many cereal grasses[END_REF][START_REF] Bureau | Stowaway: a new family of inverted repeat elements associated with the genes of both monocotyledonous and dicotyledonous plants[END_REF] . But despite the high abundance of DNA transposons in grass genomes, little is known about their level of activity and their overall impact on genome evolution. This was mostly due to the lack of suitable data sets for comparative analyses. With the recent sequencing of 11 rice genomes in the framework of the Oryza Map Alignment Project (OMAP [START_REF] Wing | The oryza map alignment project: the golden path to unlocking the genetic potential of wild rice species[END_REF] ), data sets for such studies became available. In this study, we compared the two rice species Oryza sativa and Oryza glaberrima which diverged B600,000 years ago [START_REF] Wang | The genome sequence of African rice (Oryza glaberrima) and evidence for independent domestication[END_REF] . These two species are closely enough related to allow reliable alignment of most of the genomes and yet distant enough to have numerous transposable element (TE) polymorphisms [START_REF] Roffler | Genome-wide comparison of Asian and African rice reveals high recent activity of DNA transposons[END_REF][START_REF] Roffler | The making of a genomic parasite-the Mothra family sheds light on the evolution of Helitrons in plants[END_REF] .

DNA transposons have the curious ability to move in the genome by inserting into and excising from genomic DNA. When they excise from the genome, they leave double-strand breaks (DSBs) that have to be repaired by the cell. Previous studies have shown that, this can lead to deletions and/or insertions of 'filler' sequences at the site of the DSB [START_REF] Yang | A rice Tc1/mariner-like element transposes in yeast[END_REF][START_REF] Roffler | Genome-wide comparison of Asian and African rice reveals high recent activity of DNA transposons[END_REF][START_REF] Buchmann | Interspecies sequence comparison of Brachypodium reveals how transposon activity corrodes genome colinearity[END_REF] , depending on the repair pathway. Sometimes, re-arrangements at the excision site can be so extensive that excisions are difficult to identify [START_REF] Roffler | Genome-wide comparison of Asian and African rice reveals high recent activity of DNA transposons[END_REF][START_REF] Buchmann | Interspecies sequence comparison of Brachypodium reveals how transposon activity corrodes genome colinearity[END_REF] (Supplementary Note 1). Thus, the previous studies have established that transposons leave a variety of 'scar' patterns at the site of excision. However, DSB repair is a highly complex process that involves multiple enzymes and, in some pathways, single-stranded DNA intermediates [START_REF] Storici | Conservative repair of a chromosomal double-strand break by single-strand DNA through two steps of annealing[END_REF][START_REF] Yang | Hypermutability of damaged single-strand DNA formed at double-strand breaks and uncapped telomeres in yeast Saccharomyces cerevisiae[END_REF][START_REF] Deem | Break-induced replication is highly inaccurate[END_REF][START_REF] Fishman-Lobell | Two alternative pathways of double-strand break repair that are kinetically separable and independently modulated[END_REF][START_REF] Lydeard | Break-induced replication and telomerase-independent telomere maintenance require Pol32[END_REF][START_REF] Paull | A mechanistic basis for Mre11-directed DNA joining at microhomologies[END_REF][START_REF] Shevelev | The 3 0 -5 0 exonucleases[END_REF][START_REF] Bleuyard | Recent advances in understanding of the DNA double-strand break repair machinery of plants[END_REF][START_REF] Heacock | Molecular analysis of telomere fusions in Arabidopsis: multiple pathways for chromosome end-joining[END_REF][START_REF] West | Disruption of the Arabidopsis AtKu80 gene demonstrates an essential role for AtKu80 protein in efficient repair of DNA double-strand breaks in vivo[END_REF][START_REF] Nakagawa | Role of AtPolz, AtRev1 and AtPolZ in g ray-induced mutagenesis[END_REF] (Supplementary Note 1). Considering these complex processes, we wanted to study if and to what degree DNA transposons excisions also affect the sequences surrounding the excision site and whether they have an impact on the evolution of genes. Our data suggest that transposon excisions invoke DNA repair mechanisms that lead to high numbers of mutations around the excisions sites. The preference of DNA transposons to insert near genes in grasses therefore accelerates evolution of genes and coding regions.

Results

Transposon excisions are flanked by numerous mutations. For our analysis, we annotated 27,641 DNA transposons in the O. sativa genome; the majority of them belong to the DTT_Mariner and DTH_Harbinger superfamilies. Overall, they show a strong preference to insert close to transcription start and end points of genes (Supplementary Fig. 1). This is in agreement with previous findings that showed a preference of these elements to insert near genes [START_REF]Genome sequencing and analysis of the model grass Brachypodium distachyon[END_REF][START_REF] Bureau | Mobile inverted-repeat elements of the Tourist family are associated with the genes of many cereal grasses[END_REF][START_REF] Bureau | Stowaway: a new family of inverted repeat elements associated with the genes of both monocotyledonous and dicotyledonous plants[END_REF] (Supplementary Note 2). To identify DNA transposon polymorphisms, we compared the annotated transposon loci with their orthologs in O. glaberrima. We manually screened over 2,000 potential polymorphisms and classified 482 as insertions and 158 as excisions (Table 1; Supplementary Tables 1 and2; Supplementary Note 3). The polymorphic transposons belong to five different superfamilies of which DTT_Mariner and DTH_Harbinger elements comprise the majority (Table 1). Here, we made particular efforts to ensure that indeed orthologous loci were compared (Methods, Supplementary Fig. 2 Supplementary Note 4).

Interestingly, we found that excisions often go along with the introduction of numerous nucleotide substitutions and small insertions and deletions (InDels) in sequences flanking the transposons, with some flanking regions containing over 10 times more mutations than the genome on average (example in Fig. 1). To quantify this effect, we analysed the 12 kb flanking each polymorphic transposon and added up all nucleotide substitutions and (InDels) relative to the transposon insertion/excision site. The resulting plot shows that the overall frequency of nucleotide substitutions and InDels increases in an exponential manner towards the TE excision site to at least four-fold on average, compared with randomly picked genomic sequences (Fig. 2). Numbers of nucleotide substitutions and InDels are increased up to a distance of 3 kb from the excision point (Fig. 2). In contrast, transposon insertion sites have many fewer mutations in flanking regions, showing only a small increase in nucleotide substitution frequency in their neighbourhood (Figs 2; 3, see below).

A proposed model for how transposon excisions induce mutations. Considering findings on DSB repair from yeast [START_REF] Storici | Conservative repair of a chromosomal double-strand break by single-strand DNA through two steps of annealing[END_REF][START_REF] Yang | Hypermutability of damaged single-strand DNA formed at double-strand breaks and uncapped telomeres in yeast Saccharomyces cerevisiae[END_REF][START_REF] Deem | Break-induced replication is highly inaccurate[END_REF][START_REF] Fishman-Lobell | Two alternative pathways of double-strand break repair that are kinetically separable and independently modulated[END_REF][START_REF] Lydeard | Break-induced replication and telomerase-independent telomere maintenance require Pol32[END_REF][START_REF] Paull | A mechanistic basis for Mre11-directed DNA joining at microhomologies[END_REF][START_REF] Shevelev | The 3 0 -5 0 exonucleases[END_REF] and Arabidopsis [START_REF] Bleuyard | Recent advances in understanding of the DNA double-strand break repair machinery of plants[END_REF][START_REF] Heacock | Molecular analysis of telomere fusions in Arabidopsis: multiple pathways for chromosome end-joining[END_REF][START_REF] West | Disruption of the Arabidopsis AtKu80 gene demonstrates an essential role for AtKu80 protein in efficient repair of DNA double-strand breaks in vivo[END_REF][START_REF] Nakagawa | Role of AtPolz, AtRev1 and AtPolZ in g ray-induced mutagenesis[END_REF] , we propose a molecular mechanism that explains the high numbers of mutations flanking transposon excisions in rice (Fig. 2c): in the first step, the transposons excise from the genome, leaving a DSB for the cell to repair. After transposon excision, 3 0 overhangs are produced by exonucleases (Fig. 2c, step 1). The 3 0 overhangs then anneal using micro-homologies of a few bp (Fig. 2c, step 2), or through an intermediate generated by invasion of a foreign strand (Supplementary Note 1). Subsequently, the single-stranded DNA segments are used as templates for the synthesis of a new second strand, which is the step that introduces numerous mutations (Fig. 2c, step 3). We propose that DNA replication is analogous to that described for DSB-induced replication in yeast [START_REF] Yang | Hypermutability of damaged single-strand DNA formed at double-strand breaks and uncapped telomeres in yeast Saccharomyces cerevisiae[END_REF] . Here, mutations are introduced by translesion synthesis, possibly involving a homologue of DNA polymerase zeta (which is involved in error-prone DNA repair in yeast [START_REF] Yang | Hypermutability of damaged single-strand DNA formed at double-strand breaks and uncapped telomeres in yeast Saccharomyces cerevisiae[END_REF] ) and by a DSB-induced replication complex that has deficiencies in DNA polymerase delta fidelity and mismatch repair, analogous to that described in yeast [START_REF] Deem | Break-induced replication is highly inaccurate[END_REF] . Possibly, Rev1 polymerase also contributes to erroneous DNA repair [START_REF] Nakagawa | Role of AtPolz, AtRev1 and AtPolZ in g ray-induced mutagenesis[END_REF] . The end product of the repair process are sequence segments flanking the transposon excision which are riddled with nucleotide substitutions and small InDels (Fig. 1 and Fig. 2c, step 4). The length of the segment containing the mutations depends on the size of the 3 0 overhang produced in the initial repair step. In yeast, these overhangs can be several kb in size [START_REF] Storici | Conservative repair of a chromosomal double-strand break by single-strand DNA through two steps of annealing[END_REF][START_REF] Yang | Hypermutability of damaged single-strand DNA formed at double-strand breaks and uncapped telomeres in yeast Saccharomyces cerevisiae[END_REF][START_REF] Fishman-Lobell | Two alternative pathways of double-strand break repair that are kinetically separable and independently modulated[END_REF] , and this is expected to be similar in plants, due to the high conservation of DSB repair pathways [START_REF] Shevelev | The 3 0 -5 0 exonucleases[END_REF] . Indeed, our data support this notion, since the observed average nucleotide substitution frequency levels off B3 kb away from the excision site (Fig. 2a,b). Step 2: the 3 0 overhangs anneal using microhomologies. To keep it simple, we only represent single-strand annealing (SSA [START_REF] Sampath | Genome-wide comparative analysis of 20 miniature inverted-repeat transposable element families in Brassica rapa and B. oleracea[END_REF][START_REF] Ziolkowski | Genome sequence comparison of Col and Ler lines reveals the dynamic nature of Arabidopsis chromosomes[END_REF] ) here. Alternatively, the strands could also be connected via synthesis-dependent strand annealing (SDSA [START_REF] Sampath | Genome-wide comparative analysis of 20 miniature inverted-repeat transposable element families in Brassica rapa and B. oleracea[END_REF][START_REF] Ziolkowski | Genome sequence comparison of Col and Ler lines reveals the dynamic nature of Arabidopsis chromosomes[END_REF][START_REF] Joly-Lopez | Diversity and evolution of transposable elements in Arabidopsis[END_REF] ), where the two strands are connected by 'filler' sequences (which were found in some cases, not shown).

Step 3: new strands are synthesized by a replication complex that has deficiencies in DNA polymerases fidelity and mismatch repair.

Step 4: the final repair product is rich in nucleotide substitutions and small insertions and deletions.

TE insertions suggest repair patterns similar to excisions. Interestingly, we also found a slight increase in the number of mutations close to TE insertion sites (Figs 2,3). The fact that we observed this for DNA transposons as well as for retrotransposons suggest that the underlying molecular mechanism may be the same for both classes (Fig. 3). When TEs insert into the genome, they produce a staggered cut with 5' overhangs [START_REF] Van Luenen | The mechanism of transposition of Tc3 in C. elegans[END_REF] . The insertion therefore results in an intermediate where the TE is ligated to short single-stranded segments, the subsequent repair of which produces the TSD (Supplementary Fig. 3). We propose that this intermediate can, in some cases, become the target of 3'-5' exonucleases which expose longer segments of single-stranded DNA (Fig. 3b). Repair of these single-stranded stretches would then engage the same error-prone replication complex as proposed for transposon excisions (Fig. 2c). However, the proposed model would then also require that TSDs themselves should, in many cases, not be perfect repeats, but contain more substitutions than would be expected from the overall mutation rate of the genome. We tested this hypothesis by analysing insertion sites of 192 long terminal repeats (LTRs) retrotransposons from three different families in O. sativa (Supplementary Table 3). Due to their replication mechanism, the two LTRs at the ends of the retrotransposon are identical at the time of insertion. The 'age' of a retrotransposon can therefore be estimated based on the differences the LTRs have accumulated over time [START_REF] Sanmiguel | The paleontology of intergene retrotransposons of maize[END_REF] (Supplementary Fig. 3b). By comparing substitution rates in LTRs with those in TSDs, we found that TSDs contain on average almost five times more substitutions than LTRs (Supplementary Table 3). These data suggest that second strand synthesis following a TE insertion is carried out by the same error-prone replication complex as proposed for excisions (Fig. 2), but that the single-stranded segments are on average either shorter or produced only in rare cases.

Excisions associate with elevated mutation rates in genes.

Because DNA transposons preferably reside in gene promoters, we expected that these regions should evolve at a particularly high rate. Indeed, we found that the 2,000 bp upstream of genes consistently contain 20-29% more nucleotide substitutions than intergenic sequences from the same chromosomal region (Fig. 4; Table 2). Because the genomes of the closely related O. sativa and O. glaberrima are B99.5% identical on average, the differences in sequence conservation between promoters and intergenic sequences are small, but the large sample size assures that they are highly significant (P value o2.2E À 16). Intergenic regions in rice are mostly comprised of class 1 retrotransposons which are believed to be largely free from selection pressure. It is therefore intriguing that DNA repair following transposon excisions In both cases, nucleotide substitution frequency increases slightly towards the insertion point. This indicates that insertions are also associated with small numbers of mutations in their flanking sequences. Furthermore, this result is evidence that events classified as DNA transposon insertions probably do not contain many precise excisions. (b) Proposed mechanism for errorprone DNA repair following TE insertions (see also Supplementary Fig. 3).

Step 1: the TE inserts into the genome by producing a staggered cut, resulting in a TE that is ligated to the genomic DNA via single-stranded segments.

Step 2: 3 0 -5 0 exonucleases expose large stretches of single-stranded DNA.

Step 3: second strands are synthesized by a replication complex that has deficiencies in DNA polymerases fidelity and mismatch repair (the same as described in Fig. 2c).

Step 4: the final TE insertion is flanked by segments rich in nucleotide substitutions and small insertions and deletions.

Promoter regions Intergenic regions

1/3 1/3 1/3 1/3 1/3 1/3 Bin1 Bin2 Bin3 Bin4 Bin5
Size-normalized chromosome Promoters are on average 20-29% less conserved than intergenic sequences from the same chromosome bin. To calculate sequence conservation in intergenic regions, we isolated segments that are located in the middle of intergenic sequences which are at least 10 kb in size (that is, the distance between the end of one gene and the start of the next one is over 10 kb).

apparently leads to increased mutation rates of promoters to a degree that they evolve more rapidly than selectively neutral sequences. Interestingly, sequence conservation is generally lower in the centromeric and pericentromeric regions of chromosomes than in distal regions (Fig. 4), for which we have no explanation at this point.

The preference of DNA transposons to reside in up-and downstream regions of genes also implied that the 5 0 and 3 0 ends of coding sequences (CDS) should show an overall higher substitution rate than their central parts. Thus, we aligned CDS of closest homologues from O. sativa and O. glaberrima and studied overall sequence conservation as well as distributions of nucleotide substitutions along the aligned CDS. Overall, most CDS from O. sativa and O. glaberrima CDS are 499.5% identical. However, the distribution of sequence identities trails off with some CDS being o97% identical (Supplementary Fig. 4). We expected that CDS which are 499.5% identical have not experienced transposon excisions in their vicinity, while genes with lower sequence identity could be those that have accumulated mutations due to a nearby transposon excisions. Indeed, we found that genes with lower than median sequence identity ranging from 98 to 99.4% show a 427% higher number of substitutions in their 5 0 and 3 0 regions than in the central part of the CDS (Fig. 5a; Supplementary Table 4), while genes with higher levels of sequence conservation do not show this pattern (Supplementary Fig. 5). Here, we only considered nucleotide substitutions in synonymous sites to exclude effects of differing selection pressures in different parts of the genes.

SNP accumulations predict the presence of excision sites. Since we predict that DNA repair following transposons excision is responsible for high numbers of mutations in their flanking sequences, regions containing above average numbers of mutations should, in turn, often contain transposon excision sites. Thus, we inspected sequence alignments from O. sativa and O. glaberrima that covered genes plus 3 kb of their flanking regions, and selected 50 segments that contained regions with local SNP accumulations ('high-SNP' set, examples in Supplementary Fig. 6). As a control, 50 segments with an overall low SNP density, similar to that of the genome-wide average were used (examples in Supplementary Fig. 6). The 100 alignments were manually searched in detail for the presence of polymorphic transposons and other insertions and deletions (InDels). In the high-SNP data set, we identified 16 TE excisions, while in the control data set we only identified two excisions, a highly significant enrichment (Table 3). Interestingly, the high-SNP data set was also significantly depleted in transposons insertions, with 16 insertions identified in the high-SNP and 27 in the control data set (Table 3). This complements the above findings that transposon insertions only in rare cases are associated with SNP accumulations in their flanking regions (Fig. 3). We also surveyed InDels and repeat slippages (that is, differences in numbers of repeat units in micro-and minisatellites), since they can also result from DSB repair and thus could also be responsible of the introduction of SNPs. Here, we found no significant differences between the high-SNP and control data set. Although there are obviously several different causes for SNP accumulations, we identified transposon excisions as likely the main difference between regions that contain high numbers of SNPs and those which do not. Thus, these data show that local SNP accumulations can be used as search criterion for the identification of TE excisions.

Increased mutation rates in genes are common in grasses.

Because all grass genomes sequenced so far are rich in DNA transposons, we predicted that we would find increased mutation rates also in genes from other grasses. We therefore compared closest gene homologues from wheat and barley, two species which diverged B8 Myr ago [START_REF] Middleton | Sequencing of chloroplast genomes from wheat, barley, rye and their relatives provides a detailed insight into the evolution of the Triticeae tribe[END_REF] . Indeed, the 5 0 and 3 0 regions of the genes show a 420% higher number of substitutions than the central part of the genes (Fig. 5b). We also analysed maize where many genes are present in duplicates because maize is a relatively young polyploid that underwent a whole-genome duplication 5-10 Myr ago 2 . Thus, a comparison of such intragenomic closest homologues is analogous to a comparison of genes between two species. Here we found an even stronger effect, with 5 0 and 3 0 regions showing almost 30% more substitutions than the central part of genes (Fig. 5c). For both, the wheat/barley and the maize intragenomic CDS comparisons, the effects are statistically highly significantly (Supplementary Table 4). Considering that rice, maize, wheat and barley represent three different major clades of the grasses, our data strongly indicate that the described higher mutation rates in genes and regulatory sequences is common to all grasses. Interestingly, we did not find elevated mutation rates in genes in representatives of dicotyledonous plants ('dicots') such as Arabidopsis, Brassica, poplar and soybean (example in Fig. 5d; Supplementary Fig. 7; Supplementary Note 5). A de novo search for class 2 elements in these dicot genomes revealed that they contain at least 10-100 times fewer small DNA transposons than grasses (Supplementary Fig. 8; Supplementary Note 5). This result is in agreement with recent findings [START_REF] Sampath | Genome-wide comparative analysis of 20 miniature inverted-repeat transposable element families in Brassica rapa and B. oleracea[END_REF] . Furthermore, DNA transposons in Arabidopsis were found to be similarly active to those from rice [START_REF] Ziolkowski | Genome sequence comparison of Col and Ler lines reveals the dynamic nature of Arabidopsis chromosomes[END_REF][START_REF] Joly-Lopez | Diversity and evolution of transposable elements in Arabidopsis[END_REF] , but their much lower numbers may diminish their impact, even if they have the same mutagenic effect per individual transposition event. Thus, these data strengthen the correlation even more between the presence of DNA transposons and increased mutation rates of genes. 29-31). So far, most widely accepted is their role in altering gene expression. For example, a TE-mediated increase in expression level of the tb1 gene in maize resulted in plants with fewer branches, a fundamental step in maize domestication [START_REF] Tsiantis | A transposon in tb1 drove maize domestication[END_REF][START_REF] Clark | A distant upstream enhancer at the maize domestication gene tb1 has pleiotropic effects on plant and inflorescent architecture[END_REF] . We did indeed find that the presence of transposons is associated with higher levels of DNA methylation sites, suggesting an effect on transcription (Supplementary Fig. 9; Supplementary Note 6). However, the main contrast to previous studies is that our data show that transposon activity is associated with higher mutation rate and therefore may directly change coding sequences and regulatory regions by introducing nucleotide substitutions and InDels during DNA repair. We propose that error-prone repair of excision sites can introduce many mutations hundreds or even thousands of base pairs away from the sites. This would have the profound result that, even if the excision changes only a few base pairs at the actual transposon site [START_REF] Yang | A rice Tc1/mariner-like element transposes in yeast[END_REF][START_REF] Roffler | Genome-wide comparison of Asian and African rice reveals high recent activity of DNA transposons[END_REF][START_REF] Buchmann | Interspecies sequence comparison of Brachypodium reveals how transposon activity corrodes genome colinearity[END_REF] , the entire genomic region accumulates mutations as a result of error-prone strand synthesis (Fig. 2; Supplementary Note 7). Most importantly, we show that this could affect thousands of genes in the species studied, and we provide evidence that this phenomenon is common to the vast family of the grasses with its over 10,000 species. Our data thereby also indicate that the highly successful types of non-autonomous DNA transposon elements that are associated with higher mutation rate and could therefore drive the accelerated evolution of genes only evolved after the separation of monocotyledons and dicotyledons B145-300 Myr ago [START_REF] Kawai | The deep phylogeny of land plants inferred from a full analysis of nucleotide base changes in terms of mutation and selection[END_REF][START_REF] Zimmer | Dating the early evolution of plants: detection and molecular clock analyses of orthologs[END_REF] . We previously showed that about 3% of the DNA transposons in rice have moved within the past 600,000 years, indicating that these elements are highly active [START_REF] Roffler | Genome-wide comparison of Asian and African rice reveals high recent activity of DNA transposons[END_REF] . Since DNA transposons are present in tens of thousands of copies in grasses [START_REF] Schnable | The B73 maize genome: complexity, diversity, and dynamics[END_REF][START_REF]Genome sequencing and analysis of the model grass Brachypodium distachyon[END_REF] , most genes will experience transposon excisions in their proximity at some point and therefore may accumulate particularly high numbers of mutations over time. Consequently, this may explain the stronger mutation rate gradient we found in more distantly related grasses such as wheat and barley (Fig. 5).

Discussion

Data on how TEs contribute to gene evolution has been somewhat anecdotal (examples in refs

In plants and animals, a dominant DSB repair pathway is nonhomologous end joining (NHEJ), where broken ends are directly joined, leading often to small deletions or insertions of 'filler' sequences [START_REF] Bleuyard | Recent advances in understanding of the DNA double-strand break repair machinery of plants[END_REF][START_REF] West | Disruption of the Arabidopsis AtKu80 gene demonstrates an essential role for AtKu80 protein in efficient repair of DNA double-strand breaks in vivo[END_REF] . Thus, NHEJ can explain certain repair patterns that were previously found at the immediate site of transposon excisions [START_REF] Yang | A rice Tc1/mariner-like element transposes in yeast[END_REF][START_REF] Roffler | Genome-wide comparison of Asian and African rice reveals high recent activity of DNA transposons[END_REF][START_REF] Buchmann | Interspecies sequence comparison of Brachypodium reveals how transposon activity corrodes genome colinearity[END_REF] . However, NHEJ does not require processing of the broken ends into single-stranded DNA. But our data strongly suggests that the repair pathway must involve single-stranded intermediates. Thus, our models are based on other known repair pathways. For this, we rely strongly on findings in yeast, where DNA repair processes are extremely well studied. We consider this legitimate, since most DSB repair pathways were probably established very early in eukaryote evolution. Indeed, practically all genes involved in DSB repair in yeast have homologues in plants, suggesting that DSB repair processes are virtually identical in plants and fungi [START_REF] Bleuyard | Recent advances in understanding of the DNA double-strand break repair machinery of plants[END_REF][START_REF] Heacock | Molecular analysis of telomere fusions in Arabidopsis: multiple pathways for chromosome end-joining[END_REF][START_REF] West | Disruption of the Arabidopsis AtKu80 gene demonstrates an essential role for AtKu80 protein in efficient repair of DNA double-strand breaks in vivo[END_REF] . Furthermore, studies on Arabidopsis mutants showed that many of these genes are involved in the same processes as in yeast [START_REF] Bleuyard | Recent advances in understanding of the DNA double-strand break repair machinery of plants[END_REF] . For example, the yeast genes Mre11, Rad50 and Xrs2 which are required for micro-homology mediated end joining (the type on which our models are based) were shown to be involved in the same processes in Arabidopsis 20 . These findings are especially relevant for our model of DNA repair following TE insertions (which requires replication-independent 3 0 -5 0 exonucleases for the extension of single-stranded regions), because the Mre11 exonuclease The sequences of the test set were chosen based on the presence of regions with high numbers of SNPs. Gene-containing regions that had a SNP density similar to that of the genome overall served as a control. Differences between test and control sets were tested with a w 2 -test.

produces single-stranded DNA intermediates during DSB repair in yeast [START_REF] Paull | A mechanistic basis for Mre11-directed DNA joining at microhomologies[END_REF] .

We propose that the activity of DNA transposons is a major driving force in the evolution of grasses, because DNA repair following transposon excisions may specifically accelerate evolution of genes. Our findings may, in part, explain the phenomenal evolutionary success of the grasses, a very large group of plants that contains the most important crops such as rice, maize, wheat, sorghum and barley which are the basis of most food consumed by humankind.

Methods

Survey of DNA transposon distribution relative to rice genes. A total of 101 sequences of DTT_Mariner and DTH_Harbinger transposons from rice were obtained from the TREP database (wheat.pw.usda.gov/ITMI/Repeats/). They represent 19 DTT_Mariner and 25 DTH_Harbinger families. The 101 sequences were mapped with blastn to the O. sativa genome (version 6) using an in-house Perl script. The cutoff for blast hits was 50 bp and 80% sequence identity. If multiple TE families mapped to the same location, the one with the strongest blastn hit was chosen. To analyse their position relative to genes, the TE annotation was then cross-matched with the gff format gene annotation of the rice genome. We used the annotated transcription start and end points as anchor points and generated a data set of the positions of all annotated TEs within 5 kb upstream of the transcription start point and 5 kb downstream of the transcription end point for each gene. Furthermore, positions of TEs inside the gene were recorded. We selected genes larger than 4 kb and recorded TE positions within 2 kb from each end of a gene. For simplicity, only genes in forward orientation were used. The final dataset included data for 4,994 genes. Sequences covered by TEs were added up for all genes, resulting in a final coverage plot that reflects the overall distribution of TEs relative to genes (Supplementary Fig. 1).

Identification of transposon polymorphisms. We used an alignment of B60% of the genomes of O. sativa and O. glaberrima described in our previous study [START_REF] Roffler | Genome-wide comparison of Asian and African rice reveals high recent activity of DNA transposons[END_REF] to identify insertions larger than 50 bp. Insertions were screened for homology with TE sequences by blastn against the TREP database (wheat.pw.usda.gov/ITMI/ Repeats/). Using an in-house Perl script, TEs with the highest homology were mapped onto the O. sativa/O. glaberrima alignments to facilitate visual inspection and to classify the polymorphism as transposon insertion or excision. Over 2,000 polymorphisms were screened, yielding the 482 insertions and 158 excisions (Table 1; Supplementary Tables 1 and2).

Test for orthology of the analysed loci. To ensure that the aligned sequences from O. sativa and O. glaberrima indeed come from orthologous loci, we mapped the sequences used for the alignments back onto both genomes. That is, the sequences from O. sativa were first mapped back to the O. sativa genome and then mapped on to the O. glaberrima genome. The same was done vice versa with the corresponding O. glaberrima sequence. We split the aligned 24 kb regions into segments of 1,000 bp and mapped each segment by blastn to the genome it came from as well as to the genome of the other species. This was done because blast alignments are often fragmented due to the presence of low-complexity sequences or TE insertions in one or the other species. Therefore, one cannot expect a long sequence from one species to produce a similarly long blast hit in another. We therefore rather assigned each locus a score for how many of the segments map in the putative orthologous region in the other genome as a quantitative assessment of how strong the evidence for true orthology is for a particular locus.

For each 1,000 bp segment, we recorded the positions of the top blast hit in the genome it came from as well as to the genome of the other species. We required that the top blast hit produced an alignment of at least 600 bp. Thus, some segments could not be mapped due to the presence of low-complexity sequences that are filtered out in the blastn search. Furthermore, one expects that not all segments map unambiguously to the orthologous locus in the other genome. This can, for example, be due to a large retrotransposon insertion in one species. The segments covering that retrotransposon would have no counterpart in the orthologous locus in the other species and therefore map elsewhere in the genome. The genomic region where the majority of the segments map was considered the putative ortholog. Furthermore, since we ran the analysis in both directions, we required that sequences from both species had to identify each other as the closest homologue. All analysed loci fulfilled these criteria. Additionally, as Supplementary Fig. 2 shows, all except two loci are located in perfect colinear order along the chromosomes.

Distinguishing transposon insertions and excisions. We defined a TE polymorphism as an insertion if one species contained the TE plus the duplicated target site (TSD) on both sides, while the other species only contained one copy of the target site. Excision are more difficult to define as they can go along with various re-arrangements [START_REF] Roffler | Genome-wide comparison of Asian and African rice reveals high recent activity of DNA transposons[END_REF][START_REF] Buchmann | Interspecies sequence comparison of Brachypodium reveals how transposon activity corrodes genome colinearity[END_REF] . In general, we defined an excision by the absence of the TE in one species, with the pattern differing from that of an insertion. We distinguished different types of excisions: (i) in a perfect excision, as previously defined [START_REF] Buchmann | Interspecies sequence comparison of Brachypodium reveals how transposon activity corrodes genome colinearity[END_REF] , one species contains the TE with flanked by the two units of the TSD while the other species does not contain the TE but both copies of the TSD. (ii) Excisions with deletions were defined as the TE plus some flanking sequences being absent in one species. To distinguish these events from random deletions that by chance removed the TE plus flanking regions, we requested that one breakpoint of the excision be within 3 bp of one end of the TE (we considered it unlikely that a random deletion would have one of its borders so close to the end of a TE). (iii) Excisions with fillers were defined as events where the TE in one species is replaced with a completely unrelated sequence in the other. Fillers can range from a few bp to several kb. Also here, we requested that end of the filler sequence be within 3 bp of one end of the TE. Filler insertions were often found combined with deletions as described in (ii). Additional methodological considerations on distinguishing transposon excisions from insertions are provided in Supplementary Note 3.

Quantification of mutations flanking polymorphic TEs. For all identified insertions and excisions, 12 kb of the flanking sequences were extracted from the O. sativa and O. glaberrima genome-wide alignment. We selected all alignments where 47,000 bases could be aligned (due to large insertions and deletions and/or colinearity breaks, usually o12 kb were actually aligned). This selection resulted in 206 sequence alignments for excisions and 438 for insertions. The transposon excision/insertion site was used as anchor point (that is, position zero) from which all nucleotide substitutions and InDels were recorded. Sequence polymorphisms were added up for all alignments relative to the TE excision/insertion site. For the graphical representation (Fig. 2a,b), nucleotide substitutions and InDel densities were calculated by a running average. Survey of LTR retrotransposon insertions. Consensus sequences of LTRs from the O. sativa retrotransposon families RLG_Cara, RLG_Houba and RLG_hopi were used in blastn searches against the O. sativa genome. LTRs of the same family wich were found in the same direction and o14 kb apart were considered candidates for full-length elements. These including the 5 bp flanking sequences (corresponding to the TSD) were extracted from the genome. All candidate elements were visually inspected by DotPlot against a reference sequence of the respective retrotransposon family, to ensure that indeed full-length elements were selected (instead of, for example, two solo-LTRs that just happen to be located near each other). All LTR pairs of the individual copies were aligned with the programme WATER (emboss package, emboss.sourceforge.net/) to determine the number of substitutions between LTRs. From this, the average sequence conservation of the LTRs for each retrotransposon family was calculated (we excluded LTR pairs where sequence homology was over two standard deviations lower than that of the entire family, since such events could be results inter-element recombination). Analogously, the TSD sequences of all copies were aligned. The total number of mismatches in TSDs was then compared to that in LTRs. A w 2 -test was used to test if the two values differed from each other (Supplementary Table 3).

Comparison of promoters from O. sativa and O. glaberrima. Information on start and end point of genes was extracted from the gff format annotation of the rice genome. As start and end point of genes we used transcription start and end points. Here, we used rice genome version 5, because our previously published genome alignment of O. sativa and O. glaberrima 9 was done with this version. We defined the region from the transcription start point to 2 kb upstream of it as promoter region. Alignments were accepted when 4600 bp in this 2 kb region could be aligned between O. sativa and O. glaberrima. For comparison, alignments of intergenic sequences were used. Here, we isolated segments that are located in the middle of intergenic sequence that are at least 10 kb in size (that is, the distance between the end of one gene and the start of the next one is over 10 kb). Because sequence conservation along chromosome varies (Fig. 3), chromosome arms were divided into three equally sized bins for comparison of promoter and intergenic sequences. Data for promoters and intergenic sequences were analysed separately for each chromosome bin. To test whether the data sets for the individual bins differ from each other, the wilcox.test programme from RStudio (rstudio.com) was used.

Comparison of CDS of genes. Repositories where CDS of different species were obtained are listed in Supplementary Table 5. CDS for O. glaberrima were deduced from aliment with O. sativa CDS and are available upon request. Closest homologues from different species or, in the case of maize, homeologs that originated from a whole-genome duplication were identified by bi-directional blastn searches. Only homologues which had each other as the top blastn hit were used for comparison. Bi-directional closest homologues were aligned at the protein level using the programme WATER from the EMBOSS package (emboss.sourceforge.net). The aligned protein sequences were back-translated to ensure that corresponding codons were aligned. We considered only alignment positions corresponding to the third codon base for Ala, Gly, Leu, Pro, Arg, Ser, Thr and Val. For those amino acids which all have six possible codons (Leu, Arg and Ser), we used only the codons starting with CT, TC and CG, respectively (that is, the codons in which the third base can be exchanged without causing an amino acid change). To normalize the different sizes of genes, the aligned CDS were split into five equally sized bins. To obtain sufficiently high numbers of synonymous substitutions, we used only gene pairs where 41,500 bp of the CDS could be aligned. For each bin of each gene, we calculated the number of synonymous substitutions per kb. Finally, we compiled the data for the five bins for all genes. To test whether the data sets for the individual bins differ from each other, the wilcox.test programme from RStudio was used.

De novo identification of small DNA transposons in dicots. DNA transposons are characterized by the presence of terminal inverted repeats which serve as binding site for transposase enzymes [START_REF] Wicker | A unified classification system for eukaryotic transposable elements[END_REF] . The initial step of de novo identification was to screen chromosomal segments in windows of 1,000 bp, which overlap by 500 bp. The 1,000 windows were aligned with the programme WATER from the EMBOSS package against themselves in reverse orientation. Outputs were parsed and visually inspected for the presence of inverted repeats longer than B15 bp and over B70% identity. The candidate sequences (inverted repeat and the sequences between them) were excised from the 1,000 bp. The candidate TEs were then used in blastn searches against the respective genome. Sequences with multiple hits were considered true DNA transposons. The de novo detection was done on one entire Arabidopsis chromosome, 2 Mbp of poplar linkage group 1 and 500 kb of rice chromosome 10 (Supplementary Fig. 8).

Comparative analysis of DNA methylation. Data on methylation sites in O. sativa and O. glaberrima were kindly provided by Detlef Weigel and Claude Becker (Max Planck Institute for Developmental Biology, Tu ¨bingen, Germany). These data sets will be published elsewhere (personal communication, Detlef Weigel and Claude Becker). Sequence segments of 4 kb spanning the polymorphic transposon in O. sativa and O. glaberrima were extracted from the chromosomes. Methylated sites were flagged and the sequence segments were aligned with the programme Water (emboss package, emboss.sourceforge.net/). Since we found that practically no methylation sites were conserved between the two species, methylation states were compared by simply counting the numbers of methylated sites in the sequences segments from the two species. The ratio of the number of methylation sites in O. sativa and O. glaberrima was then calculated for each transposon locus. For comparison, a second segment 2,000-4,000 bp downstream of the transposon was extracted.

Statistics. Wilcoxon rank sum test was used to test whether substitution rates in different bins of size-normalized genes differ from each other. Sample sizes depended on how many bi-directional closest homologues could be identified between species. Sample sizes are provided in Fig. 5. To test if SNP accumulations can be used to predict transposon excisions, results of 50 candidate sequences were compared with those of 50 control sequences. The sample size of 50 was used to meet the commonly used small sample size criteria. A w 2 -test was used to test for differences between test and control sets. To test if substitution rates in TSDs differ from those in LTRs, 192 full-length LTR retrotransposons were isolated from the rice genome. Sample size was determined by the copy number of retrotransposons. A w 2 -test was used to test for differences between substitution rates in TSDs and LTRs.

Data availability. Repositories where CDS of different species were obtained are listed in Supplementary Table 5. The genome sequence of O. glaberrima can be obtained from Gramene (ensembl.gramene.org). The authors declare that all other data supporting the findings of this study are available within the manuscript and its Supplementary Information Files or are available from the corresponding author upon request (such as original software and sequence aliments of genomic and CDS sequences).

Figure 1 |

 1 Figure 1 | Example for a DNA transposon excision with numerous nucleotide substitutions in its flanking region. A DTH_Harbinger transposon was excised from the genome of O. sativa (Osat) while it remained present in O. glaberrima (Ogla). In this particular event, the transposon was excised almost perfectly, only losing 2 bp of the target site duplication and replacing one of them with a mismatching base. The 211 bp upstream and 120 bp downstream of the excision contain 25 mutations (InDels 41 bp are counted as one mutation), resulting in o93% sequence identity and thus making the mutation rate over 15 times higher than for the genome overall. Outside of the region with the mutations, O. sativa and O. glaberrima sequences are identical, reflecting the overall genome-wide sequence conservation of B99.5%. The segments shown correspond to O. sativa chromosome 1 position 23,814,561-23,815,081 and O. glaberrima chromosome 1 position 16,579,166-16,580,116.

Figure 2 |

 2 Figure 2 | Mutations in sequences flanking transposon insertion and excision sites. (a) Frequencies of nucleotide substitutions relative to transposon insertion/excision sites in rice. For the plot, 438 sequence alignments carrying transposon insertions (blue line) and 206 alignments carrying excisions (red line) were compiled. As control, 340 alignments of randomly picked orthologous sequences from O. sativa and O. glaberrima were used (grey line, see methods). Nucleotide substitution frequencies were calculated in a 400 bp sliding window with a 40 bp sliding step. (b) Insertion/Deletion (InDel) frequency calculated in a 1,000 bp sliding window with a 100 bp sliding step. (c) Proposed mechanism for error-prone DNA repair following the excision of DNA transposons. Step 1: after transposable element (TE) excision, 3 0 overhangs are generated by exonuclease (blue).Step 2: the 3 0 overhangs anneal using microhomologies. To keep it simple, we only represent single-strand annealing (SSA[START_REF] Sampath | Genome-wide comparative analysis of 20 miniature inverted-repeat transposable element families in Brassica rapa and B. oleracea[END_REF][START_REF] Ziolkowski | Genome sequence comparison of Col and Ler lines reveals the dynamic nature of Arabidopsis chromosomes[END_REF] ) here. Alternatively, the strands could also be connected via synthesis-dependent strand annealing (SDSA[START_REF] Sampath | Genome-wide comparative analysis of 20 miniature inverted-repeat transposable element families in Brassica rapa and B. oleracea[END_REF][START_REF] Ziolkowski | Genome sequence comparison of Col and Ler lines reveals the dynamic nature of Arabidopsis chromosomes[END_REF][START_REF] Joly-Lopez | Diversity and evolution of transposable elements in Arabidopsis[END_REF] ), where the two strands are connected by 'filler' sequences (which were found in some cases, not shown). Step 3: new strands are synthesized by a replication complex that has deficiencies in DNA polymerases fidelity and mismatch repair.Step 4: the final repair product is rich in nucleotide substitutions and small insertions and deletions.

Figure 3 |

 3 Figure 3 | Mutations in sequences flanking transposable element insertions. (a) Frequencies of nucleotide substitutions relative to insertion sites of DNA transposons and Gypsy retrotransposons sites in rice.In both cases, nucleotide substitution frequency increases slightly towards the insertion point. This indicates that insertions are also associated with small numbers of mutations in their flanking sequences. Furthermore, this result is evidence that events classified as DNA transposon insertions probably do not contain many precise excisions. (b) Proposed mechanism for errorprone DNA repair following TE insertions (see also Supplementary Fig.3).Step 1: the TE inserts into the genome by producing a staggered cut, resulting in a TE that is ligated to the genomic DNA via single-stranded segments.Step 2: 3 0 -5 0 exonucleases expose large stretches of single-stranded DNA. Step 3: second strands are synthesized by a replication complex that has deficiencies in DNA polymerases fidelity and mismatch repair (the same as described in Fig.2c). Step 4: the final TE insertion is flanked by segments rich in nucleotide substitutions and small insertions and deletions.

Figure 4 |

 4 Figure 4 | Sequence conservation along Oryza sativa and Oryza glaberrima chromosomes. Data from all 12 chromosomes were compiled and chromosome sizes were normalized by dividing chromosome arms into three equally sized bins (x-axis). The y-axis depicts sequence identity of orthologous sequences. For each chromosome bin, promoter regions (the 2,000 bp upstream of the transcription start point, red box plots) are compared with intergenic sequences from the same bin (blue box plots).Promoters are on average 20-29% less conserved than intergenic sequences from the same chromosome bin. To calculate sequence conservation in intergenic regions, we isolated segments that are located in the middle of intergenic sequences which are at least 10 kb in size (that is, the distance between the end of one gene and the start of the next one is over 10 kb).

Figure 5 |

 5 Figure 5 | Substitution frequencies in synonymous sites showing that grass genes have higher mutation rates in their 5 0 and 3 0 regions. To normalize the different CDS sizes, genes were divided into five equally sized bins and frequencies were normalized to nucleotide substitutions per kb for each bin. The bold line inside the box is the median value, while mean values are indicated with numbers. (a) Comparison of 442 closest homologues from O. sativa and O. glaberrima. (b) Comparison of 2,314 pairs of closest homologues from wheat and barley. (c) Comparison of 428 pairs of intragenomic closest homologues in maize that originate from a whole-genome duplication. (d) Comparison of 4,133 pairs of closest homologues from A. thaliana and A. lyrata.

Table 1 |

 1 Transposition events identified and manually curated in the comparison of the two rice species O. sativa and O. glaberrima.

	Superfamily	Insertions	Excisions
	DTH_Harbinger	241	71
	DTT_Mariner	137	64
	DTM_Mutator	77	20
	DTA_hAT	23	1
	DTC_CACTA	4	2
	Total	482	158

Table 2 |

 2 Mean sequence conservation of promoter and intergenic sequences in different chromosome bins of O. sativa and O. glaberrima.

	Chromosome bin	Promoter	Random	Difference* (%)
	1	98.22	98.62	28.99
	2	98.03	98.35	19.39
	3	97.8	98.17	20.22
	4	98.06	98.39	20.50
	5	98.33	98.58	17.61

*Difference in sequence divergence between promoter and intergenic (random) sequences.

Table 3 |

 3 Test for predictability of presence of TE excisions based on SNP frequencies in O. sativa and O. glaberrima.

	Polymorphism	Test set	Control set	P value
	TE excisions	16	2	0.0003
	TE insertions	16	27	0.026
	Repeat slippage	20	17	0.53
	InDel	25	15	0.18

NATURE COMMUNICATIONS | 7:12790 | DOI: 10.1038/ncomms12790 | www.nature.com/naturecommunications

Acknowledgements

We thank Detlef Weigel and Claude Becker for providing methylome data for O. sativa and O. glaberrima. This material is based on work supported by the Swiss National Foundation grant #31003A_138505/1 to T.W., by the US National Science Foundation under grants #0321678, #0638541, #0822284 and #1026200 to Y.Y., and R.A.W. and the Bud Antle

Author contributions

T.W. designed the study, analysed the data and wrote the paper. S.R. helped design the study, created software, analysed the data and provided critical input in writing the manuscript. R.A.W., Y.Y., S.R., M.C., P.R.M. and A.Z. produced the genome sequence. K.F.X.M., G.H. and O.P. produced the genome annotation.

Additional information

Supplementary Information accompanies this paper at http://www.nature.com/ naturecommunications Competing financial interests: The authors declare no competing financial interests.

Reprints and permission information is available online at http://npg.nature.com/ reprintsandpermissions/


Dr C Fauvelot
email: cecile.fauvelot@unibo.it

C Lemaire

S Planes

F Bonhomme

Inferring gene flow in coral reef fishes from different molecular markers: which loci to trust?

Keywords: coral reef fishes, allozymes, introns, mtDNA, gene flow, neutrality

Contrasting results are usually reported in the literature regarding the factors influencing observed structuring of genetic variability. The goals of this study were, for five coral reef fishes in French Polynesia, (1) to infer the theoretical variance of single locus FST estimates expected under neutrality in order to exclude outlier loci before inferring gene flow and (2) to test thereafter whether species laying pelagic eggs effectively disperse more than species laying benthic eggs in this system. For this purpose, a total of 952 individuals from five species belonging to two families (Chaetodontidae and Pomacentridae) were screened among populations sampled within a 60-600km spatial range for intron length polymorphism at 11 loci in order to illuminate contrasting results previously published on allozymes and mitochondrial DNA (mtDNA) control region polymorphisms. Statistically speaking, among the five species, four loci (three allozymes and one intron) were identified as outliers and discarded before interpretation of genetic differentiation in terms of effective dispersal. Biologically speaking, our results suggest that the observed genetic structure is not significantly related to the reproductive strategy of coral reef fish in the island system we analysed and that observed random genetic differentiation accommodates Wright's island model in all five species surveyed. Overall, our study emphasizes how cautious one has to be when trying to interpret present-day genetic structure in terms of gene flow while using a limited number of loci and/or different sets of loci.

INTRODUCTION

Studying spatial partitioning of genetic variation is viewed as a robust method for inferring gene flow among discrete populations, conditional on taking into account all processes shaping the present genetic structure [START_REF] Hewitt | The genetic legacy of the Quaternary ice ages[END_REF]. For coral reef fishes, contrasting results have been found relative to the relationship between present genetic structure and larval duration [START_REF] Waples | A multispecies approach to the analysis of gene flow in marine shore fishes[END_REF][START_REF] Doherty | Gene flow and larval duration in seven species of fish from the Great Barrier Reef[END_REF][START_REF] Shulman | Early life histories, ocean currents, and the population genetics of Caribbean reef fishes[END_REF][START_REF] Rocha | Adult habitat preferences, larval dispersal, and the comparative phylogeography of three Atlantic surgeonfishes (Teleostei: Acanthuridae)[END_REF][START_REF] Bay | The relationship between population genetic structure and pelagic larval duration in coral reef fishes on the Great Barrier Reef[END_REF], reproductive strategy [START_REF] Shulman | Early life histories, ocean currents, and the population genetics of Caribbean reef fishes[END_REF][START_REF] Planes | Evidence of limited gene flow in three coral reef fish species in the lagoon of New Caledonia[END_REF], species swimming abilities [START_REF] Waples | A multispecies approach to the analysis of gene flow in marine shore fishes[END_REF][START_REF] Riginos | Larval spatial distributions and other early life-history characteristics predict genetic differentiation in eastern Pacific blennioid fishes[END_REF][START_REF] Purcell | Weak genetic structure indicates strong dispersal limits: a tale of two coral reef fish[END_REF] or geographic distances [START_REF] Planes | Genetic structure of Dascyllus aruanus populations in French Polynesia[END_REF][START_REF] Planes | A genetic metapopulation model for reef fishes in oceanic islands: the case of the surgeonfish, Acanthurus triostegus[END_REF][START_REF] Planes | Isolation by distance and vicariance drive genetic structure of a coral reef fish in the Pacific Ocean[END_REF][START_REF] Purcell | Weak genetic structure indicates strong dispersal limits: a tale of two coral reef fish[END_REF]. [START_REF] Fauvelot | Understanding origins of presentday genetic structure in marine fish: biologically or historically driven patterns?[END_REF] have shown, using allozymes, that the observed significant structuring of genetic variation among populations of eight coral reef fishes in a fragmented island system was not correlated with geographic distances at the studied spatial scale (ca 1500 km) nor associated with species reproductive strategy and larval duration. On average, only a third of the polymorphic loci contributed to the overall genetic differentiation. Conversely, the analysis of the genetic variation at the mitochondrial DNA (mtDNA) control region sequences in a subset of individuals revealed genetically homogeneous populations in all seven species under study [START_REF] Fauvelot | Reductions in the mitochondrial DNA diversity of coral reef fish provide evidence of population bottlenecks resulting from Holocene sea-level change[END_REF].

Such puzzling discordant patterns of genetic variation obtained at different loci are increasingly described [START_REF] Karl | Balancing selection at allozyme loci in oysters: implications from nuclear RFLPs[END_REF][START_REF] Pogson | Genetic population structure and gene flow in the Atlantic cod Gadus morhua: a comparison of allozyme and nuclear RFLP loci[END_REF][START_REF] Lemaire | Do discrepancies between microsatellite and variation reveal differential selection between sea and lagoon in the sea bass (Dicentrarchus labrax)?[END_REF][START_REF] De Innocentiis | Allozyme and microsatellite loci provide discordant estimates of population differentiation in the endangered dusky grouper (Epinephelus marginatus) within the Mediterranean Sea[END_REF][START_REF] Dufresne | Differential patterns of spatial divergence in microsatellites and allozymes alleles: further evidence for locus-specific selection in the acorn barnacle, Semibalanus balanoides?[END_REF][START_REF] Arnaud-Haond | Large discrepancies in differentiation of allozymes, nuclear and mitochondrial DNA loci in recently founded Pacific populations of the pearl oyster Pinctada margaritifera[END_REF][START_REF] Gomez-Uchida | Allozyme and AFLP analysis of genetic population structure in the hairy edible crab Cancer setosus from the Chilean coast[END_REF][START_REF] Nielsen | Evidence of microsatellite hitch-hiking selection in Atlantic cod (Gadus morhua L): implications for inferring population structure in nonmodel organisms[END_REF] and lead to question the interpretation made of genetic differentiation (FST) estimates in terms of effective dispersal. A key issue is that, under neutrality, all loci should theoretically be affected similarly by drift and migration [START_REF] Lewontin | Distribution of gene frequency as a test of the theory of selective neutrality of polymorphisms[END_REF][START_REF] Slatkin | Gene flow and the geographic structure of natural populations[END_REF], so that FST estimates should be homogeneous among loci, and multilocus FST estimates may be used to infer gene flow among populations. Yet, erroneous gene flow estimates may result from the interpretation of allelic distribution over multiple loci for two main reasons. One lies in the underlying high variance of the migration/drift process itself and its sensitivity to departure from equilibrium [START_REF] Whitlock | Indirect measures of gene flow and migration: FST does not equal 1/(4Nm+1)[END_REF]. The other resides in selection or other non-neutral mechanisms (acting directly or through hitch-hiking) differentially affecting loci [START_REF] Ll | Population structure and human evolution[END_REF][START_REF] Slatkin | Gene flow and the geographic structure of natural populations[END_REF]. Indeed, disruptive selection will result in higher FST estimates than expected under neutrality, inflating the overall genetic differentiation. On the other hand, loci experiencing balancing selection will have allele frequencies more similar than expected under neutrality, reducing the FST estimates and further giving the false view of high gene flow.

In order to properly infer gene flow among populations, it appears therefore essential to (1) estimate genetic differentiation based on reasonable numbers of individuals (50-100 individuals; Ruzzante, 1998), loci (15-20 loci;[START_REF] Nielsen | Evidence of microsatellite hitch-hiking selection in Atlantic cod (Gadus morhua L): implications for inferring population structure in nonmodel organisms[END_REF] and alleles (either few loci with many alleles or many loci with a few alleles; [START_REF] Kalinowski | How many alleles per locus should be used to estimate genetic distances?[END_REF] and (2) discard loci behaving non-neutrally before the estimation of genetic differentiation. Several neutrality tests are available in the current literature that may allow detection of selected loci at ecological timescales (reviews in [START_REF] Guinand | How to detect polymorphisms undergoing selection in marine fishes? A review of methods and case studies, including flatfishes[END_REF][START_REF] Beaumont | Adaptation and speciation: what can Fst tell us?[END_REF]. Among these tests, those derived from the Lewontin-Krakauer method (that is, based on multilocus comparisons and the identification of FST outliers) appear robust when allele distributions are available at multiple loci [START_REF] Beaumont | Evaluating loci for use in the genetic analysis of population structure[END_REF][START_REF] Beaumont | Identifying adaptive genetic divergence among populations from genome scans[END_REF][START_REF] Beaumont | Adaptation and speciation: what can Fst tell us?[END_REF].

The goal of the present study is, for five species of coral reef fish, to infer gene flow among insular populations sampled within a 60-600km spatial range and to test for the influence of reproductive strategies. Considering the problem of neutrality among loci, the theoretical neutral variance of single locus FST estimates in our insular system was first investigated in order to identify (and exclude from calculation) outlier loci before computing gene flow. For this purpose, the number of loci was increased to improve the power of the test. We reanalysed all previously assayed individuals from five species with a new set of loci with potentially higher mutation rates and gene diversities than allozymes, that is, intron length polymorphisms. Unlike exons, the nucleotide sequences of introns are little constrained so that introns accumulate mutations rapidly. Therefore, introns often display a high genetic variability [START_REF] Palumbi | Nucleic acids II: the polymerase chain reaction[END_REF], which make introns being targeted as markers of population variation and subdivision [START_REF] Lessa | Rapid surveying of DNA sequence variation in natural populations[END_REF][START_REF] Borsa | Genetic structure of round scad mackerel Decapterus macrosoma (Carangidae) in the Indo-Malay archipelago[END_REF][START_REF] Hassan | Lessepsian invasion without bottleneck: example of two rabbitfish species (Siganus rivulatus and Siganus luridus)[END_REF][START_REF] Friesen | Population genetic structure and conservation of marbled murrelets (Brachyramphus marmoratus)[END_REF][START_REF] Rohfritsch | Genetic structure of Indian scad mackerel Decapterus russelli: Pleistocene vicariance and secondary contact in the Central Indo-West Pacific Seas[END_REF][START_REF] Berrebi | Population structure and systematics of Opsariichthys bidens (Osteichthyes: Cyprinidae) in south-east China using a new nuclear marker: the introns (EPIC-PCR)[END_REF]. Exon-primed, intron-crossing polymerase chain reaction (EPIC-PCR) amplification is an appreciable strategy for finding DNA polymorphism in several species of different families [START_REF] Hassan | Seventeen new exon-primed introncrossing polymerase chain reaction amplifiable introns in fish[END_REF] and an alternative to the development of microsatellite loci. Length polymorphism at 11 intron loci was assessed in five coral reef fish species and the corresponding new data set was further added to the previous allozymes data set [START_REF] Fauvelot | Understanding origins of presentday genetic structure in marine fish: biologically or historically driven patterns?[END_REF] to conduct global analyses. The comparison of the distribution of genetic variation obtained from different markers using the same set of individuals from related species collected in the same area also provide a solid framework to test among mutation, migration and genetic drift as the main forces driving the distribution of genetic polymorphisms.

MATERIALS AND METHODS Sampling

Sampling was carried out as described in [START_REF] Fauvelot | Understanding origins of presentday genetic structure in marine fish: biologically or historically driven patterns?[END_REF]. Briefly, five different islands were sampled across French Polynesia based on a standard sampling effort of 50 individuals per island per species. Because the eight species selected for the original study were not all present in the five islands, and in order to compare effective dispersal among species at similar scales, we restricted further analysis to five species inhabiting four common islands: Rangiroa, Takapoto, Tetiaroa and Moorea (Figure 1). Geographic distances between islands ranged from 60km (between Moorea and Tetiaroa) to 566km (between Moorea and Takapoto). Chaetodon citrinellus is a butterflyfish strictly inhabiting lagoons while C. quadrimaculatus and Forcipiger flavissimus are both outer reef butterflyfishes. Dascyllus aruanus and Pomacentrus pavo are two damselfishes strictly inhabiting lagoons. Fishes were collected, tissue stored and DNA extracted as described in [START_REF] Fauvelot | Reductions in the mitochondrial DNA diversity of coral reef fish provide evidence of population bottlenecks resulting from Holocene sea-level change[END_REF]. A total of 952 individuals from five species (Table 1) were used for further genetic screening.

Genetic analysis -EPIC-PCR

Eleven introns loci were tested for length polymorphism for each species: intron 7 of creatine kinase; intron 2 of glyceraldehyde-3-phosphate dehydrogenase; introns 2 and 5 of growth hormone; introns 1, 2, 4 and 5 of aldolase B; introns 1, 2 and 3 of gonadotropin-releasing hormone [START_REF] Hassan | Seventeen new exon-primed introncrossing polymerase chain reaction amplifiable introns in fish[END_REF]. PCR amplification and electrophoresis follow protocols described in [START_REF] Hassan | Seventeen new exon-primed introncrossing polymerase chain reaction amplifiable introns in fish[END_REF]. Briefly, amplified PCRs (with forward primer labelled radioactively with 33P) were loaded on a 6% polyacrylamide gel in TBE buffer. The gels were run for 3 h at 50W. After drying, the gels were exposed against films of Biomax autoradiograph. Individual genotypes were directly scored from the autoradiograph films with alleles being identified according to their size relative to the most common allele, as for allozymes. To ensure scoring accuracy and consistency across gels and populations, a locusspecific size standard from a pool of individuals with known genotypes was run at both sides of each gel. With this method, we were able to segregate 1 bp length difference bands. Individuals for which no PCR products were visualized on gel were further re-amplified in a new PCR using a 2°C lower annealing temperature. Results of the amplification successes across species are presented in [START_REF] Hassan | Seventeen new exon-primed introncrossing polymerase chain reaction amplifiable introns in fish[END_REF] and detailed PCR conditions for each species and locus are available from the authors.

Data analysis

The allozymes data set was reanalysed for genetic divergences in the overall populations following the removal of populations from the data set (see sampling). Allele frequencies, number of alleles per locus (NA), observed heterozygosity (HO) and unbiased gene diversity (HE; [START_REF] Nei | Molecular Evolutionary Genetics[END_REF] were calculated within each population and over all populations for both allozymes and introns. Linkage disequilibrium between loci was tested using the procedure of [START_REF] Black | A FORTRAN program for the calculation and analysis of two-locus linkage disequilibrium coefficients[END_REF]. Deviation from Hardy-Weinberg expectations (HWE) at introns loci was measured using [START_REF] Weir | Estimating F-statistics for the analysis of population structure[END_REF] inbreeding coefficient f (equivalent to Wright's FIS), and locus conformance to HWE was assessed by permuting alleles within populations (1000 permutations). The presence of null alleles was tested using the program MICRO-CHECKER (Van Oo sterhout et al., 2004). Genetic divergences were estimated using the [START_REF] Weir | Estimating F-statistics for the analysis of population structure[END_REF] estimator. For each species, overall single locus FST values were estimated by jackknifing over populations and overall global FST estimates were obtained by jackknifing over loci. Unless specified, all computations were conducted using GENETIX [START_REF] Belkhir | Genetix 404 A windows program for population genetics[END_REF]. The theoretical distribution of single locus FST estimates expected under neutrality was computed through simulations and outlier loci identified using [START_REF] Beaumont | Evaluating loci for use in the genetic analysis of population structure[END_REF] test. The principle of this test is to estimate, conditional on gene diversity, the range of FST values expected under the combined effect of migration and genetic drift in a model approaching the studied system, and to compare the observed FST at this given locus to the simulated empirical distribution. The models for computer simulations were chosen as to fit at maximum the conditions of the studied system. Indeed, we chose the infinite allele model of evolution [START_REF] Kimura | The number of alleles that can be maintained in a finite population[END_REF], and the infinite island model of [START_REF] Wright | Isolation by distance[END_REF] as this is the best-suited available model for fragmented insular environments such as Pacific reefs. It should be noted that, for heterozygosities lower than 0.8 (which is our case), stepwise and infinite-allele models lead to similar simulated distributions of FST [START_REF] Beaumont | Evaluating loci for use in the genetic analysis of population structure[END_REF], so that the choice of the mutation model has little impact on the simulated empirical distribution. Coalescent simulations were performed to generate data sets with a distribution of FST close to the empirical distribution. Simulations were conducted assuming four sampled populations, with a metapopulation size of 100 populations (the approximate number of islands in French Polynesia), sample sizes of 50 individuals and an expected FST equal to the median estimate of FST calculated for all loci. Under the infinite allele mutation model, 20 000 HE and FST paired values were generated and the P-values for each locus were then calculated using the procedure described in [START_REF] Beaumont | Evaluating loci for use in the genetic analysis of population structure[END_REF]. Loci with unusually low or high FST estimates (that is, out of the 95% distribution of simulated FST), conditional on heterozygosity, were regarded as outlier loci.We used the currently distributed version of the test, namely Fdist2 (http://www.rubic.reading.ac.uk/~mab/software/fdist2.zip), described in [START_REF] Beaumont | Identifying adaptive genetic divergence among populations from genome scans[END_REF]. Because the test conducts simulations based on initial expected FST value set up by the user, this value was first fixed equal to the global FST (estimated over all loci). However, because this median value may be biased owing to the presence of an outlier locus, we systematically conducted again the simulations with a new global FST estimated omitting the outlier locus (or loci), and so on until no more loci were identified as outliers. Overall populations and pairwise FST were estimated for each species over all nuclear loci (allozymes and introns) once removing outliers detected by the [START_REF] Beaumont | Evaluating loci for use in the genetic analysis of population structure[END_REF] test. Significance of pairwise genetic divergence was tested using a permutation procedure (1000 permutations) in GENETIX [START_REF] Belkhir | Genetix 404 A windows program for population genetics[END_REF]. Mantel tests were performed to test for a correlation between geographical distance and multilocus FST (approach of [START_REF] Hutchison | Correlation of pairwise genetic and geographic distance measures: inferring the relative influences of gene flow and drift on the distribution of genetic variability[END_REF]. For each species, Mantel coefficient Z was computed and tested for significance through a permutation procedure with 100 permutations, under the null hypothesis that both matrices are not linearly related. To evaluate the level of correlation between genetic patterns exhibited by the different species, matrices of pairwise FST were compared with a Mantel test: Pearson correlation coefficient r was computed for each pair of species and its significance tested using a permutation procedure with 100 permutations. All tests were performed in GENETIX using the Mantel option.

RESULTS

Out of the 11 intron primer pairs tested, five loci (AldoB1, AldoB5, CK7, GnRH3-3 and Gpd2) could not be further used for the five species because they lead to multibanded (several loci amplified), monomorphic or not clearly interpretable patterns (Table 2). The remaining six loci (AldoB2, AldoB4, GnRH3-1, GnRH3-2, GH2 and GH2) provided single locus amplifications with clearly interpretable length polymorphism, and finally, 2-3 polymorphic intron per species could be analysed (Table 2). Over all populations, the number of alleles per locus ranged from 2 to 12 (for GH2 over all F. flavissimus sampled populations). Gene diversity among populations in each species was more or less even and variations in gene diversity were locus specific, such as that observed at enzymatic loci [START_REF] Fauvelot | Understanding origins of presentday genetic structure in marine fish: biologically or historically driven patterns?[END_REF], with no congruent patterns among species. Over all four populations, gene diversity at intron loci reached 0.726 for F. flavissimus at GH2, and 0.501 at enzymatic loci for F. flavissimus at Est1 (Table 3). After standard Bonferroni correction of the 5% significance level, a single locus significantly deviated from HWE; GnRH3-1 in C. citrinellus populations. The cause of such departures is likely to be the presence of null alleles because the deviations resulted from significant heterozygote deficiencies, observed in all populations in this locus only. The occurrence of null alleles was further confirmed using MICRO-CHECKER and the average null allele frequency at GnRH3-1 locus across populations was 0.165 (estimator of [START_REF] Brookfield | A simple new method for estimating null allele frequency from heterozygote deficiency[END_REF]. However, this locus was maintained for further analysis, and genotypes were not adjusted according to nulls to allow for multilocus analysis. Significant (P<0.01) linkage disequilibria were found between introns originating from the same gene (that is, GnRH3-1 and GnRH3-2 for C. citrinellus and C. quadrimaculatus; between AldoB2 and AldoB4 for D. aruanus). For C. citrinellus, the association was significant between the two most common alleles at each locus, and between rare alleles at GnRH3-1 with two of the four alleles at GnRH3-2. For C. quadrimaculatus and D. aruanus, the linkage involved rare alleles from one locus significantly associated with common alleles from the other locus. However, all loci were maintained for further analysis. For each species, over all four populations, a wide range of single locus FST was observed among nuclear loci (allozymes and introns), from negative values up to 0.198 at ADA-1 for D. aruanus. Multilocus FST estimates ranged from 0.008 for C. citrinellus to 0.026 for D. aruanus (Table 3). Discrepancies among classes of markers were observed because only at enzymatic loci, global FST estimates ranged from 0.012 for C. citrinellus to 0.031 for F. flavissimus and D. aruanus, while when looking at only introns, global FST estimates were zero for all species except one, F. flavissimus (FST = 0.008, P>0.05). However, these discrepancies were not marker related but locus related because single locus FST estimates at introns were similar to those of most of allozymes within species (Table 3). [START_REF] Beaumont | Evaluating loci for use in the genetic analysis of population structure[END_REF] test revealed that three loci demonstrated significantly higher FST values than expected under the null hypothesis of neutrality: ADA-1 in D. aruanus (P = 0.991), and ADA and GDA in F. flavissimus (P>0.999 and P=0.953, respectively) (Figure 2). One additional locus, AldoB4 in P. pavo, showed significantly lower FST values (P=0.037) than expected under neutrality. Recalculating the global FST for each species without the outliers loci did not induce further rejection of neutrality for any other locus. Excluding outlier loci detected by [START_REF] Beaumont | Evaluating loci for use in the genetic analysis of population structure[END_REF] test, global FST were significant for all five species and ranged from 0.008 (P=0.01) for C. citrinellus to 0.019 (P=0.02) for D. aruanus (Table 4). Pairwise FST estimated from all nuclear loci excluding outliers ranged from 0 to 0.045 (between Rangiroa and Moorea for P. pavo; Table 4). None of the five species expressed significant correlation between genetic divergences and geographic distances between populations (all P>0.05). The comparison of pairwise FST matrices among the five species revealed that two species showed a significant correlation of their patterns of genetic structuring throughout the four sampled islands: C. quadrimaculatus and D. aruanus (r=0.843, P=0.03; see Table 4).

DISCUSSION FIGURE 1

Location of sampling sites of the five fish species in French Polynesia (South Pacific Ocean). Details of sites are reported in [START_REF] Fauvelot | Understanding origins of presentday genetic structure in marine fish: biologically or historically driven patterns?[END_REF] and sampling effort is indicated in Table 1.

FIGURE 2

Summary of the results of the [START_REF] Beaumont | Evaluating loci for use in the genetic analysis of population structure[END_REF] neutrality test. For each locus, estimate of FST is plotted against its corresponding transformed P-value (expressed as logit(2|P-0.5|) as in [START_REF] Beaumont | Identifying adaptive genetic divergence among populations from genome scans[END_REF], corresponding to the rejection of the null hypothesis of neutrality. The vertical dotted line bar shows the critical transformed P-values used for identifying outlier loci at the 5% confidence level. Each outlier locus is denoted as follows: the two first letters refer to the species (Da, Dascyllus aruanus; Ff, Forcipiger flavissimus; Pp, Pomacentrus pavo), followed by the locus name. [START_REF] Nei | Molecular Evolutionary Genetics[END_REF]; HO, observed heterozygosity; M, monomorphic (no length polymorphism); Mo, Moorea; Multi, multiband pattern (more than one locus amplified); NA, number of alleles; nt, not tested; Poly, length polymorphism but hardly interpretable; Ra, Rangiroa; Ta, Takapoto; Te, Tetiaroa. Intron names are given as in [START_REF] Hassan | Seventeen new exon-primed introncrossing polymerase chain reaction amplifiable introns in fish[END_REF]. Significant deviation from Hardy-Weinberg proportions is denoted by * for P<0.05 and ** for P<0.01.

TABLE 3

Gene diversity (HE) over all four sampled populations and FST estimates for each scored locus (allozymes and introns) per species

TABLE 4

Pairwise genetic differentiation (FST) estimated among islands excluding outlier loci Abbreviations: Mo, Moorea; Ra, Rangiroa; Ta, Takapoto; Te, Tetiaroa. Significant values after sequential Bonferroni correction are indicated in bold.

TABLE 1

 1 Number of individuals scored at enzymatic and introns loci, indicated per species and per site

TABLE 2

 2 Amplifications results of the 11 tested

	Abbreviations:
	HE, unbiased gene diversity