Splicing deficiency may represent a critical yet underexplored form of genomic erosion in non-recombining regions. Across four phytoplankton species diverged ~333–639 million years ago, genes within U (female) and V (male) “UV” mating-type regions—non-recombining chromosomal regions that determine mating compatibility—show strikingly elevated intron retention relative to genes in other genomic regions. Long-read data reveal abundant aberrant, likely non-functional mRNA isoforms despite preserved coding potential. This preservation suggests that splicing defects arose early in UV evolution and have persisted over deep time. We propose that these defects arise from evolutionary changes in sequence composition and chromatin organization that accompany recombination suppression, such as reduced GC content, altered nucleosome occupancy, and disrupted methylation, that collectively compromise splicing fidelity. Unlike sex chromosomes, which often degenerate through gene loss, splicing-deficient UV regions in green algae retain hundreds of genes, indicating that transcript-level dysfunction provides an alternative route to functional decay. Our results identify chromatin-mediated splicing deficiency as a novel axis of genomic erosion and position algal UV systems as models for studying how recombination suppression reshapes RNA processing fidelity in essential, non-recombining genomes.

Splicing in green algae mating-type region is highly aberrant

Transcriptome analysis among four Mamiellales species revealed widespread and apparently conserved splicing misregulation within the mating-type region (Fig 1). We compared splicing patterns between previously identified mating-type regions to other regions in four Mamiellales species [6]. In Micromonas pusilla the mean intron retention frequency outside of the mating-type region was 0.011, while within the mating-type region the mean frequency was 0.33 (p < 0.0001, Mann–Whitney U). We estimate that nearly half (45%) of all transcribed introns derived from genes within the mating-type region are retained in these data. This difference (33% versus 45%) indicates that the frequency of intron retention is more common in transcripts of highly expressed genes (see also below). Despite approximately 339–630 million years since the last common ancestor of these taxa [15], this pattern is qualitatively consistent across each species (Fig 1 and S1 Table), implying that the processes that contribute to this apparent aberrant splicing have been operating since shortly after the emergence of the UV mating-type region and have been retained in these diverging lineages.

Possible mechanistic basis of aberrant splicing

Introns in mating-type genes are highly divergent from those elsewhere in the genome. They are significantly shorter with lower GC content (S1 Fig), both features that affect splicing efficiency [16,17]. Most strikingly, branchpoint sequences which are critical for spliceosome assembly, are significantly depleted in all species within the mating-type region (p < 0.0001 for all species; S3 Table). Splice-site flanking sequences also reflect the broader GC depletion in the mating-type region, though specific positions within the extended splice site consensus motifs retain their nucleotide compositions (Figs 2 and S2). Splice-site flanking sequences reflect the broader GC depletion in the mating-type region, though specific positions within the extended splice site consensus motifs retain their initial nucleotide compositions (Figs 2 and S2). For example, guanine at the fifth intron position of the 5′ splice site—which contributes to splicing efficiency by base pairing to U6 snRNA [18–20]—is significantly less frequent in mating-type introns (8%−23% reduction; p < 0.0001 for all species, Fisher’s exact test; S2 Table), yet these sites remain enriched for guanine relative to other intronic positions, suggesting that natural selection is relaxed but detectably present. In contrast, canonical splice sites (5′ GT/GC and 3′ AG) occur at similar frequencies inside versus outside the mating-type region (S2 Fig), with only two subtle exceptions: M. commoda shows slightly fewer canonical 3′ splice sites in mating-type genes (96.1% versus 98.6%, p < 0.0001), while M. pusilla shows slightly more canonical 5′ splice sites (96.6% versus 98.2%, p = 0.41E-3‌‌). These results indicate that while mating-type introns differ markedly in composition and length, natural selection maintains core splicing signals but acts less efficiently on secondary elements that affect splicing efficiency.

