Telomere-to-telomere genome assembly of Phaeodactylum tricornutum

Growth

Phaeodactylum tricornutum (Culture Collection of Algae and Protozoa CCAP 1055/1) was grown in L1 medium without silica at 18 °C under cool white fluorescent lights (75 mE m⁻² s⁻¹) and a photoperiod of 16 h light:8 h dark as described previously (Slattery et al., 2018).

DNA extraction

200 mL of culture (approximately 5 ×10⁸ cells) was spun at 3000×g for 10 min at 4 °C. The pellet was resuspended in one mL TE (pH 8.0) and added dropwise to a mortar (pre-cooled at −80 °C) pre-filled with liquid nitrogen. The frozen droplets were ground into a fine powder with a mortar and pestle, being careful to keep the cells from thawing by adding more liquid nitrogen as necessary. The frozen powder was transferred to a 15 mL Falcon tube where two mL of lysis buffer was added (1.4 M NaCl, 200 mM Tris–HCl pH 8.0, 50 mM EDTA, 2% (w/v) CTAB, RNAse A (250 µg/mL) and proteinase K (100 µg/mL)). The solution was mixed very slowly by inversion, incubated for 30 min at 37 °C (mixed very slowly halfway through incubation). Cellular debris was pelleted at 6000×g for 5 min. Lysate was transferred to a new 15 mL Falcon tube. One volume of 25:24:1 phenol:chloroform:isoamyl alcohol was added, mixing slowly by inversion. The sample was centrifuged at 6000×g for 5 min. The aqueous phase was transferred as slow as possible to a new Falcon tube. One volume of 24:1 chloroform:isoamyl alcohol was added, and mixed slowly with end-over-end inversion. The sample was centrifuged at 6000×g for 5 min. Approximately 450 uL of the aqueous phase was transferred into new 1.5 mL Eppendorf tubes. To the Eppendorf tube, 1/10 volume of 3 M NaAc pH 5.2 and two volumes (final volume) of ice-cold 100% ethanol were added, mixing slowly by end-over-end inversion. The sample was centrifuged at 16,000×g for 5 min, and washed twice with 500 uL 70% ethanol. Ethanol was decanted, and the pellet was dried for approximately 10 min by inverting on a paper towel. The pellet was resuspended in 100 uL 10 mM Tris–HCl pH 8.0, 0.1 mM EDTA pH 8.0. After resuspending overnight at 4 °C, DNA fragments less than 20 kbp were then selectively removed using the Short Read Eliminator (SRE) kit from Circulomics (Baltimore). DNA from the same extraction was used for sequencing on both the Oxford Nanopore MinION and Illumina NextSeq 550 platform.

Sequencing

An Oxford Nanopore MinION flow cell R9.4.1 was used with the SQK-LSK109 Kit according to the manufacturer’s protocol version GDE_9063_v109_revK_14Aug2019, with one alteration: for DNA repair and end-prep, the reaction mixture was incubated for 15 min at 20 °C and 15 min at 65 °C. Basecalling was performed after the run with Guppy (Version 3.6). NanoPlot (De Coster et al., 2018) was used to generate Q-score versus length plots and summary statistics. The read N50 of the unfiltered reads was approximately 35 kb (Supplemental Figure S1).

For Illumina sequencing, the Nextera XT kit was used, and a 2X75 paired-end mid-output NextSeq 550 library was prepared according to the manufacturer’s protocol, and sequenced at the London Regional Genomics Center (lrgc.ca). Reads were trimmed using Trimmomatic v0.36 (Bolger, Lohse & Usadel, 2014) in paired-end mode with the following settings: AVGQUAL:30 CROP:75 SLIDINGWINDOW:4:25 MINLEN:50 TRAILING:15. SLIDINGWINDOW AND TRAILING were added to remove poor quality base calls.

Assembly

Methylation

