Long-read genome assembly of the insect model organism Tribolium castaneum reveals spread of satellite DNA in gene-rich regions by recurrent burst events

Chromosome-scale ONT-genome assembly and repetitive element identification

The size of the reference T. castaneum genome assembly (Tcas5.2) is 165.9 Mb (Herndon et al. 2020). Removing placeholders and sequencing gaps results in a genome assembly size of 136 Mb. Considering that the experimentally estimated genome size is 204 Mb (Brown et al. 1990), it is clear that 68 Mb of the genome is missing from the Tcas5.2 reference assembly. To improve the existing assembly, especially in repetitive regions, we employed the ONT long-read sequencing platform. In total, 121 Gb of nanopore long reads were generated. The main steps of the genome assembly process are presented in Supplemental Figure S1. First, nanopore reads were divided into two fractions based on their length: short reads (<20 kb; total of 52.7 Gb; 258× genome coverage) and long reads (>20 kb; total of 36.3 Gb; 178× genome coverage) (Supplemental Table S1). Next, long reads (>20 kb) were used for the initial assembly using Canu (Koren et al. 2017). This initial Canu assembly had 1479 contigs with total length of 321 Mb and a N50 of 835.5 kb (Supplemental Table S2, left Canu statistics). The longest contig was 16.4 Mb in length. Our initial Canu assembly was ∼117 Mb larger than the previously experimentally estimated genome size of 204 Mb (Alvarez-Fuster et al. 1991). This larger genome size likely reflects many sequences that could not be linked to correct chromosomal positions and, as a result, cannot be accurately evaluated. To verify the genome size and genome repetitiveness in silico, the k-mer frequencies of the corrected reads generated during the Canu assembly process using the findGSE and CovEst programs were calculated (Hozza et al. 2015; Sun et al. 2018). findGSE estimated the genome to be 204 Mb, whereas CovEst RE showed only slightly higher estimation of 208 Mb (Supplemental Table S3), which is fully consistent with experimental estimation. The repeat ratio calculated by findGSE was 27% of genome, confirming the previously observed high genome repetitiveness (Wang et al. 2008). Considering that (peri)centromeric satDNA TCAST, with 17% genome abundance, might be problematic for proper assembly (Tørresen et al. 2019), we assumed that the additional 117 Mb in our Canu assembly (Supplemental Table S2) was primarily generated by gene-poor and highly repetitive TCAST arrays. To reduce the suspected negative effects of TCAST satDNA on the genome assembly process, we performed an additional filtering step. In this step, the contigs that did not contain at least 1000 bp of unique gene-coding sequence were removed, resulting in successful filtering of 471 out of 1479 total contigs (Supplemental Fig. S1).

These 471 sorted contigs were further used to achieve a new chromosome-scale assembly using a reference-guided approach. The latest improved T. castaneum genome assembly, Tcas5.2 (Herndon et al. 2020), was used as the reference for organizing these contigs into chromosomes. Two reasons were crucial for the selection of this approach: (1) the availability of the high-quality T. castaneum reference genome, Tcas5.2, and (2) an assembly approach that allows the input sequence to be ordered based on the mapping quality of the reference sequence. Given that the reference genome Tcas5.2 contained 3669 unresolved gaps with total size of 11 Mb and mean gap size of 3125 kb, the first step included gap filling of Tcas5.2 with TGS-GapCloser using 8.6 Gb of Canu-corrected larger reads (length, >30 kb). TGS-GapCloser was efficient in filling most of the gaps, so in total, 3607 gaps (98.3%) were closed, which increased the genome size by 10.5 Mb (Supplemental Table S4). The gap-filled Tcas5.2. reference genome assembly was further used as a template for orientating 471 Canu contigs into the 10 chromosomes using the RagTag software tool (Alonge et al. 2022). A total of 244 contigs were unambiguously assigned to the reference Tcas5.2. genome, resulting in a new ONT-based genome assembly named TcasONT (Supplemental Table S5). The rest (227 contigs) did not accurately map on the assembly, remaining assigned as unassembled sequences (Supplemental Table S6). To further validate structure of omitted contigs, we mapped them with highly abundant (peri)centromeric TCAST satDNAs and found that around half of unplaced sequences (14 Mb) are indeed mostly covered (>50%) with TCAST satDNAs (Supplemental Table S6). To gain insight into the actual representation of (peri)centromeric TCAST, we annotated the TcasONT assembly with TCAST satDNA and obtained 9446 monomers, accounting for 3.6 Mb (1.7%) of the genome.

