DNA extraction and sequencing<\/p>\n
For long-read sequencing, we began with 3-week-old plants grown in soil that had been transferred to darkness for 24\u201348\u2009h before harvesting to reduce the starch content. A total of 20\u201330\u2009g of flash-frozen rosette tissue, pooled from individuals, was ground in liquid nitrogen with a pestle and mortar. Nuclei were isolated as described for accession Ey15-2 (ref. 75<\/a>), and high-molecular-weight (HMW) DNA purified with Genomic-tips 100G (Qiagen, 10243) following the manufacturer\u2019s instructions. Ten micrograms of HMW DNA were sheared with either Megaruptor 3 (Diagenode, B06010003) or a FINE-JECT 26G\u00d7 1\u2033 needle (0.45\u2009\u00d7\u200925\u2009mm; 14-13651) to ca. 75\u2009kb, and used as input for long-read library preparation with the SMRTbell Express Template Preparation Kit 2.0 (Pacific Biosciences, 101-693-800). These libraries were size-selected with the BluePippin system (Sage Science) with a 30\u2009kb cutoff in a 0.75% DF Marker U1 high-pass 30\u201340\u2009kb vs3 gel cassette (Biozym, BLF7510). Libraries for accessions 9981 (Angit-1; CS76366) and 10002 (TueWal-2; CS76405) were sequenced on a Sequel II system (Pacific Biosciences), and the others on a Sequel I system.<\/p>\n To prepare PCR-free libraries for short-read sequencing, the genomic DNA was fragmented to 250\u2013350\u2009bp using a Covaris S2 Focused Ultrasonicator (Covaris). The libraries were prepared according to the manufacturer\u2019s instructions with either the TruSeq DNA PCR-free kit (Illumina, 20015962) or the NxSeq AmpFREE Low DNA Library kit (Lucigen, 14000-2). In total, libraries for 89 accessions (including the main 27 for which we assembled their genomes) were sequenced in paired-end mode on a HiSeq 3000 system (Illumina).<\/p>\n The ultra-HMW DNA extraction and sample preparation for optical maps were performed as described75<\/a>,76<\/a> at Corteva Agriscience using the Direct Label and Stain technology (Bionano Genomics).<\/p>\n Assembly<\/p>\n The CLR subreads were assembled with Canu v1.71 (ref. 77<\/a>). Since accessions 9981 and 10002 had been sequenced at higher coverage on a Sequel II instrument, only about 200\u00d7 genome coverage worth of reads were used for assembly. We performed two rounds of polishing on the resulting contigs of all assemblies\u2014first with the CLR subreads and Arrow v2.3.2 (https:\/\/github.com\/PacificBiosciences\/gcpp<\/a>), and then with PCR-free short reads and Pilon v1.22 (ref. 78<\/a>).<\/p>\n For scaffolding, we generated hybrid scaffolds with optical maps for eight accessions (Supplementary Table 1<\/a>) using Bionano Access v1.5 and Bionano Solve<\/a> v3.6. The assembly was performed in pre-assembly mode with parameters nonhaplotype and no-CMPR-cut, without extend-split. Based on what we learned from these hybrid assemblies, we set the parameters for in silico scaffolding of the other genomes. We scaffolded contigs >150\u2009kb with RagTag v1.1.1 (ref. 79<\/a>; scaffold -q 60 -f 10000 -I 0.6 -remove-small) using the TAIR10 reference with hard-masked centromeres, rDNAs, telomeres and nuclear insertions of organelles to prevent misplacement of contigs due to reference bias75<\/a>. All scaffolded assemblies were manually curated to specifically discard low-confidence centromere satellite-rich contigs or to invert contigs with satellite repeats at their edges, indicative of their correct orientation. These edits were implemented in the AGP files, which were converted to FASTA format with the RagTag agp2fa function79<\/a>. To detect traces of residual heterozygosity, we aligned the original long reads to their corresponding chromosome scaffolds using pbmm2 v1.3.0 with the parameters align -sort -log-level DEBUG -preset SUBREAD -min-length 5000. Unmapped reads, as well as secondary and supplementary alignments, were filtered out using samtools v1.9 (view -b -F 2308 Chr1 Chr2 Chr3 Chr4 Chr5). The resulting BAM file was then analyzed with NucFreq v0.1 (-minobed 2) to assess genome-wide coverage of primary and secondary alleles80<\/a>. AGP files, both before and after manual curation, as well as NucFreq plots, are available in the GitHub repository of this project.<\/p>\n Repeat annotation<\/p>\n Repetitive elements were annotated as described75<\/a>. We ran RepeatMasker v4.0.9 (-cutoff 200 -nolow -gff -xsmall) using a custom library that included various consensus sequences for the CEN178 (ref. 81<\/a>), 5S rDNA82<\/a>, 45S rDNA83<\/a> and telomere repeats. We annotated tRNAs with tRNAscan-SE v2.0.6 (ref. 84<\/a>) and TEs with Extensive de novo TE Annotator v1.9.7 (ref. 85<\/a>; \u2013step all \u2013sensitive 1 \u2013anno 1), a pipeline that combines several TE annotation tools (LTRharvest, LTR_FINDER, LTR_retriever, TIR-Learner, HelitronScanner and TEsorter)86<\/a>,87<\/a>,88<\/a>,89<\/a>,90<\/a>,91<\/a>,92<\/a>,93<\/a>. Finally, to understand the causes of contig breaks, we determined the type of repetitive element closest to each contig edge, considering the first 2\u2009kb from each edge in contigs >10\u2009kb.<\/p>\n Pannagram<\/p>\n Pannagram is a toolkit designed for reference-free pangenome alignment, annotation and analysis, as well as for generating diagrams46<\/a>.<\/p>\n We represent the WGA as a matrix of corresponding positions, where rows represent accessions and columns represent homologous positions. The construction of the alignment is done in a reference-free manner (see below). However, to visualize the alignment in genome browsers, columns must be sorted in some manner, for example, to correspond to the TAIR10 sequence order. Then, columns of the pangenome are used as positions in the pangenome coordinate system.<\/p>\n To perform reference-free WGA, we developed a three-step pipeline. First, we use several accessions as references and build draft pairwise alignments between each and all other accessions. This process results in several reference-based matrices of corresponding positions. Next, we intersect these matrices, selecting only those columns that are present in all reference-biased matrices, which produces reliable and reference-independent correspondences. In the final step, we resolve unaligned sequences between blocks of corresponding positions using multiple sequence alignment tools. Once the reference-free alignment is complete, it can be sorted according to the desired order of accessions. In our case, we employ an alphabetical order, with the TAIR10 genome first.<\/p>\n For the pairwise alignments between a reference genome (not necessarily TAIR10) and another accession, the focal accession genome is divided into blocks of 5,000\u2009bp, and each block is then mapped to the corresponding chromosomes of the reference genome using BLAST, with exactly one best hit retained for each block through this process. Next, the BLAST hits that are not in close proximity to each other in both genomes are removed. An additional BLAST search is performed to align corresponding unaligned sequences between remaining hits.<\/p>\n To resolve any unaligned blocks after the reference-randomization procedure, MAFFT94<\/a> is used. Blocks longer than 30\u2009kb cannot be aligned within a reasonable time using MAFFT, so they are considered to be highly diverged. We found the final unaligned regions to be primarily associated with centromeric regions, rDNA clusters, telomeres and complex regions of multiple and long insertions and deletions, which are regions that are not of primary interest in this paper.<\/p>\n Given the WGA, SNPs can simply be output as sequence differences. However, sequence differences can arise from ambiguities in local alignment and do not necessarily correspond to SNPs (Supplementary Note 7<\/a>). If we consider all sequence differences as SNPs, a pair of accessions differs at over 800,000 positions on average; however, if we restrict ourselves to isolated sequence differences, the number shrinks to 600,000.<\/p>\n Pangenome graphGraph construction<\/p>\n We constructed genome graphs for each of the five chromosomes using the PGGB pipeline29<\/a>. First, we prepared the assemblies by splitting them into chromosomes and removing all unplaced contigs. To enforce linearity for simpler analysis and comparison, we used a modified version of accession 22001 with the genome rearranged to a consensus pan-genomic order (suffix: \u2018f\u2019). We added the TAIR10 reference genome to the graph to enable anchoring and presentation of results in a reference framework.<\/p>\n We executed the PGGB pipeline (downloaded on 25 January 2024) with the following parameters: -s 10000 -p 90 -n 27. PGGB consists of the following three methods: an all-against-all alignment with wfmash (v0.12.4-5-g0b191bb), graph induction using seqwish (v0.7.9-2-gf44b402) and two rounds of pangenome ordering (odgi v0.8.3-26-gbc7742ed) followed by normalization with smoothxg (v0.7.2-11-g9970e0d). The graph was used for analyzing the pangenome and synteny, as well as detecting variation using vg deconstruct95<\/a>.<\/p>\n Similarity<\/p>\n We exploited graph properties to classify different levels of similarity between genomes. Nodes traversed in all accessions are labeled as core, nodes traversed in only one accession are private and all other nodes (>1 and <\/p>\n Synteny windows<\/p>\n Every node in the graph can be translated to its exact position for each path. This direct connection allows us to create sliding window approaches for each sample\/path using graph-based statistics. Here we used nonoverlapping windows of 300\u2009kb and calculated the average similarity (see above) of these regions. This was performed for each graph and path independently and the results were represented in a heat map.<\/p>\n Saturation analysis<\/p>\n A saturation analysis was performed using a bootstrapping approach. In each iteration, we removed a specific number of paths from our graph and performed the same pangenome categorization as above (\u2018Similarity\u2019). In addition, we added the total pangenome, which describes the total amount of sequence (core\u2009+\u2009shell\u2009+\u2009private sequence). We performed 20 different (unique) combinations for each size (number of genomes).<\/p>\n Deconstructing the graph<\/p>\n To achieve full insights into graph variation and cover all bubbles in the graph, VG deconstruct was run multiple times with each accession reference path once (vg deconstruct -a -e). After, the reported VCF (v1.54.0 \u2018Parafada\u2019) files were converted to a BED file with all important information provided. In addition, each chromosome was merged, and the genotype information was concluded and added. Bubbles were identified by the start and end positions, and all traversals within these bubbles were also reported. Scripts can be found in the repository.<\/p>\n sSVs and cSVs in the graphs were defined as follows:<\/p>\n All SVs represent indels, having one very small traversal (deletion) and a large one containing the SV sequence (insertion).<\/p>\n<\/li>\n Bubbles were identified as sSVs if the bubble was shared by all accessions in the graph (here 28), and as cSVs if not. Traversals covering the insertion are at least 15\u2009bp long and must exhibit high similarity (95% sequence). The deletion part of the bubble should be small, at most 5% of the length of the inserted sequence.<\/p>\n<\/li>\n Most cSVs correspond to bubbles that have a complex structure and\/or are sub-bubbles of larger bubbles.<\/p>\n<\/li>\n<\/ul>\n General pangenome<\/p>\n To perform a reference-free pangenome analysis, we used genome graphs built separately for each chromosome. The complete graph contains 18.3 million nodes and 20.9 million edges, with a total size of 225\u2009Mb, and has a mean compression rate of 6.75% across all chromosomes. Similar to other genome-wide analyses in this study, the large-scale reciprocal translocation in accession 22001 was masked to maximize linearity and increase resolution in the variation graph.<\/p>\n The mobile-ome<\/p>\n The mobile-ome refers to the collection of insertions and deletions that are likely to have occurred recently and are therefore not fixed in our sample. We hypothesize that each mobile event results in an SV, specifically a presence\u2013absence polymorphism at the location of the insertion or deletion. Consequently, our initial approach involves extracting all presence\u2013absence SVs and systematically decomposing them step by step. To distinguish between simple bi-allelic presence\u2013absence polymorphisms (indels) and cSV, we analyzed the lengths of alleles within the SVs. We distinguish two types based on the similarity threshold s, with s\u2009=\u20090.9 in our case. We consider a simple indel as one that contains alleles of two length types\u2014those that are shorter than (1\u2009\u2212\u2009s) of the SV length (absence allele) and those that are longer than s of the SV length (presence allele). The distinction between simple and complex presence\u2013absence polymorphisms is partially a computational construct to filter SVs and simplify further analysis. Simple indels and complex presence\u2013absence polymorphisms form a continuum, and by relaxing the similarity threshold (s\u2009A. thaliana TEs, as well as against themselves. The indels that exhibited some similarity to known TEs were divided into the following groups: is complete\u2014significant similarity to known TEs and can be classified as TEs themselves; contains complete\u2014contained regions with similarity to known TEs, but also additional sequences; is fragment\u2014contained only partially sequenced with similarity to known TEs; and contains fragment\u2014partial coverage by BLAST hits of TE segments, but also additional sequences unrelated to known TEs.<\/p>\n We consider all these indels as parts of the mobile-ome. Indels without similarity to known TEs but showing nested similarities within the indel data set (where one sequence is a subsequence of another) were considered as potential candidates for new mobile-ome elements. To investigate their potential function, we obtained all six open reading frames within each of these indels. From each translated sequence, we selected either all continuous stretches without stop codons that were longer than 100 codons or the longest stretch that exceeded 30 codons without a stop. Subsequently, we performed a BLAST search using the obtained amino acid sequences against the NCBI protein database and classified the potential proteins into four categories. If the BLAST results for an sSV contained keywords related to TE, we assigned the sequence to the TE-like category. These keywords were \u2018transcriptase\u2019, \u2018reverse\u2019, \u2018transpos\u2019, \u2018gag-\u2019, \u2018pol-\u2019, \u2018integrase\u2019, \u2018gag\/pol\u2019, \u2018gagpol\u2019, \u2018retrovirus\u2019, \u2018RNA-directed DNA polymerase\u2019 and \u2018RNA-dependent DNA polymerase\u2019. sSVs that only had BLAST hits with descriptions such as \u2018hypothetical protein\u2019, \u2018unnamed protein product\u2019, \u2018uncharacterized protein\u2019, \u2018predicted protein\u2019, \u2018PREDICTED:\u2019, \u2018putative protein\u2019 and \u2018unknown\u2019 were categorized as \u2018undefined proteins\u2019. Indels without any BLAST hits were classified as \u2018no protein\u2019. In all other cases, sSV was categorized as a \u2018defined protein.\u2019<\/p>\n Gene annotationPreliminary annotation<\/p>\n Gene annotation was mainly based on Augustus (v3.3.3)96<\/a>. Augustus-predicted gene models were trained using parameters obtained from \u2018hints\u2019 from three different sources. First, we ran BUSCO (v4.0.1)97<\/a> with -m genome option. Second, the A. thaliana reference gene annotation was projected onto each genome using Liftoff98<\/a> with the -exclude_partial and -copies options. Third, the RNA-seq data for each accession were used\u2014wiggle hints were generated using bam2wig and wig2hints, and EST hints were generated using bam2hints (all three tools provided by Augustus). Augustus was run with the following nondefault parameters:<\/p>\n \n \u2013softmasking 1<\/p>\n<\/li>\n \n \u2013species=BUSCO_retraining<\/p>\n<\/li>\n \n \u2013gff3=on<\/p>\n<\/li>\n \n \u2013extrinsicCfgFile=Custom_Config<\/p>\n<\/li>\n \n \u2013hintsfile=Liftoff_hints<\/p>\n<\/li>\n<\/ul>\n For every accession, the GFF3 output of Augustus was run through the Augustus-provided tool getAnno.pl to translate gene annotations into protein sequences. Finally, for each annotation, the Augustus output was combined and evaluated using augustus_GFF3_to_EVM_GFF3.pl (provided by EVidenceModeler99<\/a>).<\/p>\n In addition to the Augustus-generated annotations, we used two types of independent evidence for gene models: from the SNAP de novo annotation tool100<\/a> and Cufflinks transcriptome assemblies101<\/a>. Annotations produced by Augustus, SNAP and Cufflinks were combined and then subdivided into 1-Mb windows with 1-kb overlap using partition_EVM_input.pl (provided by EVidenceModeler). We ran EVidenceModeler with annotation GFF files, the assembly fasta file, the partitions and a weight matrix. We chose weights for each input based on their ability to recreate the Araport11 gene annotation. Running EVidenceModeler produced the final annotation compilation for each accession. We retained only the longest isoform for each gene using gffread102<\/a>.<\/p>\n Reconciling annotations<\/p>\n To enable comparison between the independent annotations, we used the pangenome coordinate system, reconciling discrepancies using majority voting (Supplementary Note 6<\/a>\u2014\u2018Details about reconciling annotations and gene filtering\u2019). Additionally, we compared the sequences of each gene across different accessions. If a gene showed significant variation because it was located in regions heavily influenced by SVs, we excluded it from the analysis. In total, 3,438 genes in our annotation were the result of splitting preliminary annotations and 1,020 were the result of merges. Lastly, we added 1,789 TAIR10 genes that had not been detected by our annotation pipeline (the likely reason for which is that our RNA-seq data only covered four tissues\/stages) to our annotation. For these genes, the same pangenome coordinate approach was used to map the TAIR10 annotation of the 1,789 added genes into their annotations in other genomes. Our approach generated a total of 34,153 putative genes. For the details of annotation reconciliation and filtering, see Supplementary Note 6<\/a>\u2014\u2018Details about reconciling annotations and gene filtering\u2019.<\/p>\n Ancestry analysis<\/p>\n All PC sequences from all accessions were compared using DIAMOND\u2019s blastp module103<\/a> (version 2.0.11) against the A. lyrata MN47 proteome (version 2, GenBank: GCA_944990045.1), and the best hit was considered as the A. lyrata homolog. To avoid bias due to mis-annotated genes in the A. lyrata proteome, we further applied Liftoff v1.63 (ref. 98<\/a>) to annotate all A. thaliana genes from all accessions on A. lyrata MN47 (v2, https:\/\/doi.org\/10.6084\/m9.figshare.22285444.v1<\/a>) and A. lyrata NT1 (v2, https:\/\/doi.org\/10.6084\/m9.figshare.22293196.v1<\/a>) assemblies. Next, each annotation group from A. thaliana was assigned to the A. lyrata homolog (by LiftOff or proteome similarity) that was common to at least 50% of its members, sharing at least 80% identity and covering at least 80% of the A. thaliana coding sequence. A. thaliana annotation groups were defined to be ancestrally relative to the A. lyrata gene if they were part of a colinear segment of at least two genes. To that end, all A. thaliana genes were ordered according to their relative position in the pangenome coordinate system. Each pair of consecutive genes in A. thaliana was assigned to the same colinear segment as its homologs in A. lyrata if the homologs were separated by fewer than six genes. The ancestral state was defined as \u2018similar\u2019 for cases where the genes from A. lyrata and A. thaliana were not part of the same colinear segment but shared at least 80% sequence identity over at least 80% of the length A. thaliana gene. Further details are available in Supplementary Note 6<\/a>\u2014under \u2018Genes and TEs\u2019 for TE analysis and \u2018New genes\u2019 for the origin of new genes.<\/p>\n Expression analysisRNA-seq read mapping and gene expression calculation<\/p>\n Raw RNA-seq reads from 7-day-old seedlings, 9-leaf rosettes, flowers and pollen54<\/a> were aligned either to the TAIR10 reference genome or the corresponding accession accession genome using STAR v2.7.1 (ref. 104<\/a>) with the following custom options:<\/p>\n \n \u2013alignIntronMax 6000<\/p>\n<\/li>\n \n \u2013alignMatesGapMax 6000<\/p>\n<\/li>\n \n \u2013outFilterIntronMotifs RemoveNoncanonical<\/p>\n<\/li>\n \n \u2013outFilterMismatchNoverReadLmax 0.1<\/p>\n<\/li>\n \n \u2013outFilterMismatchNmax 999<\/p>\n<\/li>\n \n \u2013outFilterMismatchNoverLmax 0.3<\/p>\n<\/li>\n \n \u2013outFilterMultimapNmax 1<\/p>\n<\/li>\n \n \u2013alignSJoverhangMin 8<\/p>\n<\/li>\n \n \u2013outSAMattributes NH HI AS nM NM MD jM jI XS<\/p>\n<\/li>\n<\/ul>\n Read alignment statistics are provided in Supplementary Table 4<\/a>. Expression levels were assessed using featurecounts from Subread v2.0.1 (ref. 105<\/a>) on each RNA-seq sample with either the TAIR10 gene annotation or the accession-specific annotations from this study. The entire locus, including exons and introns, was used for expression estimates. Expression levels were normalized by calculating TPMs, which represent the number of read counts divided by the gene length in kilobases, and then dividing the total number of counts per kilobase for all genes by 1 million.<\/p>\n Mapping to TAIR10 versus the own genome<\/p>\n To determine whether the gene expression calculation was consistent between RNA-seq mapping in TAIR10 versus accession-specific genomes, we focused on the annotation groups with a one-to-one correspondence with an Araport11 gene. For each RNA-seq sample, we obtained the Pearson\u2019s correlation coefficient between the number of exonic counts obtained from TAIR10 mapping and accession-specific mapping. We also determined the number of genes that were correctly or wrongly estimated using TAIR10 mapping. We called a gene \u2018wrong\u2019 if the counts in TAIR10 and the counts in its own genome differed by more than 30% (Ncounts_min\/Ncounts_max \u2264 0.7). Only genes with at least six counts in either calculation were analyzed.<\/p>\n Chromatin immunoprecipitation followed by sequencing analysis<\/p>\n We used chromatin immunoprecipitation followed by sequencing (ChIP\u2013seq) data from 6 accessions and sRNA-seq data from 14 accessions54<\/a>. We used STAR104<\/a> to map ChIP\u2013seq reads with these nondefault options:<\/p>\n \n \u2013alignIntronMax 5<\/p>\n<\/li>\n \n \u2013outFilterMismatchNmax 10<\/p>\n<\/li>\n \n \u2013outFilterMultimapNmax 1<\/p>\n<\/li>\n \n \u2013alignEndsType EndToEnd<\/p>\n<\/li>\n \n \u2013twopassMode Basic<\/p>\n<\/li>\n<\/ul>\n The ChIP\u2013seq data were log2-normalized to input using bamCompare (deeptools package106<\/a>) using<\/p>\n The ChIP\u2013seq coverage was estimated using bedtools map-mean107<\/a>. The ChIP\u2013seq coverage was further normalized to obtain value range similarity across accessions. For this, we applied quantile-normalization using an R function:<\/p>\n function(x) { (x-quantile(x,.20)) \/ (quantile(x,.80) – quantile(x,.20)) }<\/p>\n<\/li>\n<\/ul>\n which equalized the 20% and 80% quantile values of each ChIP\u2013seq sample. After quantile-normalization, the replicated samples were averaged.<\/p>\n sRNA-seq analysis<\/p>\n We used sRNA-seq data for 14 accessions54<\/a>. To process the sRNA-seq data, we trimmed the reads using cutadapt108<\/a>: cutadapt -a AACTGTAGGCACCATCAAT \u2013minimum-length 18. We then used STAR104<\/a> with the following nondefault options to map sRNA-seq reads to the corresponding genome:<\/p>\n \n \u2013runRNGseed 12345<\/p>\n<\/li>\n \n \u2013alignEndsType Extend5pOfRead1<\/p>\n<\/li>\n \n \u2013alignIntronMax 5000 \u2013alignSJDBoverhangMin 1<\/p>\n<\/li>\n \n \u2013outReadsUnmapped Fastx \u2013outSAMmultNmax 100<\/p>\n<\/li>\n \n \u2013outSAMprimaryFlag AllBestScore<\/p>\n<\/li>\n \n \u2013outSAMattributes NH HI AS nM NM MD jM jI XS<\/p>\n<\/li>\n \n \u2013outFilterMultimapNmax 10<\/p>\n<\/li>\n \n \u2013outFilterMatchNmin 16<\/p>\n<\/li>\n \n \u2013outFilterMatchNminOverLread 0.66<\/p>\n<\/li>\n \n \u2013outFilterMismatchNmax 2<\/p>\n<\/li>\n \n \u2013outFilterMismatchNoverReadLmax 0.05<\/p>\n<\/li>\n \n \u2013outFilterIntronMotifs RemoveNoncanonicalUnannotated<\/p>\n<\/li>\n \n \u2013twopassMode None<\/p>\n<\/li>\n<\/ul>\n We extracted 24-nt reads, calculated read coverage for each position of the genome using genomeCoverageBed (bedtools v.2.27.1), normalized it by the total number of uniquely mapped reads in each sample, and calculated 24-nt sRNA coverage for each locus of interest using bedtools map -mean function.<\/p>\n DNA methylation analysis<\/p>\n To estimate DNA methylation levels, we used published BS-seq data for 12 accessions53<\/a>. After trimming with TrimGalore (https:\/\/github.com\/FelixKrueger\/TrimGalore<\/a>) with \u2013clip_r1 10 \u2013clip_r2 15 \u2013three_prime_clip_r1 10 \u2013three_prime_clip_r2 10, reads for each accession were mapped to its corresponding genome with Bismark109<\/a> with \u2013score_min L,0,-0.5 for a relaxed mismatch threshold and the \u2013un \u2013ambiguous parameters to obtain additional unmapped and multiply-mapping reads. Methylation was called as described110<\/a>. CG, CHG and CHH methylation levels for genes and SVs in each accession were then calculated for each gene by focusing on all Cs in the specific context within the gene and calculating the ratio between the total number of methylated and unmethylated reads across all sites.<\/p>\n Mapping to TAIR10 versus own genome<\/p>\n To estimate reference bias, we mapped BS-seq data for all accessions to the TAIR10 genome and performed CG, CHG and CHH methylation level estimation in the same way as for own genomes. We then focused on annotation groups with a one-to-one correspondence with an Araport11 gene (the current annotation of the TAIR10 genome). We calculated Pearson\u2019s correlation coefficient between the methylation level estimates obtained from TAIR10 mapping and accession-genome mapping. We also estimated the number of genes that were correctly or wrongly estimated using TAIR10 mapping. For each methylation context, we called a gene \u2018wrongly estimated\u2019 if the methylation level in TAIR10 and the own genome differed by more than 50% (methlevel_min\/methlevel_max\u2009\u2264\u20090.5). For a more refined analysis of reference bias, see Supplementary Note 8<\/a>.<\/p>\n Reporting summary<\/p>\n\n
\n
\n
\n
\n
\n