Dynamic evolution of satellite DNAs drastically differentiates the genomes of Tribolium sibling species

Identification of T. freemani satDNAs

The first step was to characterize the T. freemani satellitome. For assembly-independent identification of satDNAs, we resequenced its genome using Illumina sequencing. The Illumina reads were processed by the TAREAN pipeline, which predicted satDNA candidates from short reads via graph-based clustering (Supplemental Table S1). The TAREAN-predicted sequences were mapped to the T. freemani reference genome assembly Tfree1.0 (Volarić et al. 2022). The candidate sequences that were mapped in five or more consecutive copies were classified as satDNAs. In this way, we identified 135 satDNAs, one already known from previous studies (Juan et al. 1993) and 134 new ones. We named them TfSat01-135 and listed their consensus sequences in the Supplemental Table S2.

Among the 135 satDNAs, we identified 11 that share partial sequence similarity (56.7%–73.0%), suggesting a related origin, so we grouped them into the four superfamilies, SF1–SF4 (Supplemental Fig. S1; Supplemental Table S3). The 135 satDNAs were also searched against the NCBI GenBank database (https://www.ncbi.nlm.nih.gov/genbank/) to find possible similarities with sequences from other organisms. The only similarity found for some satDNAs was with previously described satellites of the sibling species, T. castaneum (Supplemental Table S3). When compared with repetitive sequences from the Repbase collection, the 75 T. freemani satDNAs showed only segmental similarities with different mobile elements, mainly DNA transposons (Supplemental Table S3).

Analysis of the T. freemani satellitome revealed a biased DNA base composition, with 127 of 135 satDNAs showing an A+T content >60% (Fig. 1A; Supplemental Table S3). Regarding monomer length, the satDNAs range from 63 bp to 1106 bp (Supplemental Table S3) but with an evident preference for lengths between 140 and 180 bp (Fig. 1A). Five satDNAs turned out to be higher-order repeats (HORs), composed of two or three 90- to 188-bp subunits that share 63.0%–82.2% pairwise similarity (Supplemental Fig. S2).

Figure 1.

The overview of the T. freemani satellitome comprising 135 satDNAs. (A) The genomic proportions of the three most abundant satDNAs and 132 low-copy-number satDNAs, along with their A+T content and monomer length distribution. The genomic proportions are based on the TAREAN estimates, whereas A+T content and monomer lengths are from individual satDNA data listed in Supplemental Table S3. (B) The presence of satDNA arrays with five or more consecutive copies mapped on the 10 T. freemani chromosomes, fLG2-fLGX, in the Tfree1.0 genome assembly. The number of annotated copies is color-coded according to the color scale. (C) UpSet plot showing the number of the T. freemani satDNAs with annotated genes, transposable elements (TEs), and/or satDNAs within 10-kb flanking regions, based on data in Supplemental Table S6.

Open in new tabDownload PowerPointLink to figure

The 135 satDNAs comprise 38.7% of the T. freemani genome (Fig. 1A). The genome is dominated by a single, extremely abundant satDNA TfSat01, which makes up almost one third of the genome sequence. In addition, the two satDNAs, TfSat02 and TfSat03, are moderately represented (1.6% and 5.7%, respectively), whereas the remaining 132 satDNAs together account for only 0.9% of the genome and can be considered low-copy-number satellites (Fig. 1A; Supplemental Table S3).

Regarding the chromosomal distribution, the repeats of the three most abundant satDNAs were detected on almost all chromosomes, whereas the low-copy-number satDNAs showed different distribution patterns (Fig. 1B; Supplemental Fig. S3). The majority of satDNAs (101) were assigned to only one chromosome in the Tfree1.0 assembly, whereas fewer were mapped to multiple chromosomes (Fig. 1B). To verify the chromosomal distribution of the low-copy-number satellites obtained in silico, we experimentally analyzed six of them with a genome proportion between 0.01% and 0.05%. We found that satDNAs with a genome proportion higher than 0.02% had a broader chromosomal distribution (Supplemental Fig. S3A–C), whereas satDNAs with a lower genome proportion (<0.02%) revealed a distinct signal only on one chromosome pair (Supplemental Fig. S3D–F), which is generally consistent with the in silico analysis. In terms of monomer repetitiveness, 71 satDNAs were found exclusively in arrays of five or more consecutive repeats, whereas the remaining 64 also appeared in shorter arrays with fewer than five consecutive copies (Supplemental Table S4).

To examine whether the arrays of the T. freemani satDNAs are preferentially surrounded by a specific type of sequences, we analyzed their 10-kb flanking regions for the presence of other satDNAs, transposable elements (TEs), and genes annotated in the Tfree1.0 assembly. We found that 111 satDNAs have TEs, 95 have genes, and only 27 of them have other satellite sequences in their 10-kb vicinity (Fig. 1C; Supplemental Tables S5 and S6; Supplemental Fig. S4). TEs and genes are most frequently co-occurring (Fig. 1C), but similar average distances (4–5 kb) between satDNA arrays and the closest annotated element were observed for all three addressed sequence types (genes, satDNAs, TEs) (Supplemental Table S7).

Having obtained a general overview of the T. freemani satellitome, we next focused on the in-depth analyses of the most prominent satDNAs.

TfSat01, the main constituent of the T. freemani satellitome and centromere

We first reexamined the most abundant satDNA. The 166-bp sequence of TfSat01 matches the major satDNA earlier documented to comprise 31% of the genome (Juan et al. 1993), a finding also supported by our TAREAN analysis. Based on 228,060 repeats annotated in the Tfree1.0 assembly, we generated a refined consensus that differs by 5.4% from the older GenBank entry X58539 (Supplemental Fig. S5A). As detailed in the Methods, we consider that our consensus TfSat01 represents the major satellite more accurately, and we used it in all downstream analyses.

The major satDNA was previously reported to be located in centromeric regions (Juan et al. 1993). T. freemani has a typical coleopteran karyotype, 2n = 20, with nine autosome pairs and one pair of sex chromosomes (Shimeld 1989). Females have two Chromosomes X, whereas males carry an X and a small y_p in a parachute-like Xy_p association. To verify TfSat01 presence on all chromosomes, we performed fluorescence in situ hybridization (FISH) on female and male metaphase spreads (Fig. 2A). We indeed detected TfSat01 on all 20 chromosomes in females. In males, however, TfSat01 signals appeared on 18 autosomes and the X, but not on the male-specific y_p, revealing that the sex Chromosome y_p lacks the major satDNA TfSat01 (Fig. 2A).

Figure 2.

Organization of the major satDNA TfSat01 in the T. freemani genome. (A) FISH localization of TfSat01 on female (2n = 18+XX) and male (2n = 18+Xy_p) metaphase chromosomes. The first panels show black and white version of chromosomes to better visualize their contours, especially the male minute Chromosome y_p. Chromosomes are DAPI-stained (blue), and TfSat01 signals are shown in red. A white arrow points to the y_p Chromosome lacking the TfSat01 signal. Scale bar represents 3 µm. (B) Length distribution of TfSat01 arrays in the Tfree1.0 assembly. The median (11.8 kb) is marked by a solid red line, and the mean (14.8 kb) by a dashed line. (C) Box plot showing TfSat01 array length distribution across T. freemani chromosomes (fLGs). Medians are indicated by black lines within the boxes, and the number of arrays analyzed is shown in parentheses below a chromosome name. (D) Schematic illustrating the long-range organization of the multi-megabase-long regions containing TfSat01 arrays. Within TfSat01 arrays (blue rectangles), TfSat01 monomers (blue triangles) are organized into subarrays (black arrows) that differ in monomer orientation. Thus, the term “TfSat01 array” refers to an uninterrupted array of TfSat01 monomers, regardless of how many subarrays it includes. (E) Organization of inversion segments in which TfSat01 monomers change orientation within TfSat01 arrays. The graph shows the length distribution of the inversion segments from the 1793 TfSat01 arrays, with the six most frequent types indicated. The positions of truncated TfSat01 monomers in the most frequent inversion segments are shown in the schematic alignment. Detailed alignments can be found in Supplemental Figure S7A and Supplemental Data S1.

Open in new tabDownload PowerPointLink to figure

One of the main goals of this work was to uncover the long-range organization of the major satellite. The Tfree1.0 assembly is missing the assembled TfSat01 arrays on Chromosomes fLG4, fLG10, and fLGX, but TfSat01-containing multi-megabase regions are successfully assembled on seven other chromosomes (Volarić et al. 2022). Using the StainedGlass tool, we visualized these regions (Supplemental Fig. S6), spanning from 3.8 to 9.7 Mb (Supplemental Table S8). On average, the uninterrupted TfSat01 arrays within the multi-megabase regions have a median length of 11.8 kb (Fig. 2B). Although the number of mapped TfSat01 arrays varies up to 6.4-fold between different chromosomes, the median length of uninterrupted arrays remains rather consistent, ranging from 9.3 to 16.3 kb (Fig. 2C), suggesting a balanced array length across the genome.

The TfSat01 arrays revealed structural complexity extending beyond simple tandem repetition (Fig. 2D,E). We discovered that TfSat01 arrays typically consist of two or more subarrays that show inverted orientation of TfSat01 repeats (Fig. 2D). Arrays with two inversely oriented subarrays are the most numerous (71.8% of all mapped arrays), whereas, among the remaining arrays, those with an even number of inversely oriented subarrays predominate (Supplemental Table S9). The prevalence of an even number of subarrays suggests a tendency toward dyad symmetry, which could lead to the formation of a hairpin structure. We analyzed 1793 two-subarray arrays and found that the hypothetical hairpins exhibit some degree of asymmetry in their stems due to the varying number of TfSat01 monomers forming inversely oriented subarrays. The lack of perfect symmetry indicates flexibility in the organization of the subarrays.

Intrigued by the presence of dyad symmetries, we analyzed the segments where TfSat01 subarrays reverse orientation (Fig. 2E). The majority of these inversion segments fall into one of six conserved types, defined by specific lengths (26 bp, 106 bp, 132 bp, 149 bp, 187 bp, and 200 bp) (Fig. 2E). Each type consists of an abruptly terminated TfSat01 monomer followed by a reverse-oriented monomer truncated at a different position (Fig. 2E; Supplemental Fig. S7A). Despite different truncation sites, all six inversion types share two features: (1) No extraneous intervening sequences separate truncated TfSat01 copies; and (2) the reverse-oriented truncated copies show no increased sequence variability (Supplemental Fig. S7A). Being reverse complementary, truncated copies extend potential hairpin stems almost to the tip, ending in 9- to 27-bp loops (Supplemental Fig. S7B). Notably, all six inversion types are detected on three to seven chromosomes (Supplemental Table S10), and the sequences of each inversion type are highly conserved (95.0%–99.9%), with 100% identical copies present on nonhomologous chromosomes (Supplemental Data S1). These findings indicate that dyad symmetries within the TfSat01 arrays are a widespread and conserved feature of the T. freemani centromeric regions.

TfSat03 resides in the intercalary segments between TfSat01 arrays

Analysis of TfSat01 long-range organization revealed that TfSat01 arrays are separated by ∼4-kb intercalary segments (Fig. 3A,B), which predominantly harbor TfSat03, the second most abundant satDNA (Fig. 1A). We found that these intercalary segments are structured as “cassettes” (Fig. 3A), composed of: (1) a variable number of highly conserved TfSat03 copies, typically between 5 and 8 (Fig. 3C); (2) degenerate TfSat03 repeats that are truncated and/or show <80% similarity to the consensus; and (3) symmetrically arranged 200- to 300-bp outward-facing stretches abundant in (AAT)_n microsatellite (Fig. 3A). It is noteworthy that these components together also form dyad symmetry within the intercalary segments (Fig. 3A).

Figure 3.

Organization and relationships between TfSat03 and TfSat04 satDNAs. (A) Schematic representation of the intercalary segments between TfSat01 arrays, containing conserved and degenerate TfSat03 repeats and (AAT)_n microsatellite arrays. The yellow rectangle marks the central region where degenerate TfSat03 copies change orientation. (B) Length distribution of intercalary segments containing TfSat03 repeats. (C) The number of conserved TfSat03 repeats per segment in the analyzed intercalary segments. (D) FISH colocalization of TfSat01 (green) and TfSat03 (red) on male metaphase chromosomes stained in DAPI (blue). An arrow points to the y_p Chromosome lacking TfSat01 and TfSat03 signals. The bar represents 3 µm. (E) Schematic showing structural differences and pairwise similarities between the TfSat04 monomer and the three subunits of TfSat03. (F) PCA clustering of 159 TfSat04 monomer copies and TfSat03 subunits A, B, and C extracted from 300 randomly selected TfSat03 repeats. The monomers and subunits are represented by color-coded dots. (G) FISH colocalization of TfSat03 (red) and TfSat04 (green) satDNAs on male metaphase chromosomes stained in DAPI (blue). The arrow indicates the y_p Chromosome, which lacks TfSat03 but shows a TfSat04 signal exclusive to this chromosome. The bar represents 3 µm.

Open in new tabDownload PowerPointLink to figure

At the core of the “cassette,” in the midpoint of dyad symmetry, lies a ∼110-bp central region (yellow rectangle in Fig. 3A), where degenerate TfSat03 repeats change orientation, potentially forming a 13-bp loop (Supplemental Fig. S8). The sequence of the central region is so degenerate that the partial similarity to TfSat03 is barely recognizable, but the central region itself is highly conserved across intercalary segments throughout the genome, showing an average pairwise identity of 97.1% (Supplemental Data S2). By undirected graph networks, we visualized relationships between central regions of 1577 intercalary segments (Supplemental Fig. S9; Supplemental Data S3). Whereas some small chromosome-specific clusters were observed, the largest cluster with more than 200 completely identical sequences included intercalary segments from multiple chromosomes (Supplemental Fig. S9). This suggests that the ∼4-kb dyad symmetries of TfSat03-based intercalary segments, along with their preserved central region, are a recurring pattern. In other words, we discovered that the centromeric regions of T. freemani are composed of alternating arrays of TfSat01 and TfSat03 satDNAs, with the arrays of both satellites exhibiting a strong tendency towards macro-dyad symmetries.

To further corroborate the relationship between these two satDNAs in situ, we performed two-color FISH. FISH confirmed colocalization of TfSat01 and TfSat03 in broad regions of 19 T. freemani chromosomes in males (Fig. 3D). Smaller regions showing signals from only one satellite likely correspond to loci where that satellite forms longer uninterrupted arrays. Notably, neither TfSat01 nor TfSat03 signals were detected on the y_p (Fig. 3D, marked by an arrow). The fact that the y_p Chromosome lacks the two most abundant and widespread T. freemani satDNAs raises the question of the y_p satellite profile.

TfSat04, the male sex chromosome-specific satDNA

TfSat03 belongs to the SF1 superfamily, which also includes satellite TfSat04 (Supplemental Table S3), so we next addressed the relationship between the two satDNAs. Unlike TfSat03, TfSat04 is a low-copy-number satellite that comprises only 0.07% of the genome. Beyond the 80-fold difference in genome abundance, the two satDNAs differ in the structure of their repeat units (Fig. 3E). TfSat03 is a 340-bp HOR whose three subunits, A, B, and C, are 113 to 114 bp long and share pairwise similarities of 65.8%–73.0% (Supplemental Fig. S2A), whereas TfSat04 is a 112-bp repeat that corresponds to TfSat03 subunits, sharing 65.8%–72.6% similarity with them (Fig. 3E; Supplemental Fig. S1A). The principal component analysis (PCA) of all 159 TfSat04 monomers mapped in the Tfree1.0 assembly and 300 randomly subsampled TfSat03 subunits separated the sequences into four distinct clusters, reflecting their provenance (Fig. 3F). According to the PCA, TfSat04 monomers are most closely related to TfSat03_subunitA.

We performed FISH to localize TfSat04 on chromosome spreads. Because the Tfree1.0 assembly contains only a single continuous TfSat04 array annotated on Chromosome fLG7, we expected FISH to detect it on one chromosome. Indeed, only one TfSat04 signal was obtained; however, it was unexpectedly located on the y_p Chromosome. Two-color FISH confirmed that TfSat03 and TfSat04 occupy mutually exclusive positions: TfSat03 localized to 19 chromosomes but not to y_p, whereas TfSat04 was exclusively present on y_p (Fig. 3G). FISH on female chromosome spreads, where no TfSat04 signal was detected, confirmed that TfSat04 is specific to the male sex chromosome. The chromosomal localization of TfSat04 thus additionally underscored the diversification that the precursor sequence of the SF1 superfamily has undergone through the formation of TfSat03 and TfSat04 satellites.

The detection of TfSat04 on y_p prompted us, as the authors of the Tfree1.0 assembly, to reexamine the original contigs. After revision (described in the Methods), we concluded that the TfSat04-containing contig, originally assigned to the end of Chromosome fLG7, belongs to the y_p Chromosome. Based on this finding, we will update the assembly to version Tfree1.1, in which the 2.2-Mb end of fLG7 is reassigned as the y_p Chromosome.

Orthologous satDNAs between T. freemani and T. castaneum

After addressing the long-range organization of the most prominent satellites, we focused on the evolutionary trends in the T. freemani satellitome. To this end, we included orthologous satDNAs from its sibling T. castaneum in our study.

The origin of the centromeric major satDNAs of the sibling species

Despite their close relationship, T. freemani and T. castaneum possess entirely different major satellites, TfSat01 and TCAST, respectively. TCAST has been formerly characterized as the primary centromeric component in T. castaneum (Ugarković et al. 1996b; Gržan et al. 2020). To better understand the divergence between TfSat01 and TCAST, we investigated their origin.

By analyzing TfSat01 in the Tfree1.0 assembly, we found that 93.3% of annotated copies show >90% similarity to the consensus (Supplemental Fig. S5B), indicating that TfSat01 is a homogeneous satellite. Nevertheless, there was a distinct subset of repeats showing 82% similarity to the consensus (indicated by the black arrow in Supplemental Fig. S5B). It turned out that these less similar copies were not tandemly organized TfSat01 repeats but parts of the 1106-bp repeat units of the satDNA TfSat02 (Fig. 4A; Supplemental Fig. S10). Notably, TfSat02 repeats, which constitute 1.6% of the genome, are located outside the TfSat01 areas, often far from centromeric regions (Fig. 4B). This opened up the question of the evolutionary relationship between TfSat01 and TfSat02, that is, which satDNA preceded the other: (1) Was TfSat01 an ancestral satDNA whose monomeric unit invaded the TfSat02 precursor sequence? or (2) Was a fragment of the TfSat02 sequence excised and amplified into tandem arrays of the present TfSat01?

Figure 4.

Relationship between the T. freemani satDNAs TfSat01 and TfSat02 and the T. castaneum satDNA TCsat15. (A) Schematic representation of the TfSat01, TfSat02, and TCsat15 repeat units. The blue arrow represents the TfSat01 monomer, and the red and green arrows within the TfSat02 and TCsat15 units, respectively, indicate 166-bp segments homologous to the TfSat01 monomer sequence. The percent sequence similarities are indicated. (B) Genomic distribution of TfSat01 and TfSat02 repeats along Chromosomes fLG2-fLGX in the Tfree1.0 assembly. (C) PCA clustering of TfSat01 monomers and 166-bp segments from TfSat02 and TCsat15. The analysis included 10,000 randomly selected TfSat01 monomers, 2237 TfSat02 segments from T. freemani, and 2788 TCsat15 segments from T. castaneum. Each dot represents a single sequence, and the colors indicate the repeat type.

Open in new tabDownload PowerPointLink to figure

To resolve this, we searched the genome of T. castaneum for related sequences. No tandemized TfSat01-like copies were found, but we identified TCsat15, a previously characterized satDNA in T. castaneum (Gržan et al. 2023), as the ortholog of TfSat02, sharing 85.2% similarity across the entire 1106-bp monomer (Supplemental Fig. S10). Like TfSat02, TCsat15 also contains a 166-bp segment corresponding to the TfSat01 monomer (Fig. 4A; Supplemental Fig. S10), which implies that the TfSat01 sequence was embedded within the ancestral TfSat02/TCsat15 repeat. From this, we conclude that the major satellite TfSat01 is evolutionarily younger and that it emerged from TfSat02. A direct comparison of the TfSat01 consensus with the corresponding 166-bp segments of TfSat02 and TCsat15 showed that the TfSat02 and TCsat15 segments are more similar to each other than to TfSat01 (Fig. 4A). To verify this, we performed the PCA on 10,000 randomly subsampled TfSat01 monomers, all annotated 166-bp segments of TfSat02 in Tfree1.0, and all corresponding 166-bp segments of TCsat15 annotated in the T. castaneum TcasONT assembly. The PCA strongly grouped the TfSat01 monomers into a separate, distant cluster (Fig. 4C). Based on this result and the absence of tandemized TfSat01-like copies in T. castaneum, we conclude that TfSat01 was not present as a canonical satDNA in the common ancestor but rather represents a trait of the T. freemani genome (Fig. 5A).

Figure 5.

Comparison of satDNAs in the centromeric regions of the sibling species T. freemani and T. castaneum. (A) Proposed origin of the major satDNAs TfSat01 and TCAST from an ancestral satellite DNA and a transposon-like element in T. freemani and T. castaneum, respectively. Blue and red boxes show the resulting repeats in each species. (B) Schematic of long-range organization of the species-specific major satDNAs TfSat01 and TCAST, intermingled with orthologous minor satDNAs TfSat03 and Cast7, respectively, in the centromeric regions of T. freemani and T. castaneum. (C) PCA clustering of the orthologous Cast7, TfSat03, and TfSat04 repeats. The analysis included 1000 randomly selected Cast7 monomers from the TcasONT assembly, together with subunits A, B, and C from 1000 randomly selected TfSat03 repeats and all 159 TfSat04 monomers annotated in the Tfree1.0 assembly. Dots represent individual sequences, color-coded by repeat type.

Open in new tabDownload PowerPointLink to figure

To decipher the origin of the T. castaneum major satDNA TCAST, we searched for its presence in T. freemani, but no tandemized TCAST copies were found in the Tfree1.0 assembly. Instead, we identified a TCAST transposon-like element previously reported in T. castaneum (Brajković et al. 2012). This 1093-bp element, resembling a DNA transposon, contains two segments that share 83% similarity with a TCAST monomer and its truncated version in reverse orientation, partially overlapping with ∼290-bp inverted termini (Supplemental Fig. S11A). We found this transposon-like element scattered on all chromosomes in both species, except for the T. castaneum X Chromosome (Supplemental Fig. S11B). Phylogenetic analysis showed that copies from T. castaneum and T. freemani generally separated in species-specific groups (Supplemental Fig. S11C). Considering the 80% sequence similarity and wide genomic distribution in both species, we hypothesize that the TCAST transposon-like element was present in the ancestral genome and gave rise to the 360-bp TCAST satDNA, the current major satellite of T. castaneum (Fig. 5A).

The orthologous low-copy-number satDNAs

In addition to TfSat02, TfSat03, and TfSat04, we identified 11 more T. freemani satDNAs with orthologs among previously characterized T. castaneum satellites (Supplemental Table S3). The alignments of the species’ consensus sequences show pairwise similarities from 63.8% to 94.8% (Supplemental Fig. S13). To explore relationships between monomeric copies of orthologous satellites, we annotated the orthologous repeats in Tfree1.0 and TcasONT and analyzed their sequences (Supplemental Data S4). The PCA plots revealed three clustering patterns: (1) species-specific clustering (TfSat08, TfSat11, TfSat15, TfSat25) (Fig. 6A); (2) extensive mixing of repeats from both species (TfSat05, TfSat22, TfSat111) (Fig. 6B); and (3) segregation of repeats from one species into multiple clusters (TfSat07, TfSat27, TfSat42, TfSat76) (Fig. 6C).

Figure 6.

PCA plots of the T. freemani satDNA (TfSat) repeats and their orthologs among known T. castaneum satDNAs (Cast, TCsat) showing: (A) species-specific clustering, (B) extensive interspecific mixing, (C) segregation of repeats from one species in more than one cluster. Names of satDNA orthologs are indicated above each plot. Dots represent monomer copies from the T. freemani Tfree1.0 (blue) and T. castaneum TcasONT (red) assemblies. PCA is based on alignments in Supplemental Data S4, and interactive plots can be found in Supplemental Data S5.

Open in new tabDownload PowerPointLink to figure

We hypothesized that some of the remaining 120 T. freemani satDNAs may have orthologous copies in T. castaneum not classified as satellites, due to either lack of tandemization or gaps in the Tcas5.2 assembly used to define the T. castaneum satellitome (Gržan et al. 2023). By searching the more complete TcasONT assembly with all T. freemani satDNAs as queries, we identified an additional 45 satellites recognizing related copies in T. castaneum. The PCA of their aligned copies (Supplemental Data S4, S5) showed that they also follow one of the three previously described patterns.

Mechanisms of satDNA propagation

Being among the most variable sequences of eukaryotic genomes, satDNAs are subject to a very dynamic turnover. To understand the mechanisms shaping the T. freemani satellitome, we first focused on TfSat01 and TfSat03. Their intermingling across all chromosomes except the sex Chromosome y_p suggests linked spreading. The conserved inversion points within their arrays, found on nonhomologous chromosomes, further imply a joint mechanism propagating these satellites. One possible mechanism for this collinear spread could be the expansion via extrachromosomal circular DNA molecules (eccDNAs) that contain the TfSat01/TfSat03 arrays and spread them by reintegration into nonhomologous chromosomes. To screen for TfSat01 and TfSat03 in the eccDNA fraction, we performed two-dimensional (2D) agarose gel electrophoresis of total genomic DNA, followed by Southern blotting with specific probes (Fig. 7A). Besides signals from linear chromosomal DNA, both satellites appeared at arc positions, confirming their presence in the eccDNAs. Another potential mechanism for the spread of these satellites across nonhomologous chromosomes involves 3D chromosomal interactions. Indeed, in meiotic prophase I, we occasionally observed the bouquet-like formations of nonhomologous chromosomes associated via (peri)centromeric heterochromatin. Two-color FISH confirmed the colocalization of TfSat01 and TfSat03 at these regions (Fig. 7B). We thus conclude that the concomitance of the two satellites at multi-megabase-long regions might be a consequence of the joint action of both mechanisms—(peri)centromeric 3D interactions in meiotic bouquet-like configurations and reintegration of eccDNAs.

Figure 7.

Experimental evidence of TfSat01 and TfSat03 satDNA repeats in the extrachromosomal circular DNA (eccDNA) molecules and in prophase I associations of nonhomologous chromosomes. (A) 2D gel electrophoresis analysis of the eccDNAs in T. freemani. The schematic (top left) illustrates the migration patterns of linear and circular DNA forms in 2D gel electrophoresis. The ethidium bromide stained gel (top right) with 2D-electrophoretically separated T. freemani genomic DNA was Southern-blotted and hybridized with TfSat01 (bottom left) and TfSat03 (bottom right) probes. Arrows on the blots indicate eccDNAs signals. (B) Meiotic bouquet-like configurations showing heterochromatic associations of T. freemani nonhomologous chromosomes. DAPI-stained chromosomes (top) were analyzed by two-color FISH with TfSat01 (green) and TfSat03 (red) probes. An overlap of signals is shown below. The bar represents 5 µm.

Open in new tabDownload PowerPointLink to figure

Beyond TfSat01 and TfSat03, we observed evidence of expansion mechanisms in other satDNAs. Partial sequence similarities between some low-copy-number satDNAs and TEs (Supplemental Table S3) suggest that TEs contribute to their origin and spread. Orthologous satDNAs proved valuable insights, particularly those linked to the Rehavkus DNA transposons, characterized by inverted termini made of terminal inverted repeats (TIRs) and a variable number of subterminal tandem repeats. In Supplemental Figure S14, we highlight three cases. The first involves the orthologs TfSat15 and TCsat23, both of which form short arrays in the inverted termini of the Rehavkus-1_TC transposon and cluster predominantly in species-specific groups (Fig. 6A; Supplemental Fig. S14A). In the second example (Supplemental Fig. S14B), TfSat23 and its orthologous copies in T. castaneum also evolve in a concerted manner, as does the central part of their associated Rehavkus-3_TC transposon. However, the Rehavkus-3_TC copies in T. freemani typically contain 6–8 TfSat23 repeats, whereas T. castaneum Rehavkus-3_TC copies contain up to 50 repeats in inverted termini. The third and most illustrative example of how TE-integrated repeats can expand into notably long arrays is the TfSat25-TCsat12 orthologous pair (Supplemental Fig. S14C). In T. freemani, TfSat25 forms short arrays of 2–12 repeats within the inverted termini of Rehavkus-1_TC-like elements. In contrast, its T. castaneum ortholog, TCsat12, expanded into long arrays up to 80.5 kb (Supplemental Fig. S14C). Because T. freemani retains the repeats exclusively within Rehavkus-like elements, the expansion in T. castaneum likely occurred after the species diverged, supported also by species-specific PCA clustering of the repeats (Supplemental Fig. S14C).

We also observed concerted evolution in satDNAs not related to TEs. For example, the orthologous satDNAs TfSat08 and TCsat17 are located on Chromosome (f)LG3 in both species and form similarly sized arrays (339 repeats in T. freemani, 290 in T. castaneum). Synteny analysis revealed that these arrays are embedded in conserved genomic regions, with >90% similarity in the 10-kb flanking regions (Supplemental Fig. S15). However, the PCA showed species-specific clustering of repeats (Fig. 6A), suggesting that the ∼40-kb array was present in the common ancestor but diverged through gene conversion after the species split.

Some satDNAs, such as TfSat06, show intense interchromosomal homogenization, indicated by their presence on all T. freemani chromosomes without chromosome-specific clustering (Supplemental Fig. S16A). TfSat06 arrays are organized as long, sometimes >11 kb long, inverted termini of transposon-like elements resembling Rehavkus-1_TC (Supplemental Fig. S16B), suggesting transposition as a likely mechanism for their spread. Because only a single orthologous 20-repeat array was found in T. castaneum, we assume that the expansion and interchromosomal homogenization of TfSat06 in T. freemani occurred after species divergence. In contrast, TfSat22 and Cast9 exemplify a satellite that probably began interchromosomal homogenization in the ancestral genome. This is supported by a strong admixture of their orthologous copies between species, as well as interchromosomal mixing within each species (Supplemental Fig. S17).

Finally, some low-copy satDNAs show intra-array homogenization. For instance, TfSat54 on Chromosome fLG10 forms two short arrays with repeats clustered separately (Supplemental Fig. S18A). Similarly, TfSat97 shows an array-specific clustering of its repeats from two arrays on Chromosome fLG3 (Supplemental Fig. S18B). In both cases, the arrays of these very low-copy-number satDNAs do not exceed 20 tandem repeats, suggesting that they may represent early-stage satellite arrays with the potential for future amplification or loss.

The evolutionary aspects of Tribolium satellitomes

The origin of satDNAs, especially those that dominate genomes, is an intriguing but rarely answered question. Here, we revealed the origin of the genome-dominating and species-specific centromeric major satellites in the siblings T. freemani and T. castaneum. We discovered that the T. freemani major satDNA TfSat01 arose from a segment of another satellite, which is still present in both genomes. Apart from the fact that the primordial satellite still exists, it is equally intriguing why only 15% of its sequence length has evolved into the new satellite. The proliferation of a dominant satDNA from a segment of a repetitive sequence is also demonstrated by the T. castaneum major centromeric satDNA TCAST. Namely, a DNA transposon-like element harboring TCAST was found scattered in T. castaneum (Brajković et al. 2012), but we also detected the same DNA transposon-like element in T. freemani. The previous work on TCAST-like elements in T. castaneum also reported a CR1-3_TCa retrotransposon, which contains a segment of TCAST sequence (Brajković et al. 2012), and the authors pointed out that TCAST might be derived from this retrotransposon. However, we could not detect any CR1-3_TCa copies in T. freemani. Therefore, we believe that it is more likely that TCAST originated from the ancestrally present DNA transposon, whereas the CR1-3_TCa retrotransposon in T. castaneum might have captured part of the TCAST sequence by insertion into the TCAST satDNA array.

The T. freemani and T. castaneum major satDNAs apparently evolved from different types of repetitive sequences, TfSat01 from a satellite and TCAST from a DNA transposon (Fig. 5A). The common link between these source elements is their genome-wide distribution and presence on all autosomal chromosomes. We infer that the ancestral sequence TfSat02/TCsat15, from which TfSat01 arose, must have been widely scattered in the ancestral genome, because in today's genomes both widespread satellites, TfSat02 and TCsat15, often occur as short segments with just one or a few repeats (this work; Gržan et al. 2023). We hypothesize that wide genomic dispersal is, if not a prerequisite, at least a desirable property of sequences that have the potential to be the source of satDNAs. Under circumstances of high DNA turnover, dispersed elements may more easily distribute and exchange sequences with their counterparts from distant regions and/or nonhomologous chromosomes. Importantly, we found no tandemized copies of the T. freemani major satellite in T. castaneum, or vice versa. Therefore, we conclude that the proliferation of the major satDNAs in T. freemani and T. castaneum occurred simultaneously and possibly rapidly. The pronounced amplification of a sequence from the ancestor's genome into a dominant satDNA of the offspring can be explained by the satDNA library theory (Fry and Salser 1977), experimentally supported by different organisms (Meštrović et al. 1998; Silva et al. 2017). However, the satDNA library hypothesis assumes the propagation of the entire pre-existing repeat, whereas the relationship between TfSat01 and TfSat02 in T. freemani shows that it is possible that only a segment of a satellite monomer proliferates into a novel and highly abundant satDNA. This discovery expands the current understanding of the satDNA library concept. If more such examples are found in other organisms, this could eventually help to identify the features that predispose a sequence (or its segment) to proliferate and become the most copious satellite in a genome.

The low-copy-number satDNAs in the T. freemani satellitome show evolutionary statuses from species-specific to orthologous between the two siblings. In both groups, there is a large number of satellites that show partial similarities to TEs. Links between TEs and satDNAs have been observed for some time, and there is growing evidence that TEs are a prolific source of satellite sequences (for review, see Zattera and Bruschi 2022). High-quality genome assemblies of three Drosophila virilis group species showed that much of the divergence in genome composition between the sister species stems from Helitron TE-related tandem repeats (Flynn et al. 2024). Another TE-tandem repeat association in the virilis group involves satDNA repeats, which occur in large inverted termini of the foldback DNA transposon Tetris (Dias et al. 2014). TE-satDNA associations in T. freemani highlighted Rehavkus, the foldback DNA transposons that, like Tetris, carry tandem repeats at their inverted termini. An increased number of tandem repeats in the inverted termini of foldback elements enhances their transposition by strengthening transposase binding and increasing excision frequency (Liu and Wessler 2017). This may explain the benefit for foldback TEs to accommodate satellite repeats. Notably, orthologous satDNAs associated with Rehavkus in Tribolium undergo concerted evolution, but their evolutionary outcomes may differ significantly. Whereas some repeats remain cocooned in the inverted termini of TEs for a long time, some orthologs expanded from TE-residing repeats to satDNA arrays up to 80 kb in length.

TEs likely serve not only as a source of satDNAs but also as vectors for their propagation. Transposition may enable satellite sequences to spread across the genome, including euchromatic regions that were dogmatically thought to be deprived of satDNAs. In T. freemani, we found most low-copy-number satDNAs in the vicinity of genes and TEs, consistent with findings for the low-copy-number satDNAs in T. castaneum (Gržan et al. 2023). In T. castaneum, even more abundant satDNAs have spread into gene-rich regions (Volarić et al. 2024), resembling the euchromatic satDNA dynamics observed in Drosophila melanogaster and its close relatives (Sproul et al. 2020). A growing number of studies report satDNA arrays in euchromatic regions, indicating that euchromatin is not inhospitable to tandemly repetitive sequences (Cabral-de-Mello et al. 2023; Rico-Porras et al. 2024). However, functional studies are required to clarify whether the spread of satDNAs in euchromatin is functionally driven or whether these arrays are just tolerated scattered repetitive “seeds” with the potential to expand.

satDNA-rich centromeric regions in Tribolium

Consistent with the satellite-rich genomes, the centromeric regions of Tribolium abound in satellites. In T. castaneum, satDNA-rich and unusually extended centromeres, so-called metapolycentromeres, comprise almost half the length of individual chromosomes (Gržan et al. 2020), and FISH results in this work reveal a comparable situation in T. freemani. Such abundance raises questions about what enables a sequence to become a centromere-dominant satellite and whether a sequence or a structure is more important. In T. freemani, the long-range organization of (peri)centromeric satDNAs highlighted dyad symmetry as the most conspicuous feature. Dyad symmetries within satellite monomers can form secondary structures such as hairpins or cruciforms, and enrichment in predicted non-B-form DNA structures was found in centromeric repeats of diverse animal and plant species (Kasinathan and Henikoff 2018; Liu et al. 2023). It has been suggested that these structures may facilitate the loading of centromere-specific nucleosomes (Talbert and Henikoff 2022). In T. freemani and T. castaneum, short inverted repeats (<10 bp), potentially capable of forming non-B-form DNA, are present in monomers of major satDNAs (Juan et al. 1993; Ugarković et al. 1996b). In addition to short dyad symmetries, we revealed that T. freemani centromeric regions contain intermingled arrays of major TfSat01 and minor TfSat03 satDNA, both exhibiting striking macro-dyad symmetries based on their several kb long inverted subarrays. The high frequency and preservation of macro-dyad symmetries make one wonder about their biological significance. Because no stable secondary structures formed by multi-kilobase-long inverted satellite arrays have been reported so far, we hypothesize that segments in which satellite subarrays/monomers change orientation and possibly form DNA loops could be essential components of macro-dyad symmetries. The potential loops could contribute to the compaction of centromeric chromatin and/or possibly serve as centromere-specific protein-binding sites. The recent study on human centromeres has shown that the centromeric alpha satellite has an intrinsic tendency to form secondary structures, including hairpins, but also that the binding of the CENPB protein to the satellite repeats promotes formation of submicron-sized DNA loops important for centromere stability and its positioning (Chardon et al. 2022). The conserved inversion points in T. freemani centromeric satDNA arrays raise the possibility that TfSat01 and TfSat03 macro-dyad symmetries might have a similar functional role, although this remains to be tested.

Pronounced dyad symmetries also characterize the centromeric satDNAs of other Tribolium species, which have more complex monomer units. In Tribolium brevicornis, the 1061-bp major satellite comprises two ∼470-bp inverted subunits (Mravinac et al. 2005). In Tribolium madens, the 704-bp MAD2 satellite is also based on dyad symmetry (Ugarković et al. 1996a). Tribolium audax has an even more complex TAUD2 satellite with ∼700-bp-long inverted subrepeats (Mravinac and Plohl 2010). All these Tribolium major satellites span wide centromeric regions of all chromosomes. Therefore, we conclude that dyad symmetry may indeed be a structural imperative and that secondary structures in Tribolium centromeres may be more important than the primary structure. This could explain the lack of sequence similarity between the major centromeric satDNAs of congeneric species. Recent work on two closely related mouse species proposed that rapid evolution of (peri)centromeric DNAs does not obstruct chromatin packaging and chromosome segregation if the satellites adopt DNA shapes recognized by conserved architectural proteins (Dudka et al. 2025).

In mouse centromeres, a low frequency of direction changes in pericentromeric and centromeric satDNAs’ arrays has been attributed to inversion events (Packiaraj and Thakur 2024). In the T. freemani centromeres, however, the direction switches are too pervasive to be explained by inversion events alone. We believe that the frequent macro-dyad symmetries in the T. freemani centromeres may reflect not only a possible structural preference but also underlying mechanisms driving satellite proliferation. One such mechanism could involve eccDNAs, which may serve as a mobile reservoir of repeat sequences that reintegrate into the genome and promote the spread and homogenization of satellites. satDNA-derived eccDNAs have been found in various eukaryotes (Cohen et al. 2003, 2006; Navrátilová et al. 2008; Huang et al. 2021), including T. castaneum (Volarić et al. 2024). In T. freemani, we confirmed the centromeric major and minor satDNAs in the eccDNA fraction, so we suggest that the eccDNAs propagate adjacent arrays of the two satellites, promoting the dyad symmetries present therein. Our conclusion is additionally supported by preserved direction changes in both major and minor satDNA arrays, highly conserved even among nonhomologous chromosomes. In addition to eccDNAs, 3D spatial interactions in the nucleus among distant satellite arrays may facilitate satellite spread between remote loci (Sproul et al. 2020). In T. freemani, we detected clustering of satDNA-rich (peri)centromeric heterochromatin during meiotic prophase, also observed in other Tribolium species (Žinić et al. 2000; Mravinac and Plohl 2010). Nonhomologous centromere pairing during early meiosis has been documented in several organisms (Kurdzo and Dawson 2015), including associations between satDNA-rich centromeric regions of the X Chromosome and different autosomes in mouse (Spangenberg et al. 2021). In mice, high sequence identity of telocentric regions between nonhomologous chromosomes led to the conclusion that frequent recombinational exchange occurs between these chromosomes, likely facilitated by their close association during meiotic prophase (Kalitsis et al. 2006). The clustering of the nonhomologous centromeres of T. freemani in early meiosis and the conserved direction switches between satDNAs arrays at nonhomologous chromosomes could be viewed as supporting this interpretation.

Potential functional implications of dynamic satDNA evolution

From a functional standpoint, an important finding is that T. freemani and T. castaneum, despite different major satDNAs, share an orthologous centromeric organization. In addition to dyad symmetries, in both species, arrays of species-specific major satDNAs intermingle with orthologous minor satDNAs. The unrelated major satDNAs, alongside less abundant but similar minor satellites, challenge the intuitive assumption that the most abundant satDNAs are the critical sequences of functional centromeres. At this point, we hypothesize that the orthologous minor satDNAs (TfSat03 and TfSat04 in T. freemani and Cast7 in T. castaneum) may be centromere-competent. An analogy is drawn with the centromeres of Mus musculus, in which the more abundant major satellite forms the pericentromeric heterochromatin, whereas the less abundant minor satellite is centromeric (Joseph et al. 1989). Interestingly, the mouse Y Chromosome centromere is based on a highly diverged minor satellite-like sequence (Pertile et al. 2009). Likewise, the T. freemani sex Chromosome y_p lacks the major satDNA and instead contains the y_p-specific TfSat04, from the superfamily of minor satDNAs. We speculate that the y_p-specific variant of the minor satellite may have evolved through intrachromosomal homogenization, potentially due to suppressed exchange with other chromosomes, also suggested by the absence of TfSat01/TfSat03 arrays on the y_p Chromosome.

Functional centromeres are usually determined epigenetically by centromere-specific histone H3 variants, CENH3 (Mellone and Fachinetti 2021). In T. castaneum, major satDNA TCAST spans CENH3 metapolycentric domains but also extends beyond them (Gržan et al. 2020). Our preliminary testing of T. castaneum CENH3 antibodies in T. freemani indicates that T. freemani has its own species-specific CENH3, consistent with the centromere paradox hypothesis, which justifies the rapid evolution of centromeric DNA and protein components even among closely related species (Henikoff et al. 2001). Detailed studies of centromere-specific proteins, particularly CENH3, should help to decipher which of the T. freemani satDNAs, major or minor, are actually involved in centromere function.

T. freemani and T. castaneum separated 14 Mya and represent the youngest species within the >100-My-old genus (Hinton 1948; Ramesh et al. 2021). Although it is not known whether the two siblings coexist in the same habitats, they interbreed under laboratory conditions and produce infertile F₁ hybrids (Nakakita et al. 1981), indicating postmating reproductive isolation (Wade and Johnson 1994). Given that the satDNAs are by far the most abundant sequences in their genomes and thus the largest source of genomic difference between the two species, we are tempted to speculate whether they may contribute to reproductive isolation. In two congeneric catfish species, genome-wide satDNA divergence between the parental species has been linked to hybrid sterility (Lisachov et al. 2024). In Drosophila hybrids, differences in satDNA composition between parental species impede the proper clustering of pericentromeric satDNAs into chromocenters in interphase nuclei, thereby promoting hybrid incompatibility (Jagannathan and Yamashita 2021). To determine whether satDNA discrepancies underlie the reproductive isolation between T. freemani and T. castaneum, functional studies similar to those in Drosophila are required, and the hybrids will play a key role in answering this question.

In conclusion, the comprehensive study of the T. freemani satellitome and the comparison with T. castaneum revealed a dynamic evolution of satDNAs that yielded the greatest differences between the siblings’ genomes. Dissimilar levels of conservation, organization, chromosomal locations, and abundances point to different evolutionary dynamics and mechanisms of propagation to which individual satellites within a genome are subject. As for the centromeric satDNAs, we deciphered the origin of the completely different major satellites and uncovered the orthologous organization of the centromeric regions, which show macro-dyad symmetries and related minor satDNAs as a commonality between the siblings. These findings provide a foundation for future work that will address the role of Tribolium’s most prominent satDNAs in the context of functionality and speciation, in which T. freemani-T. castaneum hybrids will be of great assistance. In addition, this work serves as a good reference point for satellitome analyses of other Tribolium species, which should further improve the understanding of satDNA evolution in satDNA-rich genomes. We propose that satellite-abundant nonmodel organisms such as Tribolium can be very useful in tackling questions about functional and evolutionary implications of satDNA behavior, and our further research will go in this direction.

Insect material

The T. freemani beetles, obtained from the USDA Agricultural Research Service in 2015, were reared in whole-wheat flour at 27°C and 50%–70% relative humidity.

Genome assemblies

Repeat annotation and satDNA array analyses were performed using the two high-quality chromosome-level genome assemblies that we recently generated: Tfree1.0 (NCBI GenBank [https://www.ncbi.nlm.nih.gov/genbank/] accession number GCA_939628115.1) and TcasONT (European Nucleotide Archive [ENA; https://www.ebi.ac.uk/ena/browser/home] accession number GCA_950066185.1). The T. freemani assembly Tfree1.0 (262.9 Mb), based on Pacific Biosciences (PacBio) HiFi reads (Volarić et al. 2022), is currently designated as the T. freemani NCBI reference genome. The T. castaneum assembly TcasONT (191 Mb), based on Nanopore long-read sequencing, is the latest T. castaneum assembly, significantly improved in repetitive regions (Volarić et al. 2024).

DNA isolation and whole-genome sequencing

Genomic DNA was extracted from 50 mg of snap-frozen adult T. freemani (10 pooled males and females) using the DNeasy Blood and Tissue kit (Qiagen). DNA quantity and quality were assessed with the Qubit 2.0 DNA HS Assay (Thermo Fisher Scientific) and the Tapestation Genomic DNA Assay (Agilent Technologies), respectively. DNA was sent to the sequencing service provider Admera Health, where library preparation was done using a KAPA Hyper Prep kit (Roche). Illumina 8-nt unique dual-indices were used to mitigate index hopping. Whole-genome sequencing (WGS) on an Illumina NovaSeq X Plus 10B platform (Illumina) generated 2 × 35,516,583 paired-end reads (2 × 151 nt), totaling 10.7 Gb (∼35 × genome coverage of T. freemani).

satDNA mining using graph-based clustering

T. freemani satDNAs were identified from unassembled Illumina reads using graph-based clustering with the TAREAN pipeline (Novák et al. 2017). Raw reads were quality-checked and preprocessed using RepeatExplorer2 tools on the Galaxy server (https://repeatexplorer-elixir.cerit-sc.cz/galaxy/). The interlaced reads were randomly subsampled to reduce the size of the input data set to achieve a low genome coverage required for TAREAN analyses. Six randomly subsampled read sets (175,000–1,500,000 reads; 0.02–0.7 × coverage) were used (Supplemental Table S1). As determined in the first test analysis, a high proportion of TfSat01 satDNA in the genome considerably limits the number of reads that can be processed with TAREAN (Supplemental Table S1, T1 analysis). Therefore, subsequent data sets excluded TfSat01-containing reads (T2–T6 analyses). A custom database was created from consensus sequences of high and low putative satDNAs from all six TAREAN analyses, and duplicates were removed via BLAST. Candidate satDNAs were mapped to the Tfree1.0 assembly, and sequences forming arrays of five or more consecutive monomers were classified as satDNAs. The satDNA consensus sequences were BLAST-searched against the NCBI GenBank (Clark et al. 2016) and screened with the CENSOR tool against GIRI Repbase (Bao et al. 2015).

satDNA bioinformatic analyses

Basic analyses of satellite repeats such as monomer length, A+T content, pairwise similarities, and potential secondary structures were performed using the Geneious Prime 2023.2.1 package (Biomatters Ltd.). All other analyses were conducted as described below, with all R scripts (v4.2.2; R Core Team 2022) available as Supplemental Code.

Revision of fLG7 chromosome sequence in the T. freemani genome assembly Tfree1.0

The T. freemani reference genome Tfree1.0 (GCA_939628115.1) includes nine autosomes (fLG2-fLG10) and the sex Chromosome X (fLGX), but the sex Chromosome y_p remained unassembled (Volarić et al. 2022). Tfree1.0 places the satDNA TfSat04 at the end of Chromosome fLG7, but in this work, we established by FISH that TfSat04 is y_p-specific (Fig. 3G). To resolve this discrepancy, we re-examined the contigs used in the Tfree1.0 assembly (Volarić et al. 2022). We found that the longest continuous TfSat04 array lies on the contig ptg000052l, located distally on Chromosome fLG7 and connected to the rest of the chromosome by an (N)₁₀₀ assembly gap (Supplemental Fig. S19). The ptg000052l contig ends in TfSat04 repeats and has one of the lowest location confidence scores in the assembly (Supplemental Table S3 in Volarić et al. [2022]), so we hold that a slight similarity to TfSat03 caused the contig to be misaligned to Chromosome fLG7. After reviewing its sequence in detail, we conclude that the ptg000052l contig is a candidate for the nonassembled y_p Chromosome and support this conclusion with several findings. First, in addition to the longest arrays of TfSat04, this contig also harbors the longest array of the low-copy-number satDNA TfSat07 (Supplemental Fig. S19), and by FISH, we indeed detected the strongest TfSat07 signal on the y_p Chromosome (Supplemental Fig. S3C). Second, assemblies can be inaccurate in repetitive regions, so we examined the 24 Gb of the raw PacBio HiFi reads used for the Tfree1.0 assembly to analyze gene dosage. Gene dosages for 32 genes annotated on the ptg000052l contig and 2000 randomly selected genes annotated on different autosomes were estimated by mapping the genes at the raw reads with minimap2 (Li 2018). Notably, 32 genes annotated on the ptg000052l contig were represented in the raw reads, on average, two times less frequently compared to 2000 randomly selected genes (Supplemental Table S11). Third, regarding the contig ends, one side of ptg000052l terminates in an ∼18-kb TfSat04 array, and the other in a 1.2-kb (TCAGG)₂₄₃ array (Supplemental Fig. S19), matching the Tribolium telomeric sequence (Mravinac et al. 2011). The repetitive endings suggest that the repeats at both ends extend further, but the ptg000052l contig with its size of 2.2 Mb probably accounts for the majority of the y_p sequence, the smallest chromosome of T. freemani.

satDNA probes

To localize satDNAs on T. freemani chromosomes or to confirm their presence in eccDNA, specific DNA probes were generated. satDNA fragments were PCR-amplified from genomic DNA using primers designed with Primer3 2.3.7 in Geneious Prime and cloned into the pGEM-T Easy vector (Promega). Primer sequences and annealing temperatures are listed in Supplemental Table S12. Recombinant plasmids were transformed into Escherichia coli XL10-Gold Ultracompetent Cells (Agilent Technologies) following the manufacturer's protocol. Positive clones were selected by blue-white color screening and confirmed by Sanger sequencing (Macrogen Europe BV). Probes were PCR-labeled with biotin-16-(5-aminoallyl)-dUTP or aminoallyl-dUTP-Cy3 (Jena Bioscience), with a dUTP:dTTP ratio of 1:2, using specific primers and cloned inserts as a template.

Chromosome preparations, fluorescence in situ hybridization, and confocal microscopy

Chromosome spreads were prepared from freshly dissected T. freemani pupal gonads using the squash method. Chromosome preparations and FISH were performed as described previously (Gržan et al. 2020). Slides were examined with a Leica TCS SP8 X confocal laser scanning microscope with a HC PL APO CS2 63/1.40 oil objective, a 405-nm diode laser, and a supercontinuum excitation laser (Leica Microsystems). Fluorescence signals were captured separately for each fluorochrome and processed using LASX Office 1.4.7 28921 (Leica Microsystems), ImageJ (Schneider et al. 2012), and Adobe Photoshop CS5 (Adobe Systems), using only functions affecting an entire image equally. For each satDNA, a minimum of 10 metaphase spreads from 3 to 10 independent experiments were analyzed.

Extrachromosomal circular DNA isolation, 2D agarose gel electrophoresis, and Southern blot hybridization

To detect satDNAs in extrachromosomal circular DNA, total genomic DNA was extracted from 920 mg of T. freemani larvae and adults as previously described (Volarić et al. 2021). DNA concentration was measured using a Qubit 4 fluorometer and the Qubit dsDNA BR Assay kit (Invitrogen). To eliminate linear double-stranded DNA (dsDNA) from the isolated DNA, 20 µg of DNA were sheared by passing 30 times through a 0.33-mm syringe needle. The sample was then treated with 200 U of Exonuclease V (NEB) at 37°C overnight, which digests linear dsDNA but leaves circular DNA intact. The reaction was stopped with 11 mM EDTA [pH 8.0] and incubation at 70°C for 30 min. DNA was then purified using the Monarch PCR & DNA Cleanup kit (NEB).

For two-dimensional agarose gel electrophoresis, the first dimension was run in a 0.7% ethidium bromide (EtBr)-free agarose gel in 1 × TBE buffer at 0.7 V/cm for 18 h. The gel was then stained in 1 × TBE with 0.3 µg/mL EtBr for 1 h. A gel strip containing DNA of interest was excised, placed at a 90° angle to the direction of electrophoresis on a casting tray, and poured with 1.5% agarose containing 0.3 µg/mL EtBr. The second-dimension electrophoresis was run at 4 V/cm for 3 h in 1 × TBE with 0.3 µg/mL EtBr. Both electrophoreses were conducted at room temperature.

After 2D electrophoresis, the gel was treated with 0.25 M HCl for 30 min to depurinate DNA, followed by 0.4 M NaOH for 30 min to denature it. DNA was transferred overnight onto positively charged nylon membranes (Roche) by capillary blotting in 0.4 M NaOH. Hybridization probes for TfSat01 and TfSat03 were PCR-labeled with biotin-16-(5-aminoallyl)-dUTP (Jena Bioscience) using satDNA-specific clones. Hybridization was performed overnight at 68°C in a buffer containing 250 mM Na₂HPO₄ [pH 7.2], 1 mM EDTA, 20% SDS, 0.5% blocking reagent, and 100 ng of specific satDNA probe. Membranes were washed at 65°C in a buffer containing 20 mM Na₂HPO₄, 1 mM EDTA, and 1% SDS. Detection was carried out using the streptavidin-alkaline phosphatase conjugate and CDP-Star substrate (Roche). The signals were visualized with an Alliance Q9 Mini imaging system (Uvitec).