The final polished TcasONT has a size of 225.9 Mb, with 191 Mb assembled in 10 chromosomes and the remainder in unassembled contigs (Supplemental Table S7). Compared with the reference Tcas5.2, the new TcasONT assembly shows a 45 Mb increase of total chromosome length (Supplemental Table S7). Comparing the TcasONT and Tcas5.2 chromosomes, the elongation of chromosome lengths ranges from 10.6% to 56.3%, indicating a significant improvement of genome assembly continuity at the chromosome level (Supplemental Table S8).

The result of the completeness and quality of the new assembly revealed high levels of macrosynteny and collinearity across all chromosomes with high sequence identity between TcasONT assembly and the reference Tcas5.2 (Fig. 1A). Furthermore, we calculated that 88% of the total Tcas5.2 unplaced contig sequence length is now placed within chromosomes of the new TcasONT assembly.

View larger version:

• In this window
• In a new window

• Download as PowerPoint Slide

Figure 1.

Assessment of T. castaneum assemblies. (A) Dot-plot representation of TcasONT versus Tcas5.2 assembly comparison. The horizontal axis corresponds to the intervals along the TcasONT assembly, and the vertical axis corresponds to the intervals along the Tcas5.2 assembly. LG2–LG10, including LGX, represent chromosomes, whereas unplaced contigs contain sequences in Tcas5.2. that were not associated with any chromosome. Dots closest to the diagonal line indicate colinearity between the two assemblies. (B) Gene completeness assessment based on BUSCO analysis of Tcas5.2 and TcasONT assemblies using insect universal orthologs. The analysis was expressed as absolute numbers for complete and single-copy, complete and duplicated, fragmented, and missing genes. (C) The whole-genome-to-genome dot-plot analyses of the TcasONT (left) and Tcas5.2 (right) assemblies. Each dot represents a region of ≥1000-bp-long region, which is mapped in its entirety to another part of the genome. The dot density is proportional to the number of highly similar regions. (D) Comparison of main classes of TEs (DNA, LINE, LTR, SINE) found by RepeatMasker and tandem repeats (TRs) with monomer lengths of <50 bp, 50–500 bp, and >500 bp found by TRF in the TcasONT and Tcas5.2 assemblies. The height of the bars indicates, from left to right, the total number of TEs or TRs and the cumulative length of TEs or TRs.

We next assessed the gene completeness of the TcasONT assembly by comparative benchmarking single-copy ortholog (BUSCO) analyses between TcasONT and Tcas5.2 using insect universal orthologs from the odb10 database (Fig. 1B). Of the 1367 genes, we identified 1329, which corresponds to 97.2% of BUSCO single-copy completeness. Further, 17 genes were duplicated, 13 were fragmented, and only eight genes were missing. A considerable increase in number of complete BUSCOs (32 genes) can be observed in the TcasONT assembly compared with Tcas5.2. In addition to the BUSCO analysis, we “lifted” genes from the Tcas5.2 annotations onto the TcasONT assembly using the Liftoff package and assessed their final numbers (Supplemental Table S9). Only 48 out of 14,467 genes (0.3%) were left unmapped, mostly representing genes with no known biological functions.

To evaluate the improvement in the representation of the repetitive genome portion in the new TcasONT assembly compared with Tcas5.2, self-similarity dot-plot analyses were performed using >70% sequence identity as a criterion. The comparison of dot plots revealed that the TcasONT assembly has a higher level of self-similarity, with medium to large dark blocks representing regions of repetitive sequences in comparison to the Tcas5.2 assembly (Fig. 1C). Quantified, the number of self-similar sequences was 44,671 in the TcasONT genome assembly, whereas in Tcas5.2, it was only 2202, which means that the new TcasONT assembly has about a 20-fold higher amount of repetitive genome fraction.