Random forest modeling considering only introns within the mating-type loci confirmed that shared transcript and genomic features predict the occurrence of splicing dysfunction across all four species (AUC-ROC: 0.735–0.867). Transcripts per million, a measure of gene expression, was the strongest predictor, wherein highly expressed genes showed the highest intron retention rates across all species. Shorter introns, lower GC content, and greater distance between the intron and the transcription start sites also predicted retention (Fig 2 and S4 Table), while potential functional consequences (frameshifts, stop codons) and the presence of specific splicing features (branchpoints and splice-site adjacent nucleotides) showed little predictive power and were rarely retained in final models after applying minimal feature selection techniques. In contrast, Random Forest models trained on autosomal introns failed to achieve meaningful predictive performance (F1 = 0.00–0.17) due to extreme class imbalance where retained introns comprise only 0.6%–2.6% of autosomal introns (S4 Table). The absence of predictive power for frameshifts or stop codons may reflect active nonsense-mediated decay (NMD). Core NMD machinery (UPF1, UPF2) [21] is expressed across all species (S7 Table) and transcripts with retention induced frameshifts or stop codons would be primary degradation targets. Our measurements therefore likely capture the subset of retained introns that saturate or escape NMD. Notably, 16%–35% of variance in intron retention occurred at the gene level across all species, indicating coordinated splicing regulation within genes rather than independent intron-by-intron effects. Consistent with this coordination, intron number per gene was strongly correlated with retention rate in mating-type genes (Spearman’s ρ = 0.56–0.65, p < 0.0001) but not in autosomal genes (ρ = 0.06–0.22), suggesting that splicing defects compound across multiple introns specifically in the mating-type region.

The mating-type region of M. pusilla is substantially different from the other genomic regions with respect to the local chromatin environment [8]. MNase-based measurement of chromatin accessibility [22] suggests that genes residing in the mating-type region have substantially altered nucleosome occupancy. Relatedly, there are fewer CpG dinucleotide sites in exonic regions flanking each intron within the mating-type region consistent with lower overall GC nucleotide content (S3 Fig). Nonetheless, CG dinucleotides that are present within the mating-type region are significantly less likely to be methylated on a per-site level than are those in the autosomal regions of the genome (S3 Fig). This implies that overall GC content and methylation efficiency are critical determinants of differences in methylation patterns.

Diverse isoforms are created by widespread missplicing

Additional categories of alternative splicing are also enriched in genes in the mating-type region beyond intron retention, consistent with a broad disruption of splicing fidelity. We circumvented the vagaries of interpreting short-read RNA sequencing data by applying a long-read circularized consensus method [23] to generate 14,350 high-quality full-length isoforms for the M. pusilla reference strain. We quantified additional classes of alternative splicing events using junction coordinates extracted directly from R2C2 isoform consensus sequences. For each gene with at least two isoforms, we identified 5′ splice sites (A5SS), alternative 3′ splice sites (A3SS), exon skipping events (SE), and retained introns (RI) based on comparisons of intron–exon boundaries across isoforms of the same gene (see Methods). Because these event types are defined exclusively by splice junctions, their detection cannot be inflated by intron retention. Classical alternative splicing events (A5SS, A3SS, and SE) are rare in M. pusilla: combined, they account for approximately 100 events among mating-type genes and 190 events across 2,458 multi-isoform autosomal genes, one to two orders of magnitude fewer than retained introns (629 in the mating-type region; 1,595 autosomal). Critically, every alternative splicing event category is enriched on a per-gene basis in the mating-type region relative to autosomes (A5SS 8.2×, A3SS 6.9×, SE 2.3×, RI 2.8×; p < 0.0001 for each; S8 Table). These results indicate that the mechanism driving elevated splicing variability in the mating-type region operates broadly across event types and is not limited to intron retention.

