All LRS technologies fail in the same regions of the D. melanogaster Y Chromosome

Kim et al. (2024) recently published an extensive ONT sequencing data set of many Drosophila species, including a 400× coverage of the reference D. melanogaster strain (iso-1). The D. melanogaster data set was generated using ONT Q20+ (10.4 flow cells), which yields a ∼1% error rate. As DNA was not purified with cesium chloride centrifugation, we expected that the previously observed sequencing bias against Y-linked exons would be absent, namely, that the D. melanogaster Y-linked genes would have uniform coverage (at ∼200×; Kim et al. 2024 used males). However, as exemplified by the kl-3 gene (Fig. 1, rows 1 and 2), this was not the case; indeed, the kl-3 coverage profile in the ONT Q20+ data set was essentially identical to that observed in the earlier PacBio CLR data set (Kim et al. 2014).

Figure 1.

Coverage of the Y-linked gene kl-3 by raw reads across different sequencing technologies. The three LRS data sets exhibit similar coverage profiles, with two regions of the kl-3 coding sequence (CDS) nearly absent in the raw reads. In contrast, these regions are well represented in the Illumina reads. The numbered bars inside the bottom graph mark exon positions. The coverage peak in exon 1 is caused by an exon duplication (release 6 scaffold armY, coordinates 347715–348740) that was collapsed in the Illumina assembly (Supplemental Fig. S1). For a similar analysis of the other Y-linked genes, see Supplemental Figure S2.

Open in new tabDownload PowerPointLink to figure

Furthermore, males of the same strain were sequenced at ∼100× coverage (∼50× for the Y Chromosome) using PacBio HiFi by Shukla et al. (2025), again without cesium chloride purification. Nevertheless, the coverage profile was essentially the same, further reinforcing the conclusion that the observed sequencing bias is not linked to the DNA purification method. As a control using short reads, we sequenced iso1 males with Illumina at high coverage (∼550×) and found that all previously missing exons were faithfully represented in the raw reads, confirming the earlier results obtained by Carvalho et al. (2016) using a smaller Illumina data set. As before, the only major “anomalies” in Illumina coverage are occasional exon duplications (e.g., exon 1 of kl-3) (visible in Fig. 1 as an Illumina coverage peak), which are fairly common in heterochromatic regions (e.g., Tobler et al. 2017; Chang and Larracuente 2019). As expected, the assembly of these strongly biased LRS data sets consistently failed to assemble many exons, whereas the Illumina assembly recovered nearly all exons (Supplemental Fig. S1). All data sets yielded fragmented assemblies as some Y-linked introns are extremely long (in the megabasepair range) and filled with repetitive DNA (Carvalho et al. 2000; Kurek et al. 2000; Reugels et al. 2000), which leads to assembly breaks. As expected, the Illumina assembly is much more fragmented owing to the short read length.

Figure 1 highlights the coverage inconsistency for the kl-3 gene; coverage of Y-linked exons is irregular in most genes (Supplemental Fig. S2), and several exons are strongly underrepresented or entirely absent (missing exons) in multiple Y-linked genes (Fig. 2). On the other hand, a few Y-linked genes (e.g., Pp1-Y1 and Pp1-Y2) display uniform and high coverage in the raw reads, showing that not all Y-linked exons are equally affected. To further assess the coverage uniformity, we looked at a random sample of 10 X-linked genes and found that all have a uniform coverage profile close to the expected 200× (Supplemental Fig. S3). To summarize, most or all X-linked and some Y-linked genes are represented in LRS reads as expected by a random (Poisson) sampling of the genome, whereas most Y-linked genes show highly irregular coverage across different exons, with some exons being nearly absent. Notably, this bias is not exclusive to the ONT data set, as these observations hold for the other LRS data sets mentioned above (PacBio CLR and PacBio HiFi).

Figure 2.

Coverage of selected Y-linked genes in ONT raw reads (Kim et al. 2024). Pp1-Y1 and Pp1-Y2 exhibit uniform coverage, contrasting with other Y-linked genes with highly irregular coverage. The numbered bars at the top of the graphs mark exon positions (Pp1-Y1 and Pp1-Y2 are single-exon genes). Arrows point to the low-coverage exons (blue arrows) and normal-coverage exons (black arrows) that were investigated in detail. For the coverage data on additional Y-linked genes, see Supplemental Figure S2.