To define the improvement specifically in TE content, the TEs in both the Tcas5.2 and TcasONT assemblies were annotated using the RepBase database and RepeatMasker. The TEs were classified into four main types, including DNA elements, long interspersed nuclear elements (LINEs), long terminal repeats (LTRs), and short interspersed nuclear elements (SINEs). The TcasONT assembly revealed an increase of 11,297 DNA TEs, 12,268 LTR TEs, and 27,553 LINE elements (Fig. 1D; Supplemental Table S10) compared with Tcas5.2. Next, we investigated the abundance of these elements in both genomes. Among the four TE categories, in terms of cumulative sequence occupancy, there is a notable increase in TcasONT compared with Tcas5.2 for three of them: DNA TEs (5.2 Mb), LINE (14.5 Mb), and LTR (1.8 Mb) (Fig. 1D; Supplemental Table S10). Taken together, compared with Tcas5.2, the TcasONT assembly is enriched with 21.5 Mb of TEs (Supplemental Table S10).

Additionally, we assessed the content of another class of repetitive DNA sequences, tandem repeats (TRs), detecting them by the Tandem Repeat Finder (TRF), which was successfully used in previous studies of T. castaneum (Pavlek et al. 2015). Detected TRs were classified into three groups based on their monomer length (<50 bp, 50–500 bp, and >500 bp). The analysis revealed a total of 35.3 Mb TRs in the TcasONT assembly compared with 9.1 Mb in Tcas5.2 (Supplemental Table S10), indicating a significant improvement in TR content. On closer inspection, the results show that the number of TR elements detected in the TcasONT assembly was doubled in all three classes (Fig. 1D), whereas the 50–500 bp and >500 bp TRs, which represent “classical” satDNAs, contributed significantly to the cumulative genome length (31.7 Mb in TcasONT vs. 7.3 Mb in Tcas5.2) (Fig. 1D).

From all above, the 45 Mb size difference between the Tcas5.2 and TcasONT genome assemblies is mainly because of the enrichment of repetitive regions (45.8 Mb) (Supplemental Table S10). The repetitive fractions that were most enriched in the TcasONT assembly were TEs (21.5 Mb) and TRs (26.2 Mb). If we consider that the (peri)centromeric TCAST satDNAs account for only 3.6 Mb, it is evident that the TcasONT assembly has an increase of 22.6 Mb in TRs outside of (peri)centromeric regions, which mainly include an increase in “classical” satDNAs with >50-bp-long monomer units.

Characterization of Cast1–Cast9 satDNAs in the TcasONT assembly

Given the demonstrated improvement in the genome assembly, especially in its repetitive regions, the new TcasONT assembly represents an excellent platform for in-depth analyses of previously disclosed abundant satDNA sequences outside the (peri)centromeres, Cast1–Cast9 (Pavlek et al. 2015).

The first step, essential for downstream analyses, was the mapping of the Cast1–Cast9 arrays in the TcasONT assembly. Because of the well-known variability of monomer sequences within a satDNA family, it was important to establish sequence similarity and sequence coverage parameters, which would ensure the detection of the vast majority of Cast1–Cast9 monomers in the TcasONT assembly. For this reason, a BLAST search of the raw nanopore data was performed with the consensus sequences of the Cast1–Cast9 monomers using two parameters: sequence coverage and similarity. The results were shown as density plots in which the color intensity correlated with the number of monomers of a particular Cast satDNA (Supplemental Fig. S2). The majority of Cast satDNAs show the highest intensity, that is, monomer aggregation in the area corresponding to >75% of sequence coverage and >70% of percentage identity. However, the pattern with two areas of high density observed in Cast2 and Cast4 is a consequence of the mutual similarity in one part of the sequences between these satDNAs (Supplemental Fig. S3). Despite their partial similarity, they were still considered as individual satDNAs. Considering density plots for all Cast satDNAs, the parameters sequence coverage >70% and percentage identity >70% were used to ensure mapping of most of the monomers of all nine satDNAs (Cast1–Cast9).