A total of 5mC methylation sites were predicted using Megalodon v2.2.1 (Oxford Nanopore Technologies) using the model res_dna_r941_min_modbases_5mC_CpG_v001.cfg from the Rerio repository (Oxford Nanopore Technologies) with Guppy 4.5.2. A default threshold of 0.75 was used as a minimum score for modified base aggregation (probability of modified/canonical base) to produce the final aggregated output. The percentage of reads methylated at the predicted locations are plotted in Supplemental Figure S2.

Workflow

We developed a sample preparation protocol that provided high-molecular weight DNA. We observed a read N50 of 35 kilobases, with the longest reads just over 300 kb, following sequencing with the Oxford Nanopore MinION platform. Of the 7.8 gigabases of raw sequence data, approximately 2.5 gigabases were from reads longer than 50 kilobases (Supplemental Figure S1). We found that chromosomes assembled using standard approaches were often mis-assembled around telomeres, or were fragmented and only contained 1 telomere. To correct each contig, we used the unique ultra-long telomere reads as described in the Supplemental File S1 and in Fig. 1. This approach was used to manually identify a tiling path for each chromosome until each chromosome was contiguous from telomere to telomere, and validated by a tiling overlapping read path.

Telomere-to-telomere assembly comprises previous scaffolds

We ultimately obtained 25 telomere-to-telomere chromosome assemblies that recruit 98% of long reads, and these chromosomes comprise all previously proposed chromosomes from Bowler et al. (2008), as well as circularized chloroplast and mitochondrial genomes. The median coverage for unfiltered long reads across the nuclear genome was 202X, while median coverage for the chloroplast and mitochondrion were approximately 6201X and 528X, respectively. This was calculated in 1000 base windows using mosdepth (Pedersen & Quinlan, 2018).

A key feature of this updated assembly is the consistency with previous sequencing efforts (Bowler et al., 2008). Previously, 25 centromere sequences were identified (Diner et al., 2017), suggesting that there were fewer than the proposed 33 chromosomes. This agrees with our conclusion of 25 nuclear chromosomes. We independently resolved the location of all the previously proposed partial chromosomes without internal inconsistencies in Fig. 2 (i.e., scaffolds with only one telomere were resolved into a telomere-to-telomere chromosome).

Estimating the number of chromosomes using ultra-long reads

Previous studies have suggested that P. tricornutum has a minimum of 33 chromosomes using pulsed-field gel electrophoresis (Filloramo et al., 2021). Our orthogonal, reference-free method using network graphs of telomere-containing overlapping ultra-long reads revealed 25 chromosomes.

We used two properties of telomeres for this: first, telomeres on linear chromosomes can be identified by unique subtelomeric sequences, and second, that telomere-containing DNA fragments will begin or end with a telomere, representing the start or end of a chromosome. After aligning all telomere-containing reads and retaining only alignments with greater than 95% query coverage, we used iGraph to create network graphs, which resulted in two classes of independent graphs. The first class had 85 independent graphs, each with approximately 50 nodes (i.e., 50 ultra-long reads in each graph), and the second class had eight graphs with approximately 100 nodes (Supplemental Figure S3). In a diploid organism we expect four telomeres per chromosome if we assume that each chromosome has two haplotypes; i.e., (maternal + paternal) ×haplotypes. Under this assumption, 85 independent graphs with approximately 50 nodes represents 21.25 telomeres. Some chromosomes will not have diverged sufficiently, meaning there will be only two telomeres with twice the sequencing coverage per chromosome (maternal + paternal). The remaining eight graphs with 100 nodes each therefore gives a further four chromosomes.

With this logic we estimate 25.25 chromosomes exist in P. tricornutum, which agrees very closely with our final assembly of 25 chromosomes. The additional 0.25 chromosome may be explained by mitotic recombination (Bulankova et al., 2021). Using the features of ultra-long reads at the ends of linear DNA elements (i.e., eukaryotic chromosomes) thus enables an orthogonal method for estimating the number of chromosomes in a reference-free manner.