Open in new tabDownload PowerPointLink to figure

Finally, even the Pp1-Y1 and Pp1-Y2 genes, which have uniform and high coverage (160× and 100×, respectively), are below the expected coverage of ∼200×. This discrepancy is probably caused by polytenic tissues in adult males, in which heterochromatic regions are known to be severely underreplicated (Yarosh and Spradling 2014). In line with this, Flynn et al. (2020) observed that D. virilis adult whole flies contain 40% less heterochromatic satellites than fully diploid tissues (represented by imaginal discs from larvae), strongly suggesting that the underrepresentation of heterochromatin in polytenic tissues systematically lowers sequencing coverage of these regions in adult flies. This phenomenon differs from the “missing exons” we are investigating here and does not cause assembly problems.

A detailed investigation of the missing exons

Our results show that certain regions of the D. melanogaster Y Chromosome are recalcitrant to LRS. But what exactly is inside these regions? To investigate this, we used the ONT Q20+ data set from Kim et al. (2024) as its huge coverage allowed us to recover a small number of surviving reads from the missing exons. Given the consistency of the bias across platforms, we presume that the same phenomenon is happening with the shallower PacBio CLR and PacBio HiFi data sets (Kim et al. 2014; Shukla et al. 2025). First, we used a simple BLASTN search to recover the few reads covering the missing exons (Fig. 2, blue arrows) and assembled them using a variety of approaches. Commonly used programs such as Canu (Koren et al. 2017) and Flye (Kolmogorov et al. 2019) either failed entirely or yielded very small contigs. The most effective approach we found was to tweak several parameters of miniasm (Li 2016); this yields draft assemblies of all missing exons, with contigs sizes ranging from 12 kb to 102 kb (Supplemental Methods). With these draft assemblies, we could now examine the sequence composition of these regions. As shown in Figure 3, all Y-linked exons, both those with normal coverage and those with low coverage in LRS data sets, are embedded in highly repetitive DNA (∼95% of repetitive DNA). However, there is an important difference between them: Exons that are missing or underrepresented in LRS data sets are closely associated with simple satellite sequences (e.g., (AAGAA)_n, (AATATA)_n), whereas exons with normal coverage are devoid (or nearly so) of satellite DNA in their vicinity. Instead, these regions are embedded in transposable elements (TEs). These consistent association patterns suggest that simple satellite DNA is the major factor disrupting LRS sequencing. For the sake of simplicity, we classified the kl-5 exon 11–13 region as a “normal coverage” region, but it actually is an intermediate case: It has satellite blocks, but they are not in close proximity to the exons (Fig. 3), the raw read statistics (discussed below) are only mildly disturbed, and it has fairly high coverage (Fig. 2).

Figure 3.

Sequence composition of selected regions of the D. melanogaster Y Chromosome, along with a control euchromatic region (white gene). Coding exons are shown as numbered black rectangles (unlabeled when pseudogenes), transposable elements in pink, simple satellites in blue (with darker shades indicating higher AT-richness), and single-copy sequences as a thin black line. The top six regions all have very low read coverage, whereas the bottom five have approximately normal coverage. All regions are represented at approximately the same scale. In the region containing ORY exons 1–2, we could not determine which copy is functional, so we labeled both. The most abundant satellites were identified as described in the Supplemental Methods.

Open in new tabDownload PowerPointLink to figure

How does simple satellite DNA interfere with LRS?

A careful examination of the few surviving reads from the missing exons suggested several potential mechanisms underlying this sequencing bias. First, the base quality values drop precipitously as sequencing enters the satellite blocks (Supplemental Fig. S4). Because ONT sequencing software filters low-quality reads and puts them into a separate file by default, we initially reasoned that the missing Y-linked reads might be found in this low-quality subset. However, after examining these files, we found no enrichment of missing exon reads (Supplemental Table S1), ruling out read quality as the primary explanation for the missing exons.