To disclose the organization of the Cast1–Cast9 arrays, we analyzed the relationship between monomer distance in the arrays for each Cast (Supplemental Fig. S4). More specifically, we wanted to determine whether the arrays were organized in tandem continuously (consisting exclusively of monomers of a particular Cast satDNA) or had mixed tandem organization (with a different sequence incorporated into the arrays). The results showed that arrays of most satDNAs show no steep increases in the curve, a pattern representing a typical satDNA array organization consisting of continuously arranged monomeric variants. In contrast, Cast2, Cast5, and Cast7 exhibited a disturbed tandem contiguity demonstrated by sharp increases in average array length at certain distances. The sharpest increase is visible in Cast2, indicating the presence of a specific sequence within the arrays that required further investigation. In addition to monomeric Cast2 arrays, we found Cast2 monomers also organized as a part of a new, longer repeat unit of ∼1270 bp in length (Supplemental Fig. S5A,B). This new family of TRs is termed Cast2′, and in the downstream analysis, these two types of Cast2 arrays were examined separately. A more detailed introspection of Cast5 arrays showed that the observed profile is the result of occurrence of previously described R66-like sequences (Pavlek et al. 2015) dispersed within continuous Cast5 arrays (Supplemental Fig. S5C). In addition, analyses of Cast7 arrays also reveal a mixed organization in which Cast7 monomers often occur in association with (peri)centromeric TCAST satDNAs but with low sequence length and similarity (Supplemental Figs. S4, S5D).

Considering the difficulty of correctly assembling satDNA sequences, the next question was whether the array lengths and arrangement of the Cast1–Cast9 satDNAs in the TcasONT assembly correspond to the real situation of the genomic loci containing the Cast1–Cast9 repeats. Namely, repetitiveness of satDNAs can lead to assembly collapse, which can result in the number of monomer units (array length) in a genome assembly being lower than in the real genome. Given that raw reads provide a realistic picture of what is originally present in the genome as they do not undergo an assembly process, we performed a comparative analysis of Cast1–Cast9 across the individual reads and across the new genome assembly. For this analysis, it was crucial to define optimal parameters for array detection to avoid fragmentation owing to short gaps. Therefore, considering monomer length and mixed array organization of specific Cast satDNAs, such as Cast2, Cast5, and Cast7, the best window length that ensures detection of the maximum number of arrays was evaluated for each Cast satDNA. The analysis showed that a minimum of three consecutive monomer units should be a criterion for correct array characterization for each Cast satDNA. The array patterns between the data sets of individual reads and the TcasONT assembly show remarkable similarity for almost all Cast satDNAs (Supplemental Fig. S6). The only exception is Cast6, whose pattern indicates a suppression of the representation of arrays in the assembly regardless of size. Therefore, Cast6 was not included in the downstream analyses that required a substantial number of assembled arrays.

To quantify the enhanced detection of Cast1–Cast9 regions in the TcasONT genome assembly, we compared both their content and array length with Tcas5.2 (Table 1). In terms of overall enrichment in the Cast satDNA portion, the TcasONT assembly shows an increase of 7.2 Mb in the representation of analyzed Cast satDNAs compared with Tcas5.2. Among them, Cast1, Cast2′, and Cast5 emerged as particularly prominent satDNAs with an increase of 0.7 Mb, 2.8 Mb, and 2.6 Mb, respectively. The total content of Cast satDNAs in the new TcasONT genome was 8.8 Mb (4.6% of the genome), which is remarkably close to the previous experimental estimate of the genome abundance for these satDNAs (>4.3%) (Pavlek et al. 2015). Among them, Cast2′ and Cast5 showed the highest genome abundance with 1.5% for each; Cast8 and Cast9 were the lowest (<0.1%); and the rest remained in the moderate genome abundance level (0.1%–0.4%). The results also show an exceptional increase in the number of arrays in the TcasONT assembly for Cast2′, Cast5, and Cast7. In Cast2′, for example, the number of arrays has increased by more than 9×. In addition to the number of arrays, the average increase in array length (expressed as mean values) can also be observed in the TcasONT genome assembly for all Cast satDNAs. A significant increase in the average array length of up to 3× was observed for the Cast2′, Cast5, and Cast6 satDNAs compared with Tcas5.2. The length of an average array for Cast satDNAs in the TcasONT assembly ranges from 1.2 to 31 kb. The highest average array size of ∼31 kb has been found for Cast6, whereas Cast1 has the longest array of ∼112 kb.