Analysis of full length isoforms in M. pusilla produced detailed insights into mRNA processing at the transcript level. Using a linear regression model, isoform number increases with both expression and intron count (0.27 isoforms per log-expression unit, 0.76 per intron). Genes in the mating-type region produce ~1.9 more isoforms than autosomal genes, with predicted means of 1.83 (non-mating-type region genes) versus 3.72 isoforms for mating-type region genes at average expression and intron number (p < 0.0001, OLS; Fig 3). Additionally, genes within the mating-type region show more evenly distributed isoform expression, with less dominance by a single major isoform, compared to autosomal genes (p < 0.0001, GLM; see Methods). This pattern suggests that purifying selection is less efficient at suppressing non-canonical isoforms in the mating-type region.

We also compared other features of transcript diversity including transcript start sites and polyadenylation sites for genes within and outside of the mating-type regions. This analysis revealed selective disruption of 3′ end regulation whereby genes in the mating-type region exhibit significantly more distinct transcript termination sites compared to genes outside the region (39.4% increase; p < 0.0001, Poisson GLM). This indicates that like splicing, termination and polyadenylation may be similarly disrupted. However, genes within and outside the mating-type region display comparable diversity in transcription start sites across isoforms (p = 0.563, GLM). We speculate that while the accumulation during evolution of AT-rich sequences may contribute to increased alternative polyadenylation, mispriming of oligo-dT primers during library construction could artifactually contribute to this pattern (although see below). Together, these findings suggest that co-transcriptional and post-transcriptional RNA processing—but not transcription initiation—are perturbed in the recombination-suppressed mating-type region.

Transcripts derived from genes within the mating-type region are less likely to be functional than those in other regions of the genome. In the absence of direct experimental validation, it is challenging to determine if a given transcript produces a functional protein. We therefore indirectly estimated functionality by asking what fraction of full-length transcripts are derived from isoforms that contain at least 80% of the annotated CDS in the correct reading frame without premature termination. Strikingly, transcripts from genes within the mating-type region showed a significantly lower fraction meeting this threshold (p < 0.0001, GLM, Fig 3B). To further assess the effect of alternative isoforms on protein function, we annotated functional protein domains and ontologies for each isoform. Isoforms derived from genes in the mating-type region also more frequently exhibit functional domain disruptions and changes in functional ontology compared to genes in the rest of the genome (S5 Table). For example, for genes in the mating-type region, ~87% of isoforms contained at least 90% of annotated functional domains on average, compared to >92% for genes in autosomal regions (p < 0.001, Mann-Whitney U; S5 Table). These patterns suggest that the elevated alternative splicing observed in the mating-type region is primarily deleterious, though we cannot exclude the possibility that some alternatively-spliced isoforms are neutral or positively selected.

Mispriming during cDNA preparation is unlikely to explain increasing missplicing in the mating-type region

A potential confounder for intron retention estimates in AT-rich regions is the capture of immature transcripts through oligo-dT mispriming, but several lines of evidence argue against this explanation. First, R2C2 isoforms were processed with Mandalorion v4.6.0 using an A-content filter (Acutoff = 0.5) that discards any read whose 20 bp downstream genomic window exceeds 50% adenine. Second, the GC content of the M. pusilla mating-type region (~42%) is comparable to that of the human genome, where oligo-dT mispriming is not typically considered a dominant artifact. Third, as described above, the alternative-splicing events beyond intron retention are enriched in the mating-type region, consistent with widespread misprocessing. Finally, to directly test whether downstream A-content predicts intron retention, we fit a logistic regression model within the 164 mating-type genes that contain both IR and non-IR isoforms (n = 1,393 isoforms). The per-isoform percentage of adenine in the 20 bp window downstream of the transcript termination site showed no significant association with intron retention (p = 0.31; S5 Fig). Together, these results demonstrate that internal priming artifacts of immature transcripts are unlikely to account for the widespread splicing deficiency we observe.

Fig 4. Proposed model of genomic erosion in the Mammiellales mating-type region.