We also noticed partial exon duplications (i.e., pseudogenes) occurring in tandem with their functional copies. In some cases, these pseudogenes are in reverse orientation in relation to their functional copies (e.g., read SRR26246282.1135981 has two copies of kl-5 exon 14 in reverse orientation). Because single-strand DNA (SSD) forms at some point in all LRS technologies, these reverse-oriented copies could theoretically generate hairpin structures that might interfere with sequencing. However, this could not be a general explanation because most missing exons lack reverse-oriented pseudogenes. Furthermore, in most cases, the reverse-oriented sequences seem to be sequencing artifacts not present in the genome (Supplemental Fig. S5).

Two additional observations proved to be more fruitful. First, we found that reads from the missing exons were noticeably shorter (Fig. 4). Indeed, there is a statistically significant difference in size between them and reads from normal-coverage Y-linked regions (P = 10⁻⁴; nested ANOVA on log-transformed values). This observation strongly suggests that simple satellite sequences disturb the read traversal across the pores (read elongation for short), by either slowing it or causing a premature termination. Regardless of the precise mechanism, this effect alone would reduce the sequencing coverage in the affected regions, making exons near satellites disproportionally underrepresented. We will deal with the second observation in the next section.

Figure 4.

Raw read size in selected regions of the D. melanogaster genome. The six Y-linked regions on the left (blue) have very low read coverages, whereas the next four Y-linked regions and the control euchromatic white gene (black) have approximately normal coverage. The difference in read size between the two groups of Y-linked regions is statistically significant (F_1,8 = 93.2; P = 10⁻⁴), as are the within-group differences (F_8,708 = 8.6; P < 10⁻⁵ ; nested ANOVA on log-transformed values). Note that this figure likely underestimates the association between satellite DNA and smaller read sizes because the low-coverage regions CCY exon 2 and Ppr-Y exon 3 contain a large amount of nonsatellite sequence, whereas the normal-coverage region kl-5 exons 11–13 contain satellite DNA (Fig. 3).

Open in new tabDownload PowerPointLink to figure

Satellite DNA as a barrier to read initiation

Sequencing library preparation is expected to be blind to DNA composition, and hence, all genome regions should be randomly sampled by the reads (Fleischmann et al. 1995). Consider, for example, exon 3 of the kl-5 gene: It should be sampled in both orientations (sense and antisense, or F and R) in equal numbers, namely, no strand bias. Furthermore, the location of the exon within the reads should be random, following a uniform distribution.

Two key findings indicate that the reads covering the missing exons are not a random sample of the genome. First, the exon orientation within the reads (F/R) is strongly skewed in most missing exons, deviating from the 50:50 ratio expected by theory and observed in the four Y-linked regions with normal coverage (Table 1, columns 3–4). Given the small sample sizes for the missing exons, which reduces the statistical power of individual tests, we combined the P-values (column 4) from all six missing exons using Fisher's method (Sokal and Rohlf 1995). This combined analysis rejects the null hypothesis of random exon orientation (χ² = 25.10, 12 d.f., P = 0.014); that is, there is strand bias in these regions. In contrast, applying the same procedure to the four normal-coverage regions yields a nonsignificant P-value (χ² = 3.37, 8 d.f., P = 0.91), confirming that forward/reverse exon orientation follows random expectations in these regions (i.e., no strand bias).

A second sign of nonrandomness in the missing exons reads is that exons seem to be frequently located at the same position within reads instead of being randomly distributed. In other words, different reads seem to start at nearly identical genomic positions. Supplemental Figures S6 and S7 illustrate this “stereotyped” pattern for the kl-3 exon 15–16 and the kl-5 exon 14 regions, respectively, with similar trends observed in most missing exons. We statistically tested this apparent departure from a uniform distribution as follows. First, we used BLASTN to obtain the distance between the start of each read and the target exon. If all reads have the same size and if read starts are random, the null hypothesis for exon positions would be a simple uniform distribution within the range [1, read size]. However, read sizes are variable, and in this case, it seems reasonable to suppose that the proper null hypothesis for n reads would be a composite of n uniform distributions, one for each read size. We confirmed this supposition through simulations, which also validated the statistical test for uniform distribution, described below (Supplemental Fig. S8; Supplemental Table S2). We then compared the observed distribution of exon positions against the null hypothesis of uniform distribution, analyzing F and R reads separately owing to the previously noted strand bias. The result using the Anderson–Darling test, which is more sensitive for small sample sizes than the Kolmogorov–Smirnov test (Razali and Yap 2011), is shown in Figure 5 and Table 1, columns 5–8 (for the results for the Kolmogorov–Smirnov test, see Supplemental Table S3). We found that most “missing exon” regions strongly depart from randomness (i.e., the hypothesis of uniform distribution of exon positions is rejected), whereas in reads from the normal-coverage regions Ppr-Y exon 4, Pp1-Y1, and Pp1-Y2, exon positions are random (i.e., follow a uniform distribution). The kl-5 exon 11–13 region is again an exception, but note that its deviations from randomness are mild (Fig. 5; Supplemental Fig. S9). As we commented before, this region has intermediate characteristics between normal- and low-coverage regions, most likely because it contains satellite DNA but not in close vicinity to the exons (Fig. 3).