The graph of the Cast1–Cast9 monomer copy number in TcasONT compared with Tcas5.2 also clearly shows the presence of significant enrichment (2–10×) in the majority of satDNAs analyzed (Fig. 2A). From the above, the new TcasONT assembly contains the representative sample of Cast arrays and can be taken as a credible platform for further genome-wide analysis of Cast satDNAs.

View larger version:

• In this window
• In a new window

• Download as PowerPoint Slide

Figure 2.

Distribution evaluation of Cast1–Cast9 satDNAs in TcasONT assembly. (A) Comparison of monomer copy numbers for Cast satDNAs annotated in the reference Tcas5.2 and the TcasONT genome assemblies (see also Supplemental Table S11). The assemblies were mapped with Cast1–Cast9 satDNA monomeric consensus sequences published previously (Pavlek et al. 2015). (B) Comparison of monomer copy numbers for Cast satDNAs annotated in the reference Tcas5.2 and the TcasONT genome assemblies (see also Supplemental Table S11). The assemblies were mapped with Cast1–Cast9 satDNA monomeric consensus sequences published previously (Pavlek et al. 2015). (C) Circular representation of the TcasONT genome assembly with annotated Cast1–Cast9 satDNAs, genes, and transposable elements (TEs). The boxes indicate the most prominent regions where the accumulation of Cast satDNAs correlates positively with gene density and negatively with the TEs density. The position of the (peri)centromere on each chromosome is indicated by a red arrow.

Cast1–Cast9 satDNAs chromosome distribution and genomic environment

The distribution of Cast1–Cast9, including Cast2′ satDNAs on T. castaneum chromosomes, was further investigated by mapping all Cast1–Cast9 monomers on the assembled TcasONT LGs. The two main patterns were observed. The first pattern is characterized by Cast satDNAs located on a subset of chromosomes (Fig. 2B). For example, Cast1 dominates on LG3, LG6, and LG10, whereas it is underrepresented on other chromosomes. Similarly, Cast6 dominates on LG3 but is also present, albeit to a much lesser extent, on LG6, LG8, and LG9. In contrast, the other Cast satDNAs such as Cast2′, Cast5, and Cast7 show a more uniform pattern characterized by similar content on almost all chromosomes (Fig. 2B).

Next, we analyzed the distribution of Cast1–Cast9 along chromosomal length (Supplemental Fig. S7). Because it is already known that TCAST builds up huge blocks of (peri)centromeric heterochromatin (Gržan et al. 2020), the presumed localization of the (peri)centromere on each chromosome is marked. Detailed examination of Cast1–Cast9 distribution disclosed an exceptional spread of Cast satDNAs along the chromosomes with no tendency to cluster in proximal regions of the (peri)centromere. Only Cast7 occurs preferentially in regions enclosed by (peri)centromeric satDNA TCAST. In addition, some Cast satDNAs show a tendency to cluster in distal chromosomal regions. The relatively limited chromosomal distribution observed in Cast6 and Cast8 could be the result of an underrepresentation of Cast6 in the assembly or a recent expansion as shown for Cast8 (see below).

The distribution of Cast1–Cast9 satDNAs on chromosomes compared with genes and TEs was also analyzed (Fig. 2C). Interestingly, Cast1–Cast9 satDNAs are often part of gene-rich regions and rarely overlap with TEs (the most prominent regions with these features are highlighted in Fig. 2C). The only coincidence of the occurrence of Cast satDNA and TEs was found in the (peri)centromeric regions (indicated by red arrows), which tend to accumulate TEs, and also corresponds to the occurrence of Cast7 satDNA.