Schematic comparing autosomal regions (left) with mating-type regions (right), illustrating how recombination suppression initiates a cascade of molecular changes. Suppression of recombination eliminates GC-biased gene conversion (gBGC), allowing GC→AT mutational bias to reduce GC content from ~60% to ~40%. This compositional shift reduces CpG dinucleotide availability and methylation efficiency, disrupts nucleosome positioning (elevated MNase coverage indicates more open chromatin), and shortens introns through deletion bias. These chromatin and sequence changes collectively impair splicing fidelity: ~35% of mating-type introns are retained (thick orange bars in gene models) compared to ~3% in autosomal regions, producing abundant truncated, dysfunctional transcripts despite preservation of coding sequences. Expression levels (number of gene boxes) are elevated in the mating-type region due to altered chromatin accessibility, further exacerbating splicing defects through kinetic mismatches between transcription and co-transcriptional splicing.

https://doi.org/10.1371/journal.pbio.3003823.g004

Software development

We used LLM support to develop initial scripts and analyses. Each script was tested and ultimately validated by the authors prior to inclusion in this manuscript.

Comparative intron retention analyses

We identified publicly available short-read RNAseq data for four species in the Mamiellales family through metadata search on the sequence read archive (S1 Table). We further restricted the search to species whose genomes are chromosome-level assemblies with a mating-type region defined in [6]. We aligned RNAseq data to each genome using STAR [36] version 2.7.10b and provided the appropriate gtf file during index construction, and the “--quantMode GeneCounts” option to obtain estimates of gene expression, but used otherwise default parameters. We used rMATs [37] version 4.3.0 to quantify intron retention providing the parameters “--allow-clipping --novelSS” but otherwise default in running this analysis.

We fit a linear mixed-effects model of intron retention ratios, using gene ID as a random effect, to estimate the variance in intron retention explained by differences among genes and to calculate the intra-gene correlation.

Branchpoint-sequence detection

We extracted introns from annotated CDS features in GTF files and retrieved the corresponding 50-nt upstream branchpoint windows from the reference genome, outputting unique sequences in FASTA format for downstream analyses. These sequences were then analyzed with DREME [38] to identify enriched short motifs. In each case, the most significantly enriched motif is consistent with expected branchpoint sequences in these species.

Methylation and chromatin accessibility data

For M. pusilla, we complimented gene expression estimates and intron retention rates with MNase coverage and CG methylation fractions obtained from [22] under Gene Expression Omnibus accessions “GSM1134620” and “GSM1134621”, respectively.

Random forest modeling

We trained a Random Forest classifier to predict intron retention events using a comprehensive set of genomic and transcriptomic features. Because intron retention ratios are highly skewed and difficult to model accurately using regression, we framed the problem as a binary classification task. The classification target was defined as introns with retention ratios ≥0.025 (misspliced) versus <0.025 (properly spliced). Features included transcription levels (TPM), intron length, GC content, distance from the transcription start site, canonical splice site identity, frame disruption status, and one-hot encoded nucleotide sequences flanking the 5′ and 3′ splice sites (excluding the canonical dinucleotide positions). To prevent data leakage, we performed gene-level splitting, randomly partitioning genes into 80% training and 20% test sets such that all introns from a given gene were assigned to the same partition.

We optimized Random Forest hyperparameters (number of trees, maximum depth, minimum samples per split/leaf, and maximum features per split) using randomized search with 30 iterations and 3-fold grouped cross-validation on the training set, using ROC-AUC as the optimization metric. Class imbalance was addressed using balanced class weights. To identify a minimal feature set without sacrificing predictive performance, we implemented an automated feature selection procedure. We ranked features by their importance scores from the full model, then incrementally added in order of importance while tracking cumulative ROC-AUC on the held-out test set. We determined the optimal feature set using elbow detection via the Kneed algorithm. We produced partial dependence plots to visualize the marginal effect of each selected feature in the resulting minimal models on predicted intron retention probability. We also trained Random Forest models on autosomal introns using the same procedure. The extreme rarity of retention events in autosomal regions prevented effective model training, yielding low F1 scores (S4 Table).