Figure 5.

Random genome sampling and read start positions. Under random genome sampling, read start points are expected to follow a uniform distribution (red lines). Exon positions within reads (blue dots) serve as a proxy for read start points in the genome. (Left) The kl-5 exons 3–10, a low-coverage region, show a strong deviation from uniformity (Anderson–Darling test: P < 10⁻³). (Center) Pp1-Y2, a normal-coverage region, shows very good agreement with the expected distribution (P = 0.767). (Right) kl-5 exons 11–13, another normal-coverage region, show intermediate characteristics. Note the fairly good agreement; the P-value of the Anderson–Darling test (P < 10⁻³) reflects the much higher statistical power in normal-coverage regions owing to their much larger number of reads. For the remaining regions see Supplemental Figure S9. Figure 5 provides a formal statistical test for the stereotyped read patterns described in Supplemental Figures S6 and S7.

Open in new tabDownload PowerPointLink to figure

The most likely explanation for this nonrandomness both in F/R exon orientation and exon position within reads is that some genomic regions surrounding missing exons have a much lower probability of serving as successful read initiation sites in ONT sequencing. As a consequence, read starts would be clustered in some genomic regions and depleted in others instead of being evenly distributed. Unless these “forbidden” genomic regions are equally present on both sides of the exon, one of the strands would be preferentially sampled; namely, there would be strand bias.

Given our previous observations (Fig. 3), the most likely culprit for this read initiation suppression is satellite DNA. We tested this hypothesis by counting for each missing exon region how many reads started within satellite blocks and how many started in other sequence types (single-copy, TEs, or the exons themselves). We then compared with a binomial test the observed number of satellite-initiating reads to their expected frequency; the latter is simply the amount of satellite DNA in the region (Fig. 3). As shown in Table 2, there is nearly complete avoidance of satellite DNA as a sequencing starting point. Among the 87 reads covering the missing exons, only one initiated within a satellite block, despite the fact that satellite DNA accounts for 37% to 86% of these regions. This is the main cause of the missing exons: When an exon is located between two large satellite blocks, the only chance of getting sequenced is when reads starts in the small “permissive” regions nearby (TE or single copy), pointing toward it. Reads starting outside the satellite blocks cannot reach the exon because it is too distant, and the satellite DNA cripples the read extension. Ppr-Y exon 3 (Fig. 6) illustrates this pattern very well: It is flanked by permissive TEs on the right side and a nonpermissive satellite block on the left, and reads only initiate in the permissive regions, never within the satellite DNA.

Figure 6.

Sequencing bias in the Ppr-Y exon 3 region. This region has permissive transposable elements on the right side (pink) and a nonpermissive satellite block on the left (blue). The exon itself (shown in red) is flanked on both sides by two short single-copy regions (∼190 bp each; thin black line). A total of 18 reads cover the exon (black arrows), showing a strong strand bias (three forward reads:15 reverse reads) and nonuniformly distributed start points. Fifteen reads (top) originated on the right side, starting at scattered points within the TEs and extending through the exon and partially into the satellite block. Only two reads originated on the left side, both starting in the small 190 bp single-copy region between the satellite block and the exon and extending toward the exon and the TEs. The 18th read (SRR26246282.1589819, shown at the bottom) appears to originate within the satellite block but is actually a chimeric read. It begins with ∼3 kb of Chromosome 3L (a permissive sequence, shown in green) before transitioning into the end of the satellite block, reinforcing the pattern of satellite avoidance as a read start site. This avoidance explains several key observations: the low exon coverage, the strand bias, and the nonrandom exon position within reads (Table 1, columns 3–8), as read start sites are clustered outside the satellite block. We found additional chimeric or rearranged reads in other missing exons, suggesting that the scarcity of surviving reads increases the relative frequency of rare sequencing artifacts (Supplemental Fig. S5).