For a more detailed analysis of Cast1–Cast9 satDNAs arrays, 50 kb flanking regions surrounding Cast1–Cast9 arrays were examined for the presence of genes and TEs. More specifically, all genes in the flanking regions of Cast1–Cast9 arrays were counted, and the base gene content of the whole-genome assembly was calculated and used to assess the relative gene density in the vicinity of Cast arrays (Fig. 3A). From the scaled values based on the total number of arrays, it is evident that Cast satDNAs have median gene content values that are above the median of the genome, except for Cast1, Cast6, and Cast7. Moreover, the arrays Cast2′, Cast3, Cast5, and Cast8 are surrounded by a significantly larger number of genes (Kolmogorov–Smirnov test, P < 0.01) than the rest of the genome (Supplemental Table S12). In addition, these Cast satDNAs have distributions above the third quartile range of the genome, indicating the presence of arrays embedded in extremely gene-rich regions. The only Cast satDNA that has a lower number of adjacent genes is Cast7 (Fig. 3A; Supplemental Table S12). This correlates with the observation that Cast7 arrays occur in an intermingled organization with (peri)centromeric TCAST satDNA (Supplemental Fig. S7). The similar pattern around Cast6 arrays could be the result of the lower total number of Cast6 arrays in the assembled genome. In contrast to the gene-rich regions near the Cast arrays, analyses of Cast satDNAs distribution in relation to TEs showed a significantly lower number of TEs near Cast satDNA arrays compared with the whole genome (Kolmogorov–Smirnov P < 0.01) (Fig. 3B).

View larger version:

• In this window
• In a new window

• Download as PowerPoint Slide

Figure 3.

Distribution of genes and TEs in the vicinity of Cast arrays. (A) The number of genes was calculated based on the distribution of gene content in the assembled genome, with the black line representing the median number of genes in regions of same size (100 kb windows, 5 kb sliding frame) and the red dotted lines representing the first and third quartiles (Q1 and Q3) of the genome-wide distribution. Results of the Kolmogorov–Smirnov tests are represented in Supplemental Table S12. (B) The number of TEs was calculated using the same window size as for the genes, with black and red lines representing the same values. All Cast arrays had significantly (Kolmogorov–Smirnov P < 0.01) fewer TEs in their vicinity compared with the genome-wide distribution, highlighting the fact that they do not overlap with TE presence.

Next, we aimed to determine if there is a correlation between the arrays’ length and the density of nearby genes. For this purpose, Cast1–Cast9 arrays were divided into three classes according to their length (small arrays, <1000 kb; moderate arrays, 1000–10,000 kb; and long arrays, >10,000 kb), and gene density was calculated. The results for each Cast satDNAs are presented in Supplemental Figure S8. Interestingly, most Cast satDNAs have the same profile with the highest gene density being associated with arrays of intermediate length (1000–10,000 bp). In addition, a significant number of genes were observed in the Cast5 arrays, especially near large arrays (>10 kb). In the examined 50 kb windows around satDNA monomers, gene counts show similar trends (Supplemental Fig. S8), without distinct graph patterns, meaning that the arrays are uniformly embedded in gene-rich regions.

Genome dynamics and mechanisms of propagation of Cast satDNAs

We considered that analyzes based on the comparison of the mutual variability of monomers and arrays, together with their position on the chromosomes, can be informative for monitoring the genomic dynamics of Cast satDNAs. Taking into account that the data set (monomers and arrays) was large and the overall variability of monomers within the family relatively low (Supplemental Fig. S2), we noted that phylogenetic analyses of a huge number of similar sequences would not be sufficiently sensitive to disclose trends in Cast satDNAs genome dynamics. Therefore, we decided to perform the principal component analysis (PCA) of monomer alignments and to sequence similarity relationships between the arrays for each of Cast satDNAs (Fig. 4A). First, a genome-wide database of Cast monomers was created and annotated with their actual chromosomal positions. PCA results obtained this way showed a similar trend for almost all Cast satDNAs characterized by a scattered pattern of monomers. This is particularly evident for Cast1, Cast2, Cast2′, Cast3, and Cast4 satDNAs, in which regions of extremely high density of monomers from different chromosomes become visible. The only exception with a pattern characterized by an accumulation of monomers from the same chromosome was observed in Cast6, satDNAs with long homogenized arrays.