M. pusilla culture

The isolate of M. pusilla used here for RNA and chromatin analysis was obtained from the Roscoff Culture Collection. RCC834 is the same as CCMP1545, initially sequenced by Worden and colleagues [7] and widely used as a reference for this species.

We grew cultures in low salinity (0.8×) artificial seawater with L1 supplements (except that Na₂SiO₃ was not added) [39], made using the L1 kit purchased from the National Center for Marine Algae/Bigelow (cat. No. MKL150L). We maintained small (10 ml) cultures in vented T25 flasks at 18 °C under a 14h light:10 h dark illumination cycle at 220 microEinsteins (μE).

To obtain cell pellets for RNA extraction, we transferred approximately 2 mL of the RCC834 culture into eight separate wells and maintained them in the above conditions for 2 days. On the third day, at 3 hours into the light cycle (1,020 lux, 14 μE), we removed and pelleted each sample by centrifugation for 4 min at 5,000g and 4 °C, immediately flash froze the pellets in a dry ice-ethanol bath, and stored them at −80 °C.

Full-length isoform sequencing with R2C2 and analysis with mandalorion

We extracted RNA from approximately 1E8 cells of M. pusilla strain RCC834 (also known as CCMP1545) using a mixture of Trizol and chloroform, followed by purification with an RNeasy extraction kit. We prepared cDNA libraries through two sequential reactions. First, we performed reverse transcription on 200 ng of total RNA using Smartscribe reverse transcriptase, FS buffer, DTT, indexed oligo-dT primers, and a template switching oligo (TSO), incubating at 72 °C for 3 min and 42 °C for 60 min. Next, we amplified the resulting cDNA by PCR using KAPA HotStart master mix with primers targeting the oligo-dT and TSO sequences.

We circularized the cDNA libraries via Gibson assembly (NEBuilder HiFi) using a short, pre-annealed double-stranded DNA splint overlapping the cDNA ends. To remove uncircularized molecules, we digested the reaction with ExoI, ExoII, and Lambda Exonuclease (all NEB) for 16 hours at 37 °C, then heat-inactivated the enzymes for 20 min at 80 °C. We cleaned the reaction with SPRI beads at a 1:0.85 ratio.

We used the clean, circularized library as a template for rolling circle amplification (RCA) with Phi29 polymerase (NEB) and a random hexamer primer for 18 hours at 30 °C, followed by heat inactivation for 10 min at 65 °C. We then debranched the Phi29 reaction using T7 endonuclease for 2 hours at 37 °C and cleaned and concentrated the product with a NEB Monarch DNA clean and concentrator column. We size-selected the resulting R2C2 library for DNA fragments larger than 3 kb using ProNex beads (Promega) and quantified DNA concentration with a Qubit fluorometer.

S6 Table lists the oligonucleotides used for library construction and circularization. We sequenced the multiplexed R2C2 library on an ONT PromethION R10.4 flow cell and basecalled reads with the Dorado basecaller using model dna_r10.4.1_e8.2_400bps_sup@v5.0.0.

The resulting raw reads were converted into consensus reads and demultiplexed using C3POa (v3.2). Based on these demultiplexed consensus reads, Mandalorion (v4.6) was used to identify and quantify high-confidence isoforms.

We refined the initial Mandalorion isoform calls using “sqanti3_qc” [40,41] using default parameters with the addition of “--isoAnnotLite” to attempt to liftover CDS features from the original gff3 file.

Junction-based alternative splicing classification for R2C2 isoforms