Open in new tabDownload PowerPointLink to figure

Are there “missing exons” in non-Y-linked D. melanogaster genes?

We addressed this question by screening all D. melanogaster protein-coding genes for evidence of missing exons in two assemblies based on PacBio HiFi reads, which have lower coverage and seem more sensitive to the effect of satellite DNA (Fig. 1). The PacBio HiFi reads were assembled with hifiasm and verkko (Cheng et al. 2021; Rautiainen et al. 2023); we also examined two ONT assemblies generated with Flye (Kolmogorov et al. 2019), with different read length cutoffs (1 kb and 45 kb; the coverages are 400× and 100×). Briefly, we checked for each mRNA, the proportion of its sequence that is present in the assemblies (Supplemental Methods). The rationale for testing multiple assemblies is that autosomal genes have twice the coverage of Y-linked genes, which may mask sequencing biases on the autosomes.

We found that 21 genes (out of 13,986) exhibited missing exons in at least one of the four assemblies (Supplemental Table S4). Among the 21 genes, 10 are Y-linked, 10 are autosomal and heterochromatic (from Chromosomes 2L, 2R, 3R, and 4), and one is autosomal euchromatic. We then examined each gene region for read coverage profiles and the presence satellite blocks, following the same approach used for the Y-linked genes (e.g., Figs. 2, 3). As shown in Supplemental Figures S10–S20, three autosomal genes (including the euchromatic one) were false positives, as they showed high and uniform coverage and lacked satellite blocks in the vicinity of the exons. All three were flagged in the PacBio HiFi verkko assembly, which suggests that this assembler may be less efficient than hifiasm when handling Drosophila highly repetitive regions. The remaining eight autosomal genes (JYalpha, Marf1, DIP-lambda, CG42402, Myo81F, Pzl, Gpa2, and klhl10) had satellite blocks near the exons either within introns or in flanking intergenic regions, which were associated with clear drops in read coverage. Two autosomal genes (Gpa2 and klhl10), along with the Y-linked CG41561, were completely missing from three of the four assemblies, being recovered only in the HiFi hifiasm assembly. This is a striking demonstration of how severe the issue of long-read sequencing bias can be. Taken together, these findings show that in D. melanogaster the close proximity of coding exons to large, simple satellite blocks (and the ensuing sequencing bias) is not restricted to the Y Chromosome.

All cases of sequencing bias we detected so far involve simple satellites, with repeat units of up to 8 bp. We looked for blocks of the complex satellite 1.688 (monomer size, ∼ 359 bp), which is known to occur close to genes (Kuhn et al. 2012), and found that it does not cause a consistent and significant sequencing bias (Supplemental Fig. S21; Supplemental Results).

Non-B DNA and sequencing bias

It is known that satellite DNA and other repetitive sequences can adopt alternative DNA structures (non-B DNA) such as left-handed Z-DNA, three-strand triplexes (H-DNA), four-stranded guanine quadruplexes (G4 DNA), hairpins, etc. (Matos-Rodrigues et al. 2023 and references cited therein). It is also known that non-B DNA interferes with PacBio sequencing, although the known effects are fairly subtle (Weissensteiner et al. 2023). Hence, it is possible that non-B DNA is involved in the missing exons phenomenon. As a preliminary test of this hypothesis, we searched for these structures in the low-coverage and normal-coverage regions (Fig. 3), using the “non-B DNA Motif Search Tool (nBMST) program (Cer et al. 2013). As shown in Figure 7 and Supplemental Figure S22, “direct repeats” and “mirror repeats” are the only motifs that are consistently abundant in low-coverage regions and rare in normal-coverage regions. The potential significance of these findings will be dealt with in the Discussion.

Figure 7.