View larger version:

• In this window
• In a new window

• Download as PowerPoint Slide

Figure 4.

Evolutionary analysis of satDNA monomers and their junction regions. (A) Principal component analysis (PCA) of aligned Cast1–Cast9 satDNAs monomer distance matrices. Dots represent a monomer unit. Monomers were colored according to the chromosome of origin. The most prominent and recent exchange events between different chromosomes are visible in regions with high density of points with different colors. (B) Heatmap visualization of sequence similarity of 2 kb regions around Cast1–Cast9 arrays throughout the TcasONT assembly. The legend indicates that the sequence similarity of the surrounding regions increases, with colors going from blue to red. (C) Graph networks of Cast3, Cast2, and Cast1 arrays based on their sequence similarity relationship. Cast3 shows a cluster containing all arrays. Cast2 shows a pattern with six different clusters containing arrays mixed from different chromosomes. Cast1 shows no significant clustering of arrays, except small divergent clusters. (D) Sequence logo of junction regions of Cast1, Cast3, and Cast9 arrays in which an enrichment of poly(T) and poly(A) tracts can be recognized.

Second, after a database of Cast arrays with corresponding chromosome annotation was created, comparative analyses of the arrays were performed based on their mutual sequence variability for each Cast satDNA. Sequence similarity relationships between the arrays for each of Cast satDNAs are represented by graph networks (Fig. 4C; Supplemental Data, Graph Networks), except for Cast6 and Cast7, which were not analyzed because of the small number of arrays. The distance and interconnectivity between dots correlate with the sequence similarity between arrays, with closer dots indicating higher similarity of monomers between different arrays. Our assumption is that the relationship between the arrays reflects their spread throughout genome. Three graphs, Cast1, Cast2, and Cast3, were selected as illustrative examples (Fig. 4C), whereas the other Cast satDNAs arrays networks are presented in the Supplemental Data, Graph Networks. The Cast3 shows the pattern of array dynamics with a dominant cluster containing relatively closely related arrays from all or almost all chromosomes, and the similar pattern applies to Cast4, Cast8, and Cast9. In addition to the dominant mixed clusters, Cast4, Cast8, and Cast9 also have some small distant clusters (up to 10 arrays) from the same chromosome. The second pattern includes Cast2 with six relatively distant array clusters with related arrays from different chromosomes, indicating multiple genome events of intense array expansion during Cast2 evolution. Interestingly, only one of them represents intrachromosomal expansion, whereas the rest show very extensive interchromosomal exchange. A similar pattern with several distant clusters containing related sequences from different chromosomes is also observed in Cast2′ and Cast5. Finally, the Cast1 network shows clusters with higher distance between arrays, with some array sets completely separated owing to their sequence divergence. The Cast1 network also shows two separate subgroups of sequences representing different variants, one of them (Fig. 4C, labeled) directly connected to the TE, Polytron (Supplemental Fig. S9).

Interestingly, these three different patterns of Cast satDNA propagation events correlate with the average lengths of the arrays. For example, satDNAs for which only one expansion event can be observed (Cast3, Cast4, Cast8, and Cast9) have a relatively short array length (mostly around 4000 bp). SatDNAs with several expansion events, as seen in Cast2, Cast2′, and Cast5, have a moderate array length of ∼15,000 bp. Finally, Cast1, for which no recent expansion centers were observed, intends to have very long arrays (up to 112 kb).

