{'fraction_non_alphanumeric': 0.030822592723976194, 'fraction_numerical': 0.05281620748098222, 'mean_word_length': 4.7105559306774705, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Chromosome-level assemblies of multiple Arabidopsis 1 thaliana accessions reveal hotspots of genomic 2 rearrangements 3 4

We report chromosome-level, reference-quality assemblies of seven Arabidopsis thaliana accessions selected across the geographic range of this model plant. Each genome assembly revealed between 13-17 Mb of rearranged and 5-6 Mb of novel sequence introducing copy-number changes in ~5000 genes, including ~1,900 genes which are not part of the current reference annotation. By analyzing the collinearity between the genomes, we identified ~350 hotspots of rearrangements covering ~4% of the euchromatic genome. Hotspots of rearrangements are characterized by accession-specific accumulation of tandem duplications and are enriched for genes implicated in disease resistance and secondary metabolite biosynthesis. Loss of meiotic recombination in hybrids within these regions is consistent with the accumulation of rare and deleterious alleles and incompatibility loci. Together this suggests that hotspots of rearrangements are governed by different evolutionary dynamics as compared to the rest of the genome and facilitate rapid responses to the ever-evolving challenges of biotic stress.


rare and deleterious alleles and incompatibility loci. Together this suggests that 23
hotspots of rearrangements are governed by different evolutionary dynamics as 24 compared to the rest of the genome and facilitate rapid responses to the ever-evolving 25 challenges of biotic stress. 26 The first complete assembly of a plant genome, the reference sequence of A. thaliana 27 (Col-0), was based on a minimal tiling path of BACs sequenced with Sanger technology and 28 was released in the year 2000 1 . Since then, this reference has been widely used for 29 identification of genetic variants, which typically relied on short-read based resequencing or 30 reference-guided assembly 2-8 . Although millions of small variants have been identified, the 31 identification of large genomic rearrangements remained challenging. In contrast, reference-32 independent, chromosome-level assemblies promise accurate identification of all sequence 33 differences independent of their complexity 9 . So-far, however, there are only few reports on 34 2 whole-genome de-novo assemblies for A. thaliana available and the assemblies have not 35 been thoroughly compared to each other 10-13 . 36 Using deep PacBio (45-71x) and Illumina (56-78x) whole-genome shotgun sequencing, 37 we assembled the genomes of seven accessions from geographically diverse populations 38 including An-1 (Antwerpen, Belgium), C24 (Coimbra, Portugal), Cvi-0 (Cape Verde Islands), 39 Eri-1 (Eringsboda, Sweden), Kyo (Kyoto, Japan), Ler (Gorzów Wielkopolski, Poland) and Sha 40 (Shahdara, Tadjikistan) (Supplementary Table 1). The assembly of Ler was already described 41 in a recent study and used for the development of a whole-genome comparison tool 9 , however, 42 as it was generated in the same process it is integrated in this study as well. These accessions 43 (together with the reference accession Col-0) were initially selected as the founder lines of 44 Arabidopsis Multi-parent Recombination Inbreeding Lines (AMPRIL) 14 population. The contig 45 assemblies featured N50 values from 4.8 -11.2 Mb and chromosome-normalized L50 46 (CL50) 15 values of 1 or 2 indicating that nearly all chromosome arms were assembled into a 47 few (typically one to five) contigs only (Table 1 and Supplementary Table 2). We generated 48 chromosome-level assemblies based on the homology of these contigs to the reference 49 sequence and validated the resultant contig order of two of the assemblies with three different 50 genetic maps, where we did not find evidence for a single mis-placed contig (Supplementary 51 Table 3). The seven chromosome-level assemblies consisted of 43 -73 contigs and reached 52 a total length of 117.7 -118.8 Mb, which is very similar to the 119.1 Mb of the reference 53 sequence ( Fig. 1 and Table 1) and even included parts of the highly complex regions of 54 centromeres, telomeres and rDNA clusters (Supplementary Table 4 and 5). The remaining 55 unanchored contigs had a total length of 1.5 -3.3 Mb and consisted almost entirely of repeats, 56 which agrees with gaps between the contigs, which most of them (84.5% -98.0%) were 57 introduced due to repetitive regions (Supplementary Table 6). 58 Between 99.1% and 99.4% of the reference genes 16 could be aligned to each of the 59 seven assemblies. Almost all of the remaining genes were truly missing in the genomes as 60 we could confirm by mapping their Illumina short reads against the reference sequence 61 (Supplementary Table 7), suggesting that the assemblies covered almost all of the genic 62 regions. In agreement with this, we annotated 27,098 to 27,574 protein-coding genes in each 63 of the assemblies, which is similar to the 27,445 genes annotated in the reference sequence 16 64 ( Table1, Supplementary Table 8 and 9). 65 The lack of contiguity of short sequencing reads makes resequencing-based analyses 66 mostly blind for large rearrangements. In contrast, the high contiguity of these new assemblies 67 now enables the comprehensive description of complex structural rearrangements including 68 inversions, translocations and duplications (Fig. 2a). By comparing each of the new 69 assemblies against the reference sequence using the whole-genome comparison tool SyRI 9 , 70 we found between 102.2 and 106.6 Mb of syntenic regions and between 12.6 and 17.0 Mb of 71 3 rearranged regions in each of the genomes. The rearrangements included 1.5 -4.2 Mb (33 -72 46) inversions, 1.0 -1.7 Mb (364 -566) intra-chromosome translocations, and 0.9 -1.3 Mb 73 (365 -626) inter-chromosome translocations. Apart from balanced structural variation, we 74 also found 7. 150) of duplication loss and gain variation ( Fig. 2b and 75 Supplementary Table 10). Similar to sequence variation, rearrangements were not evenly 76 distributed along the chromosomes, but were enriched in pericentromeres (Supplementary  77   Table 11). Their lengths ranged from a few dozen bp to hundreds of kb and even Mb scale 78 (Fig. 2c). Among the largest differences were inversions including a 2.48 Mb inversion on 79 chromosome 3 of Sha (Supplementary Fig. 1  inversions and translocations were generally lower as compared to duplications, where more 86 than 50% of them were shared by at least two accessions suggesting differences in the 87 selection pressures acting on balanced and non-balanced variation (Fig. 2e). 88 In each of the pairwise comparisons, we identified 5.1 -6.5 Mb accession-specific 89 sequence (Fig. 2b). Using these regions and their overlaps for a pan-genome analysis 19,20 , 90 we estimated a pan-genome size of ~135 Mb and a core-genome size of ~105 Mb (Fig. 2f). 91 A similar analysis with genes modelled a gene pan-genome size of ~30,000 genes illustrating 92 that one reference genome is not sufficient to capture the entire sequence diversity within A. 93 thaliana ( Supplementary Fig. 3). 94 Genomic rearrangements have the potential to delete, create or duplicate genes resulting 95 in gene copy number variation (CNV). Based on gene family clustering of all genes in all eight 96 accessions 21 we found 22,040 gene families with conserved copy number, while 4,957 gene 97 families showed differences in gene copy numbers ( Fig. 2g and Supplementary Table 14). 98 Almost 99% of the copy-variable gene families had 5 or less copies, while only less than 10% 99 of them showed more than two different copy numbers across the eight accessions 100 ( Supplementary Fig. 4). Among the copy-variable genes we found 1,941 non-reference gene 101 families including 891 gene families present in at least two of the other accessions (Fig. 2h). 102 Around 23% of the non-reference gene families featured orthologs in the closely related 103 genome of Arabidopsis lyrata and, according to RNA-seq read mapping, 26%-40% of them 104 showed evidence of expression in the individual accessions (Supplementary Table 15). The 105 remaining 1,050 non-reference gene families, which were evenly distributed across the 106 accessions ( Fig. 2g), with the only exception of Cvi-0, where we found nearly twice as many 107 (214) accession-specific genes which is in agreement with the divergent ancestry of this relict 108 4 accession 4,22 . 109 The high contiguity of the assembled sequences enabled the first analysis of the 110 conservation of collinearity between different Arabidopsis genomes. For this we introduced a 111 new concept called Synteny Diversity p syn , which is similar to Nucleotide Diversity 23 , however, 112 instead of measuring average sequence differences it measures average fraction of non-113 collinear genome pairs in a given population. p syn values can range from 0 to 1, where 1 refers 114 to the complete absence of collinearity between any of the genomes and 0 to regions where 115 all genomes are collinear. p syn can be calculated in any given region, however, the annotation 116 of synteny still needs to established within the context of the whole genomes to avoid false 117 assignments of homologous but non-allelic sequence (Methods). 118 We calculated p syn in 5-kb sliding windows across the genome using pair-wise 119 comparisons of all eight accessions (Fig. 3a). As expected, p syn was generally high in 120 pericentromeric regions and low in chromosome arms. Overall, this revealed around 90 Mb 121 (76% of the genome) where all genomes were syntenic to each other, while for the remaining 122 29 Mb (24%) the collinearity between the genomes was not conserved. This, for example, 123 included a region on chromosome 3 (ranging from Mb ~2.8 -5.3), where p syn was increased 124 to ~0.25 due to the 2.48 Mb inversion specific to Sha (Fig. 3a, arrow labelled with (a)). 125 Unexpectedly, however, some regions featured p syn values even larger than 0.5 indicating 126 that not only two, but multiple independent, non-syntenic haplotypes segregate implying that 127 these regions are more likely to undergo or conserve complex mutations as compared to the 128 rest of the genomes. Overall, we found 576 such hotspots of rearrangements (HR) with a total 129 size of 10.2 Mb including 351 regions in the gene-rich euchromatic chromosome arms with a 130 total length of 4.1 Mb (or 4% of euchromatic genome) (Supplementary Table 16). 131 Even though HRs in euchromatic regions included more transposable elements and less 132 genes as compared to conserved regions, they still contained significant numbers of genes, 133 many of which occurred at high and variable copy-number between the accessions (Fig. 3b, 134 c). For example, a HR on chromosome 4, which overlapped with the RPP4/RPP5 R gene 135 cluster 24,25 , displayed five to 15 intact or truncated copies of the RPP5 gene within the eight 136 genomes ( Fig. 3d and Supplementary Table 17). The different gene copies were primarily 137 introduced by an accumulation of forward tandem duplications and large indels (Fig. 3e), a 138 pattern, which was shared by many HR regions ( Supplementary Fig. 5). While the individual 139 duplications were not conserved between the haplotypes, the borders of HRs were typically 140 well conserved (Fig. 3f). This suggested that either different selection regimes introduced 141 clear-cut borders between highly HRs and their surrounding regions, or that increased tandem 142 duplication rates were exactly limited to the HR regions. Such a local increase of mutation 143 rates could potentially be enabled by non-allelic homologous recombination (NAHR) which 144 5 could be triggered by the high number of local repeats in these regions [26][27][28] . In any case the 145 accession-specific accumulation of duplications suggested that the HR regions are only 146 partially linked to their genomic vicinity. To test this, we analysed the linkage disequilibrium 147 (LD) within 1,135 genomes of the 1001 Genomes Project 4 across the HR regions. LD was 148 high in the regions around the HR and was also high within HRs, however, when calculated 149 across the HR border LD was significantly lower supporting the idea that HRs are not linked 150 to the surrounding haplotypes but that they are segregating as large non-recombining units 151 ( Fig. 4f). To test if HRs are indeed depleted for meiotic recombination, we overlapped them 152 with 15,683 crossover (CO) sites previously identified within Col-0/Ler F2 progenies 29,30 . Only 153 64 of them partially overlapped with non-syntenic regions while all other COs co-located 154 syntenic regions (Fig. 4b), suggesting that HRs are almost completely silenced for COs (P-155 value < 2e-16) and that they follow different evolutionary dynamics as compared to the rest 156 of the (recombining) genome. 157 While lack of meiotic crossovers coupled with increased duplication rates can lead to 158 rapid generation of new haplotypes, reduced meiotic recombination also has been linked to 159 the accumulation of new (deleterious) mutations 31 . In agreement with this HRs also showed 160 an accumulation of single nucleotide polymorphisms with low allele frequencies and 161 potentially deleterious variation as compared to other regions in the genome (Fig. 4c, d and 162 Supplementary Fig. 6). Moreover, reduced recombination combined with geographic isolation 163 can provide the basis for the development of alleles, which are incompatible with distantly 164 related haplotypes leading to intra-species incompatibilities 32 . To test this, we searched the 165 location of nine recently reported genetic incompatible loci 33 (DM1-9) and found that all except 166 of one overlapped with HRs, while DM3, the locus which did not overlap with HRs, was closely 167 flanked by two HRs . In addition, we also checked the 168 locus of a recently published single-locus genetic incompatibility 34 and found that it was also 169 residing in a HR region ( Supplementary Fig. 13). 170 The excess of tandem duplications, accumulation of accession-specific gene-copy 171 numbers and increased sequence diversity paired with the lack of meiotic recombination and 172 the accumulation of deleterious alleles were reminiscent of the patterns that have been 173 described for R gene clusters [35][36][37][38][39][40][41][42] . In fact, the 808 reference genes in HRs were significantly 174 enriched for genes involved in defense response, signal transduction and secondary 175 metabolite biosynthesis (Fig. 4e) suggesting a reoccurring role of HRs in the adaptation to 176 biotic stress. As biotic challenges are constantly changing as pathogens evolve, it has been 177 proposed that the accumulation of new gene duplicates could increase the genomic diversity 178 and thereby enable rapid genomic changes supporting the genomic response of plants 179 against pathogens 28,43-45 . 180 The availability of these high-quality, chromosome-level genomes assemblies now 181 6 enables the simultaneous access to all these regions, which can now be studied not only 182 locally, but in the context of entire genomes. For example, all but one of the 47 NB-LRR R-183 gene clusters were completely reconstructed in all of the seven new assemblies (Fig. 5,184 Supplementary all rearrangements, what could be essential as many of the HRs were only covered with a 186 few or even no allelic markers and could easily be missed with genome-wide selection or 187 diversity-scans ( Supplementary Fig. 14 and 15). The genome-wide estimation of Synteny 188 Diversity might be one way to identify such regions which have the potential to act as the 189 genetic origin of rapid responses to the ever-changing environmental challenges. 190

Plant material and whole genome sequencing 192
The seeds of all seven accession were received from Maarten Koornneef (MPI for Plant 193 Breeding Research), and were grown under normal greenhouse conditions. DNA preparation 194 and next generation sequencing was performed by the Max Planck Genome center. The DNA 195 samples were prepared as previously described for the Ler 9 , and sequenced with a PacBio 196 Sequel system. For each accession, data from two SMRT cells were generated. Besides,197 Illumina paired-end libraries were prepared and sequenced on the Illumina HiSeq2500 198 systems. 199

Genome assembly 200
PacBio reads were filtered for short (<50bp) or low quality (QV<80) reads using 201 SMRTLink5 package. De novo assembly of each genome was initially performed using three 202 different assembly tools including Falcon 47 , Canu 48 and MECAT 49 . The resulting assemblies 203 were polished with Arrow from the SMRTLink5 package and then further corrected with 204 mapping of Illumina short reads using BWA 50 to remove small-scale assembly errors which 205 were identified with SAMTools 51 . For each genome, the final assembly was based on the 206 Falcon assembly as these assemblies always showed highest assembly contiguity. A few 207 contigs were further connected or extended based on whole genome alignments between 208 Falcon and Canu or MECAT assemblies. Contigs were labelled as organellar contigs if they 209 showed alignment identity and coverage both larger than 95% when aligned against the 210 mitochondrial or chloroplast reference sequences. A few of contigs aligned to multiple 211 chromosomes and were split if no Illumina short read alignments supported the conflicting 212 regions. Assembly contigs larger than 20kb were combined to pseudo-chromosomes 213 according to their alignment positions when aligned against the reference sequence using 214 MUMmer4 52 . Contigs with consecutive alignments were concatenated with a stretch of 500 215

Assembly evaluation 217
We evaluated the assembly completeness by aligning the reference genes against each 218 of the seven genomes using Blastn 53 . Reference genes which were not aligned or only 219 partially aligned might reveal genes which were missed during the assembly. To examine 220 whether they were really missed, we mapped Illumina short reads from each genome against 221 the reference genome using the BWA 50 and checked the mapping coverage of these genes. 222 The genes, which were missing in the assembly, should show fully alignment coverage 223 (Supplementary Table 7). 224 Centromeric and telomeric tandem repeats were annotated by searching for the 178 bp 225 tandem repeat unit 54 and the 7 bp tandem repeat unit of TTTAGGG 55 . rDNA clusters were 226 annotated with Infernal version 1.1 56 . 227 The assembly contiguity of Cvi-0 and Ler were further tested using three previously 228 published genetic maps 57-59 (Supplementary Table 3). For this we aligned the marker 229 sequences against the chromosome-level assemblies and checked the order of the markers 230 in the assembly versus their order in the genetic map. The ordering of contigs to 231 chromosomes was perfectly supported by all three maps. Overall, only six (out of 1,156) 232 markers showed conflicts between the genetic and physical map. In all six cases we found 233 evidence that the conflict was likely caused by structural differences between the parental 234 genomes. 235

Gene annotation 236
Protein-coding genes were annotated based on ab initio gene predications, protein 237 sequence alignments and RNA-seq data. Three ab initio gene predication tools were used 238 including Augustus 60 , GlimmerHMM 61 and SNAP 62 . The reference protein sequences from the 239 Araport 11 16 annotation were aligned to each genome assembly using exonerate 63 with the 240 parameter setting "--percent 70 --minintron 10 --maxintron 60000". For five accessions (An-1, 241 C24, Cvi-0, Ler-0, and Sha) we downloaded a total of 155 RNA-seq data sets from the NCBI 242 SRA database (Supplementary Table 8). RNA-seq reads were mapped to the corresponding 243 genome using HISAT2 64 and then assembled into transcripts using StringTie 65 (both with 244 default parameters). All different evidences were integrated into consensus gene models 245 using Evidence Modeler 66 . 246 The resulting gene models were further evaluated and updated using the Araport 11 16 247 annotation. Firstly, for each of the seven genomes, the predicted gene and protein sequences 248 were aligned to the reference sequence, while all reference gene and protein sequences were 249 aligned to each of the other seven genomes using Blast 53 . Then, potentially mis-annotated 250 8 genes including mis-merged (two or more genes are annotated as a single gene), mis-split 251 (one gene is annotated as two or more genes) and unannotated genes were identified based 252 on the alignments using in-house python scripts. Mis-annotated or unannotated genes were 253 corrected or added by incorporating the open reading frames generated by ab initio 254 predications or protein sequence alignment using Scipio 67 . 255 Noncoding genes were annotated by searching the Rfam database 68 using Infernal 256 version 1.1 56 . Transposon elements were annotated with RepeatMasker 257 (http://www.repeatmasker.org). Disease resistance genes were annotated using RGAugury 69 . 258 NB-LRR R gene clusters were defined based on the annotation from a previous study 70 . 259

Pan-genome analysis 260
Pan-genome analyses were performed at both sequence and gene level. To construct a 261 pan-genome of sequences, we generated pair-wise whole genome sequence alignments of 262 all possible pairs of the eight genomes using the nucmer in the software package MUMmer4 52 . 263 A pan-genome was initiated by choosing one of the genomes, followed by iteratively adding 264 the non-aligned sequence of one of the remaining genomes. Here, non-aligned sequences 265 were required to be longer than 100bp without alignment with an identity of more than 90%. 266 The core genome was defined as the sequence space shared by all sampled genomes. Like 267 the pan-genome, the core-genome analysis was initiated with one genome. Then all other 268 genomes were iteratively added, while excluding all those regions which were not aligned 269 against each of the other genomes. The pan-and core-genome of genes was built in a similar 270 way. The pan-genome of genes was constructed by selecting the whole protein coding gene 271 set of one of the accessions followed by iteratively adding the genes of one of the remaining 272 accessions. Likewise, the core-genome of genes was defined as the genes shared in all 273 sampled genomes. 274 For each pan or core genomes analysis, all possible ( integrating the eight genomes (or a subset of them) were evaluated. The exponential 276 regression model y = A e Bx + C to was then used to model the pan-genome/core-genomes by 277 fitting medians using the least square method implemented in the nls function of R. 278

Analysis of structural rearrangements and gene CNV 279
All assemblies were aligned to the reference sequence using nucmer from the 280 MUMmer4 52 toolbox with parameter setting "-max -l 40 -g 90 -b 100 -c 200". The resulting 281 alignments were further filtered for alignment length (>100) and identity (>90). Structural 282 rearrangements and local variations were identified using SyRI 9 . The functional effects of 283 sequence variation were annotated with snpEff 71 . The gene CNV were identified according to 284 the gene family clustering using the tool OrthoFinder 21 based on all protein sequences from 285 9 the eight accessions. 286

Synteny diversity, hotspots of rearrangements and diversity estimates 287
Synteny Diversity was defined as the average fraction of non-syntenic sites found within 288 all pairwise genome comparisons within a given population. Here we denote Synteny 289 where . and / refer to the frequencies of sequence and and ./ to the average 292 probability of a position in to be non-syntenic. Note, p syn can be calculated in a given region 293 or for the entire genome. However even when calculated for small regions the annotation of 294 synteny still needs to established within the context of the whole genomes to avoid false 295 assignments of homologous but non-allelic sequence. Here we used the annotation of SyRI 296 to define syntenic regions. p syn values can range from 0 to 1, with higher values referring to a 297 higher average degree of non-syntenic regions between the genomes. 298 For the analyses, we calculated p syn in 5-kb sliding windows with 1kb step-size across the 299 entire genome. HR regions were defined as regions with p syn larger than 0.5. Neighboring 300 regions were merged into one HR if their distance was shorter than 2kb. 301 The nucleotide and haplotype diversity were calculated with the R package PopGenome 72 302 using SNP markers (with MAF > 0.05) from 1001 Genomes Project 4 . LD were calculated as 303 correlation coefficients r 2 using SNP markers with MAF > 0.05. GO     Each point corresponds to a pan-or core-genome size estimated with a particular 540 the same number of genomes. Pan-genome (blue) and core-genome (red) estimations were 542 fitted using an exponential model. (g) Gene copy number variations (CNV) between the 543 reference and each of the accessions group by differences in the gene families: n vs. n: both 544 the reference and the query genomes have the same number of genes in a gene family. 545 Others categories (1 vs. 2, 2 vs. 1, 1 vs. 3, 3 vs. 1, 2 vs. 3 and 3 vs 2)