Abundance of non-B DNA motifs in low- and normal-coverage regions. Motifs were detected using the nBMST program (Cer et al. 2013). The y-axis shows the proportion of each region's sequence occupied by non-B DNA motifs. Blue bars represent low-coverage Y-linked regions; black bars correspond to normal-coverage Y-linked regions and the euchromatic control region (white gene). Note that mirror repeats and direct repeats are consistently abundant in low-coverage regions but are rare in normal-coverage regions. In contrast, guanosine quadruplex motifs, and other motif types shown in Supplemental Figure S22, do not exhibit this pattern.

Open in new tabDownload PowerPointLink to figure

Why do reads seldom start within satellites?

The striking similarity between the PacBio and ONT coverage profiles (Fig. 1) suggests that the same underlying mechanism is responsible for low coverage in both technologies. Despite their fundamental differences—PacBio relies on a DNA polymerase and measures the fluorescence of modified DNA precursors while they are being incorporated into a nascent chain, whereas ONT measures the electrical conductance as native DNA passes through a membrane pore—one shared component stands out: T4 DNA ligase, which is used to glue sequencing adaptors in both platforms.

Evaluation of the possible role of T4 ligase is more complex than it might look at first sight. The enzyme exhibits different behaviors depending on the type of ligation (blunt ligation vs. cohesive-end ligation vs. nick ligation), and PacBio uses blunt-end ligation for gluing the SMRTbell adaptor to the target DNA (https://www.pacb.com/wp-content/uploads/2015/09/Guide-Pacific-Biosciences-Template-Preparation-and-Sequencing.pdf), whereas ONT uses a 1 bp T/A overhang ligation (SQK-LSK114; https://nanoporetech.com/document/genomic-dna-by-ligation-sqk-lsk114). Bauer et al. (2017) reported that the T4 ligase has a preference for AT/TA over GC/CG for blunt end ligation, whereas Bilotti et al. (2022) found that GC is favored in cohesive end ligations. Additionally, ligation efficiency varies significantly even among substrates with the same GC content (Fig. 1 of Bilotti et al. 2022). Given these complexities, it seems difficult to derive from the T4 ligase properties an explanation for the satellite-induced bias we observed. Finally, Jia et al. (2024) developed the LILAP method, which replaced T4 ligase by the Tn5 transposase in the main step of adaptor-target DNA ligation, and tested it on D. melanogaster (iso-1 males), allowing us to examine the effect of T4 ligase on the missing exons. As shown in Supplemental Figure S23, the sequencing bias persists in LILAP reads, albeit to a lesser extent. This reduced bias is expected because of the smaller average size of LILAP HiFi reads (∼5 kb; see comment about Illumina a few paragraphs above). One must also consider that Tn5 has its own sequence biases, which favor GC-rich regions (Kia et al. 2017), and that the LILAP protocol still includes a T4 ligase step to close the nick after the Tn5-mediated transposition. Hence, although DNA ligase remains a possible explanation for the missing exons, the limited available data do not support it as the primary cause. Perhaps the main difficulty of the T4 ligase hypothesis is that all its known preferences are short range (e.g., blunt end ligation of AT/TA vs. GC/CG ends), whereas the effect of satellite DNA on LRS technologies seems to obligatorily involve a much broader scale. For example, several satellites such as (AAAC)_n and (AAAAG)_n, which we found near missing exons, contain AT and GC basepairs and hence should at least partially satisfy T4 ligase requirements, and yet they seldom were used as read initiation points (Table 2). The observation that only highly homogeneous satellite blocks exhibit these effects further suggests that the mechanism at play extends beyond a few basepairs.

Non-B DNA might provide a broad scale mechanism for the effects of satellite DNA on sequencing, and as shown in Figure 7 and Supplemental Figure S22, “direct repeats” and “mirror repeats” are consistently abundant in low-coverage regions and rare in normal-coverage regions. Direct repeats can generate slipped-strand DNA structures and hairpins and are more relevant when the repeat is at least partially self-complementary, as happens in the (CAG)_n and (CTG)_n regions associated with some human diseases (Gacy et al. 1995; Sinden et al. 2007). Although we cannot exclude a role of these structures as a cause of “missing exons,” none of the satellites we found (Fig. 3) have strong self-complementarity. Mirror repeats (or, more precisely, homopurine–homopyrimidine mirror repeats) can form a triple-helix structure (H-DNA) when the DNA is partially denatured (usually by negative supercoiling) and the homopyrimidine single-strand folds back and associates with the duplex DNA via Hoogsteen basepairing (Mirkin et al. 1987; Hisey et al. 2024). Note that many major satellites we found at the missing exon regions are homopurine–homopyrimidine (Fig. 3). It seems that H-DNA can easily forms at mirror repeats during ONT and PacBio sequencing because both technologies produce SSD (the newly synthetized strand in PacBio and the displaced strand that is not traveling through the pore in ONT) that can fold back and associate with the duplex DNA, without supercoiling. It is also likely that a triple-helix DNA forming just ahead of the “active site” of sequencing (the DNA polymerase or the pore) would disturb or hamper the process. Indeed, H-DNA is known to cause stalling of the replication fork (Hisey et al. 2024 and references cited therein). Finally, H-DNA would nicely explain why degenerated satellites seem to be unharmful (Table 3), as a few symmetry-breaking substitutions can abolish H-DNA formation (Mirkin et al. 1987). One way to test the hypothesis that H-DNA causes the “missing exons” problem and, perhaps solve it, would be to add a SSD-binding protein or an SSD-specific nuclease to the sequencing mix. Both enzymes are commercially available, and SSD-binding proteins are indeed used to increase yield and accuracy in PCR templates prone to secondary structures (Kur et al. 2005).

Besides ligase and non-B DNA, there are other possible explanations for the “missing exons.” Satellite DNA might be more resistant to shearing, reducing the number of reads that initiate in these regions. Indeed, Illumina sequencing bias seems to stem largely from DNA fragmentation biases, as sonication or nebulization (two commonly used fragmentation methods) preferentially induces breaks in the middle of CG dinucleotides (Poptsova et al. 2014). Library preparation for LRS is always very gentle, which may facilitate a DNA fragmentation bias. It also occurred to us that most protocols for LRS library preparation use a size-selection buffer that selectively precipitates large DNA fragments; it is quite possible that satellite DNA does not behave as “normal” DNA in this respect and gets lost along with the “missing exons.”

Although we focused more on the characterization and mechanistic aspects of the “missing exons,” there are other interesting (and more directly biological) questions. The strong biases caused by satellite DNA must be related to its poorly known properties in vitro and possibly in vivo. The data shown in Figure 3 allow us to start looking at the sequence level at these mysterious regions of the Drosophila Y Chromosome in which protein-coding exons are embedded in huge blocks of intronic satellite DNA (Kurek et al. 2000; Reugels et al. 2000) and yet are properly transcribed and spliced (Fingerhut et al. 2024). A complete assembly of the Drosophila Y is bound to shed light on some of these mysteries.

The phenomenon of missing exons has obvious relevance for genomics and sequence technology; it most likely is not unique to Drosophila, and it is probably a matter of time before other cases of “missing exons” or unclosable gaps in assemblies are discovered or recognized. We hope that this study stimulates both further investigations into its ultimate cause and improvements in sequence technology and that eventually a T2T assembly of the Drosophila Y Chromosome becomes feasible.

Raw reads

The raw reads used in this work are listed in Supplemental Table S5. All data sets were obtained from adult D. melanogaster males of the reference strain iso-1. For the Illumina data set, we extracted DNA from 40 freshly collected iso-1 males using the Wizard genomic DNA purification kit (Promega A1120), following the manufacturer's recommendations. Library preparation and sequencing were performed at Macrogen (Korea), using the Illumina TruSeq Nano DNA PCR-free library protocol, 151 bp paired-end, with a 350 bp insert size. Assuming a 180 Mbp genome, the raw coverage was 548× (274× for the X and Y Chromosomes); the corresponding values for ONT were 406× (203× for the X and Y Chromosomes) and for PacBio HiF 96× (48× for the X and Y Chromosomes). The ONT sequencing runs were performed in late 2022. The original basecalling was performed using Guppy (version 6) basecaller with the [email protected] (superaccuracy) model. Prior all analysis, reads <1 kb were excluded, and the adaptors were removed using porechop_abi (Bonenfant et al. 2023) with the settings ‐‐ab_initio ‐‐no_split (we did this in part because adaptor sequences would interfere with read start point analysis) (e.g., Table 2).

Genome assemblies

Genome assemblers usually work well around 100× coverage and very high coverages can be detrimental, so we downsampled the ONT and Illumina reads for assembly purposes only. Exon coverages shown, for example, in Figures 1 and 2, were measured with all reads.

Assembly of “missing exons” contigs

Normal-coverage regions of the Y Chromosome were successfully assembled into large contigs using ONT reads (Kim et al. 2024) and Flye (Kolmogorov et al. 2019), as described above. However, all low-coverage regions were absent from this assembly, which was how we initially identified them. To reconstruct these missing exons, we used a targeted assembly approach. First, we used a BLASTN search using the CDS of each “missing exon” as the query and the ONT reads as the database and pulled all matching reads (we found 12 to 18 reads for each region). This procedure reduced assembly complexity by including only reads from the region of interest. We initially attempted to assemble each region separately, using Canu (Koren et al. 2017) and Flye (Kolmogorov et al. 2019), with very poor results (small contigs or no contig at all; it must be added that these tools were not designed for local assembly of highly repetitive regions, under shallow coverage). Using minimap/miniasm (Li 2016) and after trial and error on several parameters (including the undocumented parameter -S, which skips the last steps of miniasm), we eventually succeeded in assembling all “missing exon” regions into contigs ranging from 12 kb to 102 kb. Another helpful procedure was to remove before the assembly the reads that clearly are chimeric or rearranged (e.g., Supplemental Fig. S5). Finally, as miniasm does not have a consensus step, we performed it using racon (Vaser et al. 2017). These steps are detailed in the Supplemental Methods. The polished assembly of the “missing exons” (which should be considered draft assemblies) is available at GitHub (https://github.com/bernardo1963/missing_exons). We also deposited there the sequences containing the normal-coverage exons along with ∼ 50 kb of flanking sequence on each side; these were used in comparisons with the low-coverage regions (e.g., Fig. 3). These sequences were extracted from the Flye whole-genome assemblies mentioned above.

Statistical procedures

Most statistical tests were performed using SYSTAT 13 (nested ANOVA) or custom Python scripts based on the statistics and scipy.stats libraries (scripts are available at GitHub [https://github.com/bernardo1963/missing_exons]). For the Anderson–Darling test, we could not find a program or Python library that allows for a user-specified null hypothesis. To address this, we implemented the Anderson–Darling statistic in Python, allowing for a user-specified null hypothesis. We obtained the P-value corresponding to the Anderson–Darling statistic by calling a modified version of the program AnDarl.c (Marsaglia and Marsaglia 2004). These procedures are implemented in the missingExon_stat_1jan2025.py script, which is available at GitHub (https://github.com/bernardo1963/missing_exons).

Detection of repetitive sequences

We ran Censor (Kohany et al. 2006) locally in order to detect and classify repetitive sequences (TEs and satellite DNA) present in the raw reads and assembled contigs. We found that several simple satellite repeats that are abundant in the “missing exon” contigs (e.g., (AATATAT)_n) were not detected by Censor owing to their absence in its internal reference library (file smprep.ref). We fixed this problem by replacing the smprep.ref file with a complete, nonredundant list of all possible satellites up to 8 bp (file satellite_8_pass2.fasta, available at GitHub [https://github.com/bernardo1963/missing_exons]). Even then, we later found that Censor sometimes misidentify the simple satellites. We wrote a Python script (find_tandem_repeats_v2.py) based on a regular expression that finds any perfect tandem repeat (head-to-tail) present in a DNA sequence (McGinty et al. 2025). This approach is more reliable and was used to detect the most abundant satellites reported in Figure 3 and Table 3 (see Supplemental Methods).

Read coverage estimation

We obtained the read coverage data (e.g., Fig. 1) by doing a BLASTN search of the CDS of the target genes against databases of sequencing reads (ONT, Illumina, etc.). The output was saved in tabular format (m8) and processed using a Python custom script that reports the per base coverage and produces graphical representations of the data (read_coverage_CDS_v6.py; available at GitHub [https://github.com/bernardo1963/missing_exons]). We used WU-blast but obtained essentially the same results when using NCBI BLAST. We have not used BWA (Li and Durbin 2009) or similar read aligner programs because they all assume that the reference sequence contains all sequences present in the reads. This assumption was violated in our case, as we aligned genomic reads against a reference set of Drosophila CDS sequences rather than a complete genome assembly (which is not available for Drosophila).