Assembly quality

To assess the quality of the assembly, we used Merqury (Rhie et al., 2020) to estimate the base-level accuracy and completeness by k-mer frequency, shown in Supplemental File S3. We found that the estimated quality value (estimated log-scaled probability of error for the consensus base calls by Merqury) ranged from 27–53, depending on the chromosome. The mean quality value (QV) for nuclear chromosomes was 28.86, with chromosome 19 as an outlier at 43. The QV for all nuclear genomes except for 19 are likely lower because the chromosomes were polished using heterozygous reads. The chloroplast and mitochondrial genomes have a QV of 53 and 42, respectively. Importantly, the k-mer completeness estimate of 80% suggests that many k-mers in the Illumina reads are not represented in this genome assembly, implying significant haplotype variation. This was also the case when using the Bowler assembly as input for Merqury.

We also estimated the genome completion using BUSCO (Manni et al., 2021). Using the stramenopiles_odb10 model, we found our assembly was 95% complete, with only 3% of expected BUSCOs missing. When evaluating the chromosome scaffolds of the Bowler assembly, we found it was 96% complete with 3% of expected BUSCOs missing.

After removing Lambda spike-in reads with NanoLyse, we found that 98.12% of long reads are recruited by the assembly. When reads are filtered by removing any read that does not align over more than 90% of it’s length (i.e., query coverage is higher than 90%), the number of reads recruited drops to 74%.

Filtered long-read coverage for Chromosome 19 is inconsistent with diploid state

We observed that chromosome 19 has remarkably consistent (i.e., no drops in coverage) filtered long-read coverage relative to the other chromosomes (Fig. 2, Supplemental Figure S2). While we initially predicted P. tricornutum would have two haplotypes since it is diploid, recent work has demonstrated that while each cell has two haplotypes, many haplotypes within a population arise due to mitotic recombination (Bulankova et al., 2021). The consistency of filtered long read coverage for chromosome 19 indicates that there is only a single haplotype, whereas the other chromosomes have two or more haplotypes present, which can be inferred from inconsistent read depth at regions where haplotype divergences occur in Fig. 2 and Supplemental Figure S2. This indicates that there are not two haplotypes for chromosome 19, suggesting a different recent history for this chromosome.

5mC methylation and transposable elements

The Extensive de-novo TE Annotator (EDTA) pipeline (Ou et al., 2019) was used to predict transposable elements in the genome. We found that the majority of transposable elements are long terminal repeat (LTR) retrotransposons (3.43% of the genome was found to be Copia-type, 5.86% were unknown, while terminal inverted repeats were only 1.17% of the genome, and helitrons were 0.54% of the genome). Each LTR region is represented as a shaded blue region in Supplemental Figure S2 in blue, and density plots of the end locations are shown in the top quadrant. Chromosome 19 contained the fewest transposable elements at 50. The locations and density of LTR-retrotransposons are plotted in Fig. 3 for proposed Chromosome 3 and Supplemental Figure S2 for all other chromosomes.

Previous studies have found that some tranposable elements were hypermethylated (Veluchamy et al., 2013). Using chromosome scale nanopore methylation basecalling, we found a strong signal between many predicted LTR retrotransposons and methylation status (Fig. 3, Supplemental Figure S2). To test this, we enumerated all chromosome positions with methylated sites and transposons, and performed a Fisher’s Exact Test, resulting in a p-value of 2.2e−16.

We examined the association between LTR transposon dense regions and regions where the previous assembly failed to generate overlapping regions. We observed that scaffolds with overlapping regions (Supplemental Figure S2) generally were not assembled into full chromosomes because of ambiguity in the placement of the LTR-rich regions at the ends of the scaffolds. These are now resolved by the long-read assembly identified here. Additionally, many of the low-coverage regions of our assembly overlap with the locations of the LTR-dense regions, consistent with chromosomal rearrangements being more likely in these regions. Further investigation at these regions is required.