To gain insight into the spread of Cast1–Cast9 arrays throughout the genome, the surrounding regions of the Cast arrays were also analyzed. To address this, regions of 2 kb in length upstream of and downstream from individual Cast arrays were extracted. The extracted regions were aligned using the MAFFT program, and the sequence similarity values were calculated from the alignment for each Cast satDNA. The results showed that of the 10 Cast satDNAs analyzed, only Cast5 and Cast7 show similarity of surrounding regions for most arrays (Fig. 4B). A detailed analysis of the surrounding regions for Cast5 arrays revealed two dominant regions consisting of sequences with high similarity. One of them is the R66-like sequence, found as a scattered sequence in the Cast5 arrays (see Supplemental Fig. S5C) being also located at one end. The R66 element overlaps with Cast5 monomers, particularly with divergent Cast5 monomers. On the other side of most Cast5 arrays, the R140-like sequence was found. Detailed analysis showed that out of 150 Cast5 arrays, two-thirds of the arrays had R140-like sequences, and one-third of the arrays had R66-like sequences at the ends of the arrays. Likewise, Cast7 arrays showed that most of them are surrounded by the (peri)centromeric TCAST (Supplemental Fig. S5D). The rest of Cast satDNAs showed only partial similarities of the surrounding regions noticed only for a smaller number of arrays (e.g., Fig. 4B, Cast1 panel). For example, a high-similarity subset of Cast1 arrays showed that those nine arrays originate from the same chromosome (LG7), having the same transposon-like sequence at the end of arrays (Supplemental Fig. S9).

For proper examination of the precise Cast array junction regions, their edges need to be precisely defined. Because the monomers at the ends of the arrays have a higher variability owing to the less efficient recombination at the ends of arrays (for review, see Lower et al. 2018), we developed a k-mer similarity-based approach. The results of the analyses did not reveal any part of the monomer that would be favored at the ends of the Cast arrays (Supplemental Data, Breakpoint alignment). However, the analyses of 20 bp sequence motifs near the array boundaries showed that Cast1, Cast3, and Cast9 have regions in the vicinity of arrays with significant poly(A)/poly(T) tracts, whereas the other Cast satDNAs do not share any common motif (Fig. 4D; Supplemental Fig. S10).

Considering the observed extensive inter- and intrachromosomal exchanges of all Cast satDNAs, one of the possible mechanisms that influence this phenomenon could be the insertion of satDNA arrays mediated by extrachromosomal circular DNA (eccDNA). To test the presence of Cast satDNAs in the eccDNA fraction, we performed two-dimensional (2D) agarose gel electrophoresis followed by Southern blot hybridization. The hybridization probes were developed for the most abundant Cast satDNAs, Cast1, Cast2, Cast5, and Cast6 (0.3%–1.5% of the genome), and the eccDNA molecules containing the examined satDNAs were indeed detected (Fig. 5).

View larger version:

• In this window
• In a new window

• Download as PowerPoint Slide

Figure 5.

Detection of extrachromosomal circular DNA (eccDNA) molecules containing Cast1, Cast2, Cast5, and Cast6 satDNA. Results of T. castaneum genomic DNA separation in 2D electrophoresis are shown on the left, and Southern blot hybridization with specific probes for each Cast element is shown on the right. Signals obtained by Southern blots confirm the presence of specific Cast satDNA in chromosomal linear double-stranded DNA (linDNA) but also in eccDNA molecules.

In general, chromosomal regions with suppressed or without recombination tend to accumulate satDNAs (Stephan 1986). Because sex chromosomes represent regions with suppressed recombination, we examined the number and length of Cast1–Cast9 arrays on the X Chromosome to compared it with autosomal chromosomes. This type of analysis would be even more informative using the Y Chromosome, which is mostly nonrecombining, but unfortunately, the Y Chromosome was not available in the former Tcas5.2 nor in the new TcasONT assembly. Mapping of Cast1–Cast9 arrays to chromosomes disclosed that the X Chromosome did not have a significantly higher average number of Cast1–Cast9 arrays per megabase, compared with the others (Supplemental Table S13). The comparison of the sequence variability of Cast satDNA monomers from X Chromosome in contrast to autosomes did not show significant difference. However, the analysis of Cast1–Cast9 array length showed a statistically significant longer length of arrays on the X Chromosome compared with the others (Wilcoxon test, P < 0.05) (Fig. 6). This trend is especially prominent in Cast 2, Cast2′, Cast5, and Cast9, for which the average length of the arrays is increased even 10× at the X Chromosome (e.g., Cast2′) (Supplemental Fig. S11).

View larger version:

• In this window
• In a new window

• Download as PowerPoint Slide

Figure 6.

Distribution of Cast1–Cast9 satDNA array lengths across chromosomes, LG2–LG10 including the X Chromosome (LGX).