For each M. pusilla gene with at least two R2C2 isoforms, we classified alternative splicing events by comparing junction coordinates across isoforms: alternative 5′ splice sites (A5SS) and 3′ splice sites (A3SS) were identified as junctions sharing one boundary but differing at the other, exon skipping (SE) as junctions in one isoform spanning an exon present in another, and retained introns (RI) as junctions fully contained within an exon of another isoform. This analysis was restricted to the Mandalorion/SQANTI3-filtered isoform set. Per-region counts are provided in S8 Table.

Isoform usage analyses

We quantified isoform usage for each gene using Evenness, defined as the Shannon diversity index of transcript usage normalized by the maximum possible Shannon diversity for the observed number of isoforms. We fit linear models with log-transformed gene expression (log(TPM + 1)) and intron number as covariates, and we included mating-type region status (inside versus outside) as the main predictor. We tested the significance of each predictor, and the coefficient for mating-type region was highly significant, indicating that isoform usage is more even for genes inside the mating-type region after correcting for expression and intron number.

Disruptions to protein function

We quantified possible effects on protein function in two ways. First, we compared the length of CDS annotated for each isoform with the CDS length in the canonical gene annotation. The expectation is that shorter CDS are less likely to be functional. Second, we first used interproscan (version 5.59-91.0) [42,43] with the flag “-goterms” to annotate functional domains and ontologies in all M. pusilla proteins. Then, to assess the effect of isoforms on protein function, we compared the proportion of isoforms that contained X percent of the total functional domains present for each gene across multiple values of X. We performed the same analysis for gene ontology terms. We then compared these results between genes in the MT region and the rest of the genome.

Logistic regression for internal priming

To test whether oligo-dT mispriming on AT-rich sequences could explain elevated intron retention in the mating-type region, we fit a logistic regression of the form P(has_IR) ~ perc_A + gene_id, where perc_A is the percent adenine in the 20 bp genomic window downstream of each isoform’s transcript termination site and gene_id is a categorical fixed effect. The model was restricted to the 164 mating-type genes containing both IR and non-IR isoforms (n = 1,393 isoforms).

S3 Fig. Nucleosome occupancy and methylation patterns in M. pusilla introns and flanking exons.

(Top) MNase-based library preparation coverage within introns. (Middle) Average per CG dinucleotide modification rates in bisulfite sequencing data. (Bottom) Number of CG dinucleotide sites in the 50 bp exons flanking each intron. We used a Mann–Whitney U test to determine p-values. MNase and methylation data from [22]. Underlying data for this figure are available in the GitHub repository at https://github.com/russcd/mating-type-missplicing (file: data/Mpusilla.features.methylation.tsv).

https://doi.org/10.1371/journal.pbio.3003823.s003

(TIF)

S4 Fig. Individual intron retention patterns across the chromosome containing the mating-type region of Mamiellales species.

(Left), per intron retention frequency along autosomal regions (blue) and within the mating-type region (orange). (Right) Retention frequency for each intron plotted along each chromosome with colors as on the left panel. We obtained coordinates for mating type regions from [6]. Tick marks indicate the midpoint of chromosome in each genome assembly. Underlying data for this figure are available in the GitHub repository at https://github.com/russcd/mating-type-missplicing (file: data/{species}.features.branchpoints.tsv).

https://doi.org/10.1371/journal.pbio.3003823.s004

(TIF)

S3 Table. Consensus branchpoint sequences were identified genome-wide using DREME motif discovery.

‘Present’ = exact 6/6 nucleotide match within the 50-nt branchpoint window upstream of the 3′ splice site. Ambiguity codes: V = A/C/G, R = A/G, S = C/G. Branchpoint is underlined.

https://doi.org/10.1371/journal.pbio.3003823.s008

(DOCX)

S5 Table. Protein structure preservation in annotated isoforms.

Columns 2–4: proportion of isoforms retaining ≥X% of GO (Gene Ontology) terms assigned to the canonical protein. Columns 5–7: proportion of isoforms retaining ≥X% of annotated protein domains. These represent complementary measures of functional preservation.

https://doi.org/10.1371/journal.pbio.3003823.s010

(DOCX)