Quantification of Arabidopsis total mRNA levels from public sequencing data (GSE110125; Xiaobo et al., 2019) was performed using vast-tools v.5.1. For each Arabidopsis transcript, this tool provides the corrected-for-mappability RPKMs (cRPKMs), which represents the number of mapped reads per million mapped reads divided by the number of uniquely mappable positions of the transcript (Labbé et al., 2012). To compare the expression of wild-type (WT) seedlings treated or not with norflurazon, vast-tools_compare_expr command was employed with the option -norm, which allows a quantile normalization of cRPKMs between samples. Moreover, genes that had read counts < 50 and were not expressed at cRPKM > 5 across all replicates of at least one of the two samples compared were filtered out. Finally, those genes with a fold change of at least 2 between each of the individual replicates from each sample analyzed were defined as differentially expressed genes (NF-regulated genes; Supplemental Table 1 ). See https://github.com/vastgroup/vast-tools for details.

Vast-tools (Tapial et al., 2017) was employed to quantify alternative splicing from each individual sample of three public RNA-seq experiments: GSE110125 (Xiaobo et al., 2019), GSE130337 (Xu et al., 2020) and GSE164122 (Martín and Duque, 2021), which respectively address responses to norflurazon, lincomycin and light in Arabidopsis thaliana seedlings. This tool maps RNA-seq data to the araTha10 library, based on Ensembl Plants v31 and composed by an extended annotation with all exon-exon and exon-intron junction sequences found in the Arabidopsis thaliana genome using a large compendium of RNA-seq datasets (see Martín et al., 2021 for details). The mapping and experimental details of the samples used in the current study are summarized in Supplemental Table 2 . Then, vast-tools quantifies ES, IR, ALT5 and ALT3 of each sample analyzed, and provides the percent of inclusion (PSI; Percent of Spliced In) of each putative alternative sequence using only exon-exon (or exon-intron for IR) junction reads (Braunschweig et al., 2014; Tapial et al., 2017). In addition, vast-tools associates a quality score to each AS event based on the read coverage that sustains the PSI quantification. To define differentially spliced events in response to norflurazon, lincomycin or light, I used the command vast-tools compare adding specific filters. First, VLOW events were discarded, which represent the lowest range of read coverage (for further details see https://github.com/vastgroup/vast-tools). To improve the read coverage that supports each splicing junction and thus decrease the number of AS events discarded by this filter, replicates of each sample were pulled together (using the vast-tools merge function). The consistency of these merged samples was confirmed ( Supplemental Figure 1 ). Moreover, to evaluate up and down sequence inclusion of all ALT events, not only of the most external splice sites, the analysis was conducted with the –legacy_ALT vast-tools compare option. Additionally, the –p_IR and –min_ALT_use 25 options were also applied. The first eliminates those IR events with a significant imbalance between the two exon-intron junctions (P < 0.05; binomial test; see Braunschweig et al., 2014 for details). –min_ALT_use 25 ensures that ALT3 and ALT5 events are located in exons with a sufficient inclusion level, in particular, implies that the host exon has a minimum PSI of 25 in each compared sample. Then, those splicing events with a |ΔPSI| > 15 between the pair of samples being compared were selected as differentially spliced events (NF-, lincomycin- or light-regulated AS events; respectively found in Supplemental Tables 3 - 5 ). Furthermore, genes exhibiting differential AS were considered differentially spliced genes (AS genes). Then, the pertinent control sets were generated to conduct molecular and functional analyses of the AS events and AS genes. First, non-regulated AS events were those that passed the coverage criteria and were not differentially spliced (|ΔPSI| < 15). Genes for which all their AS events belong to the group of non-regulated events constituted the group of non-regulated AS genes; by definition, these genes are all multiexonic.

Three different publicly available RNA-seq experiments were analyzed with vast-tools to quantify the PSI values of NF-regulated AS events in three different conditions. First, GSE130337 (Xu et al., 2020), to address the molecular response to lincomycin in Arabidopsis seedlings ( Figure 1C ). GSE164122 (Martín and Duque, 2021) was used to study light responses of dark-grown seedlings ( Figure 6C ), and finally, GSE41432 (Drechsel et al., 2013), to show the PSI values in WT, upf1, upf3 and upf1upf3 seedlings ( Figure 4C ). The mapping and experimental details of the samples used in the current study are summarized in Supplemental Table 2 .

For each AS event and gene, the ΔPSI and the logarithmic to base 2-fold change were calculated comparing either WT or gun1 norflurazon-treated seedlings with untreated WT seedlings. This value represents the magnitude of change in expression or splicing in response to retrograde signals. To compare WT and gun1 values, the WT value was set to 100, to then calculate the percentage of response in the gun1 based on this equivalence. This is, gun1 value multiplied by 100 and then, divided by the WT value. The resulting percentage indicates how much of the molecular response to norflurazon occurring in the WT is kept in the gun1 mutant.

Alternatively spliced sequences were classified based on their localization within the gene structure: 5´-UTR, CDS and 3´-UTR. In addition, the predicted impact on the canonical open reading frame was mainly categorized in function of whether or not they preserve the reading frame (i.e., they have length multiple of three nucleotides) and contain in-frame stop codons predicted to trigger NMD (stop codons considered premature because they are located at least 50 nucleotides upstream of an exon-exon junction; see Martín et al., 2021 for details). These alternative possibilities generate three different scenarios, AS events that produce alternative protein isoforms, or unproductive transcripts when the alternative sequence is included or excluded. This information was retrieved from the download section of PastDB (Legacy version; http://pastdb.crg.eu/wiki/Downloads; Plant alternative splicing and transcription Data Base). Moreover, comparison of exon and intron features associated with splicing regulation was performed with Matt v1.3.0 (Gohr and Irimia, 2019). For each compared group of exons and introns events: upregulated, downregulated and non-regulated, this tool obtained and compared the exon and intron length, the GC content and the splice site strength (see Martín et al., 2021 for details regarding calculations of splicing site strength).

NF-regulated AS genes were functionally categorized according to their subcellular localization and gene ontology (GO) annotation. A single subcellular localization was assigned to each locus based on the annotation of the Arabidopsis genome (Araport 11) and percentages of each type were plotted discarding those genes with unknown subcellular localization. To identify enriched GO biological processes, molecular functions and cellular components, analyses were performed using the functional annotation classification system DAVID (Huang et al., 2007). The same strategy was used to define GO category enrichment of lincomycin-regulated AS genes.

To define how retrograde signals control splicing, the AS profiles of WT seedlings grown in the presence or absence of norflurazon were compared (GSE110125; Xiaobo et al., 2019). This herbicide inhibits carotenoid biosynthesis and causes accumulation of Mg-protoporphyrin IX, which acts as a retrograde signaling (Oelmüller et al., 1986; Strand et al., 2003). First, vast-tools (Tapial et al., 2017) was used to quantify steady-state mRNA levels in each sample (hereafter gene expression; GE), and to determine changes in expression levels of representative PhANGs (Allen et al., 2003). As expected, all PhANGs were consistently downregulated in the WT in response to norflurazon ( Supplemental Figure 2 ). Moreover, in line with the key role of GUN1 mediating the transcriptional changes induced by RS (Koussevitzky et al., 2007), this downregulation was weaker in the gun1 mutant ( Supplemental Figure 2 ). Additionally, 94.5% of the genes described as NF-regulated in this analysis ( Supplemental Table 1 ) were also identified as NF-regulated in Xiaobo et al., 2019 ( Supplemental Figure 3 ). These results thus confirm the quality of this GE quantification.

Then, vast-tools was used to define alternatively spliced genes. This tool quantifies the sequence inclusion levels for all four major AS types (IR, ES, ALT5, ALT3; Figure 1A and Material and Methods for details). For each sample and splicing event in the genome, vast-tools provides the percent of alternative sequence inclusion (using the PSI metric), which corresponds to the percentage of expressed transcripts that include the alternative sequence. The comparison between WT seedlings treated or not with norflurazon retrieved 252 differentially spliced AS events (|ΔPSI| > 15; Figure 1A and Supplemental Table 3 ) in 205 genes. As commonly observed in plants, the majority of these events were IR and ALT3 events (Martín et al., 2021). For 165/252 AS events, the alternative sequence was more included in the mRNAs in the presence of norflurazon (NF-upregulated), while for the remaining 87 AS events it was less included (NF-downregulated) ( Figure 1B ). Interestingly, AS changes were also milder in the gun1 mutant ( Figure 1B and Supplemental Figure 4 ). In fact, 66% of the NF-regulated AS events in WT plants were not differentially spliced when comparing untreated WT seedlings with NF-treated gun1 mutants (|ΔPSI| < 15), indicating the implication of this protein in triggering the downstream AS regulation. As observed in Supplemental Figure 5 , the molecular response to norflurazon in the gun1 mutant was around half of the response observed in the WT, both for GE and AS. This finding indicates that similarly to the transcriptional control that RS exerts in the nucleus, GUN1 is also key for the nuclear control of AS. Finally, to further confirm that retrograde signals impact splicing of these AS events, their PSI was compared between samples treated or not with lincomycin, an herbicide also known to activate retrograde signals by inhibiting plastid translation (Oelmüller et al., 1986). This analysis showed that inclusion of AS events upregulated or downregulated by the effect of norflurazon was similarly regulated by lincomycin ( Figure 1C ). Therefore, these data demonstrate that retrograde signals control splicing of a subset of transcripts in a GUN1-dependent manner.

Next, the proportion of AS events whose inclusion is either up- or downregulated by retrograde signals among each type of AS event was calculated. This uncovered that most of the IR events are upregulated in response to norflurazon, while for the other types of AS events this proportion is similar ( Figure 2 ), meaning that RS has a major role in enhancing the inclusion of intronic sequences. Accordingly, this pattern was also characteristic of lincomycin-treated seedlings ( Supplemental Figure 6 and Supplemental Table 4 ). Because genomic features such as the intron and exon length, their GC content or the strength of their splice sites modulate their splicing, characterization and comparison of these regulatory features between NF-regulated introns and exons with those non-regulated was conducted using the Matt software. This study revealed that NF-upregulated introns, those whose splicing fails in the presence of norflurazon (retained), are longer and have weaker 5’ splice sites ( Figure 3A ; Supplemental Figures 7, 8 ), characteristics known to difficult their spliceosomal recognition (de Conti et al., 2013). On the other hand, exons skipped in response to norflurazon (NF-downregulated) tend to be shorter and are surrounded by longer introns ( Figure 3B ; Supplemental Figures 7, 8 ), a pattern also known to hamper their recognition (Hollander et al., 2016). Also, according to their alternative nature, both subsets of NF-regulated exons have weaker splice sites ( Figure 3B ; Supplemental Figures 7 , 8 ). Moreover, the GC content of differentially spliced introns was higher than that of non-regulated introns ( Figure 3A ; Supplemental Figures 7 , 8 ), which implies a reduction of the existing difference between the GC content of introns (lower) and exons (higher), and therefore an additional complication for their recognition (Amit et al., 2012). This reduction is also characteristic of differentially spliced exons, especially of the NF-upregulated subgroup ( Figure 3B ; Supplemental Figures 7 , 8 ). Thus, this analysis detected the genomic particularities shared by the RS-regulated AS events, and revealed that retrograde signals impact the spliceosome’s capacity to properly splice introns and exons imbibed in complex genomic contexts.

Then, the molecular functions of the AS events upregulated or downregulated in response to norflurazon were addressed. As expected given the dominant proportion of the CDS region within the gene structure, the majority of NF-regulated AS events were located in this region ( Figure 4A ). However, despite this predominant presence in the CDS regions, a significant increase in the proportion of AS events located in the UTRs, both in the subsets of NF-upregulated and NF-downregulated AS events, could also be detected ( Figure 4A ; two-sided Fisher’s test; P < 0.0012). Focusing on the CDS-located AS events, their impact on the canonical open reading frames (ORFs) was next assessed. In comparison to non-regulated AS events, those events whose alternative sequence is more included in response to norflurazon (NF-upregulated) were enriched for AS events predicted to generate unproductive mRNAs when the alternative sequence is included ( Figure 4B ; i.e., events that introduce a PTC; see Material and Methods for details). Conversely, most of the NF-downregulated AS events generate unproductive mRNAs when the alternative sequence is excluded ( Figure 4B ). Therefore, 67% of the CDS-located NF-regulated AS events have the ultimate effect of enhancing accumulation of unproductive mRNAs when retrograde signals take place (repressive AS events; Figure 4B ), compared to only 25% predicted to generate alternative protein isoforms, and to a minimal group of AS events (8%) that generate unproductive mRNA isoforms in the absence of norflurazon ( Figure 4B ). Accordingly, among the three categories, repressive AS events also represent the majority of the CDS-located lincomycin-regulated AS events (53%). Because repressive AS events are mainly characterized by carrying NMD-eliciting features, their inclusion levels (PSI) in the Arabidopsis upf mutants were quantified (GSE41432; Drechsel et al., 2013). UPF1 and UPF3 genes encode for the two key NMD factors (Hori and Watanabe, 2005; Arciga-Reyes et al., 2006). This analysis globally revealed higher inclusion levels of repressive AS events in these mutants, particularly in the double upf1upf3 mutant ( Figure 4C ), meaning that their transcript isoforms, accumulated in response to RS, are degraded through the NMD surveillance pathway. In contrast, as expected for AS events predicted to generate alternative proteins, these isoforms lacked NMD regulation ( Supplemental Figure 9 ). Therefore, this result implies that RS, through the action of repressive AS events, post-transcriptionally downregulates the mRNA levels of specific transcripts. In agreement, particularly genes harboring repressive AS events are significantly enriched for genes whose mRNA levels are differentially downregulated by norflurazon ( Figure 4D ; Supplemental Figures 10 , 11 ; two-sided Fisher’s test; P < 2.5e-5; see Material and Methods for details on the GE analysis).

To gain insight into the biological functions of the 205 genes whose splicing is regulated by norflurazon, a gene ontology (GO) analysis was first conducted. This revealed significant enrichment for chloroplast- and splicing-related GO terms as well as for genes whose expression is light-regulated ( Figure 5A ). A great example of splicing-related genes are those encoding the spliceosome proteins U2AF65A and U2AF65B (Chen and Moore, 2015), and the splicing regulators RS2Z32 (Barta et al., 2010) and GRP8 (Streitner et al., 2012). Additionally, genes from different chloroplast-related categories were found. For example, PDM2, GC1 and PPL1, all of them with key roles in chloroplast development and functioning (Maple et al., 2004; Du et al., 2017; Che et al., 2020); KAC1 and KAC2, involved in the movement of chloroplast through actin filaments (Suetsugu et al., 2010); and also, TIC100, a component of the TIC chloroplast transport complex recently implicated in RS (Kikuchi et al., 2013; Loudya et al., 2022). Supporting the impact that RS exert on the splicing pattern of this type of genes, these GO terms were also enriched among the subset of genes differentially spliced in response to lincomycin ( Supplemental Figure 12 ). Next, classification of NF-regulated AS genes based on the subcellular locations of their encoded proteins revealed overrepresentation of chloroplast-located proteins ( Supplemental Figure 13 ; two-sided Fisher’s test; P < 0.00074). Strikingly, this enrichment was specific of genes harboring repressive AS events ( Figure 5B ), which represent the majority of RS-regulated AS events ( Figure 4B ). Therefore, these results indicate that RS post-transcriptionally controls the mRNA levels of chloroplast proteins.

Retrograde and light signals converge to antagonistically regulate transcription of a subset of common genes (Ruckle et al., 2012; Martín et al., 2016; Xu et al., 2016). Given this preceding information, the molecular link between retrograde and light signaling at the RNA level was investigated. First, by following the same criteria as those established for identifying RS-regulated AS events (see Material and Methods), the splicing profiles of 3-day-old WT seedlings grown either in continuous darkness or white light were compared (GSE164122; Martín and Duque, 2021). This comparison resulted in the definition of 309 light-regulated AS events ( Supplemental Table 5 ). Importantly, the GO analysis of genes harboring these AS events also revealed enrichment for chloroplast proteins ( Figure 6A ). Moreover, the study of the predicted impact on the ORFs of the 246 light-regulated AS events located in the CDS unveiled that light enhances inclusion of alternative sequences generating unproductive mRNAs when the sequence is excluded, while it downregulates inclusion of unproductive sequences when it is included ( Figure 6B ). Hence, a minority of AS events produce unproductive mRNAs in the light (repressive AS events; Figure 6B and Hartmann et al., 2016), while most do so in the presence of norflurazon (repressive AS events; Figure 4B ). This result indicates that RS- and light-mediated AS exert opposite molecular functions, respectively enhancing and diminishing the accumulation of unproductive transcripts. Interestingly, some of the NF-regulated AS genes mentioned, such as U2AF65B, GRP8 and PPL1 have also been characterized as being differentially spliced between dark and light conditions (Hartmann et al., 2016). To assess whether light globally influences NF-regulated AS events, their inclusion levels were quantified in 3-day-old WT and pifq seedlings grown in continuous darkness or light (GSE164122; Martín and Duque, 2021). NF-downregulated and NF-upregulated AS events were respectively upregulated or downregulated by light in WT seedlings ( Figure 6C ). Further validating the influence of light in the AS of these events, the inclusion levels of these AS events had the same trend in dark-grown pifq seedlings, known to mimic light responses both at the molecular and phenotypic level (Leivar et al., 2008; Leivar et al., 2009). This result firmly revealed that RS-regulated AS events are antagonistically regulated by light, which is in line with the opposite molecular functions of RS- and light-regulated AS events ( Figures 4B , 6B ), and also, with the common enrichment of chloroplast genes among the sets of RS- and light-regulated AS genes ( Figures 5A , 6A and Supplemental Figure 12 ).