A Genome-Wide Alternative Splicing Landscape Speciﬁcally Associated with Durable Rice Blast Resistance

: The rice blast, caused by the hemibiotrophic plant pathogen Magnaporthe oryzae , is a devastating disease that threatens rice crop production worldwide. The molecular interactions that underlie the rice- M. oryzae interaction have received much attention. However, genome-wide research focusing on alternative splicing (AS) has not been well-studied in rice— M. oryzae interactions. AS in plants leads to diverse proteomes without an expansion in gene numbers to regulate cellular processes during abiotic or biotic stress. The Pi21 gene negatively regulates rice resistance to M. oryzae infection, and thus the Pi21 -RNAi silenced transgenic line (#241) exhibits partial but durable resistance. We compared the AS landscape in #241 and “Nipponbare” (Nip) during interacting with M. oryzae Guy11, and the alternative 3 (cid:48) splice-site (A3SS) is the most common AS type. GO enrichment analysis of #241-speciﬁc differentially alternatively spliced genes (DASGs) revealed that WRKY transcription factors (TFs), bHLH TFs, F-box protein with leucine rich repeats, AAA-type ATPase, and protein kinase were enriched in the GO terms “response to jasmonate acid (JA)” and “ethylene (ET)” at 24 h post-inoculation (hpi). At 48 hpi, one #241-speciﬁc DASG, ubiquitin gene ( Os08g0295000 ), was predicted to be involved in endoplasmic reticulum (ER) stress. In silico analysis combined with PCR ampliﬁcation conﬁrmed that multiple isoforms are produced by Os08g0295000 and a skipped exon (SE) event results in isoform switching during interaction with M. oryzae . Protein–protein interaction (PPI) network analysis predicted that Os08g0295000 -encoding proteins may interact with SNARE protein Q9LGF8 (Uniprot ID) to cooperatively regulate rice’s response to M. oryzae . This study uncovered the AS proﬁle of rice in response to M. oryzae , which will help to explore the linkage between AS and durable rice resistance.


Introduction
Alternative splicing (AS) describes the mechanism by which multiple transcript isoforms are produced from single precursor mRNAs (pre-mRNAs) [1], and results mainly from the different combinations of exons for inclusion into the mRNAs. AS increases the proteome diversity without massive gene gain events [2]. The translation of alternatively spliced transcripts determines which domains are present in the corresponding proteins, and this can result in a diverse higher structure of proteins with different functions [3]. The spliceosome is a ribonucleoprotein complex consisting of different RNP subunits that can accurately recognize exons and introns to carry out the splicing reaction, resulting in AS [4]. Genome-wide studies have shown that AS occurs widely in eukaryotes. More than 95% of human multi-exon genes undergo AS [5]. In Drosophila melanogaster, alternative splicing results in tissue-or growth stage-specific protein isoforms with different functions in development [6]. Based on high-throughput sequencing, 61% of intron-containing Arabidopsis genes undergo AS events under multiple conditions [7], such as development and abiotic

Identification of Rice Differentially Alternatively Spliced Genes (DASGs) during Interaction with Magnaporthe oryzae
To investigate transcriptional changes and the AS scheme of defense-associated genes, we first collected the publicly-available transcriptomes of the Pi21-RNAi transgenic rice line (#241) and the susceptible cultivar "Nipponbare" that had been infected by M. oryzae Guy11 (SRA492222) [19]. It has been suggested that Pi21 negatively regulates rice blast Agronomy 2022, 12, 2414 3 of 15 resistance, and plants carrying a loss-of-function mutation of Pi21 express durable and non-race-specific rice blast resistance [18]. In this study, the transcriptome data used for AS detection included the following samples: (i) #241 infected by M. oryzae Guy11 at 24 and 48 hpi (sample names: Guy11_241_24 h and Guy11_241_48 h); (ii) Nip infected by the M. oryzae strains Guy11 at 24 and 48 hpi (sample names: Guy11-Nip-24 h, Guy11-Nip-48 h); (iii) the un-infected samples of #241 and Nip were defined as the controls (sample names: Guy11-241-0 h and Guy11-Nip-0 h). Due to the pathogen perception by rice plants that occurs within 48 h [21], selection of the 24 and 48 hpi timepoints will help us to uncover links between AS and activation/inhibition of rice innate immunity.
We designed a transcriptome assembly pipeline to detect AS events (Figure 1a). Of the trimmed clean reads, 91.58% to 93.93% mapped to the rice reference genome with HiSAT2 (step II) ( Table 1). Using the mapped reads, transcripts were assembled with StringTie. A total of 68,103 transcripts (across 39,231 genes) were yielded and 59,026 of them (across 33,902 genes) had coding potential, which were retained for the next step (step III, IV). Using gffcompare, retained assembled transcripts were compared with annotated transcripts, which produced comparison class codes for distinguishing each other following the method proposed by [22]. In total, we identified 39,224 annotated transcripts (from 32,562 reference genes), 18,650 novel isoforms (from 9929 reference genes), and 1054 novel isoforms (from 910 novel genes) (Figure 1b). Based on the expression level of the obtained transcripts/isoforms and the corresponding genes, we identified 7802 differentially expressed genes (DEGs) across all samples by DEGseq ( Figure 1a) (step V), and 2504 of them were also predicted as alternatively spliced genes (ASGs) by SUPPA2 [23], which were defined as differentially alternatively spliced genes (DASGs).  (7) alternative last exon (AL). As shown in Figure 1d, alternative donor (A3SS) was the dominant AS type in all samples, followed by retained intron (IR), alternative acceptor (A5SS), and skipped exon (SE).

Comparative Analysis Uncovered
To identify #241-specific DASGs, we compared DASGs identified in #241 and Nip. As Figure 2 showed, 226 DASGs were specific to #241 at 24 hpi (Table S2) and 215 at 48 hpi (Table S3). The AS events that occurred in the #241-specific DASGs may involve general response mechanisms that are specific to resistible cultivars. Thus, these #241-specific DASGs were absent in Nip and were retained for further analysis to explore the linkage between AS events and durable rice blast resistance.    To identify #241-specific DASGs, we compared DASGs identified in #241 and Nip. As Figure 2 showed, 226 DASGs were specific to #241 at 24 hpi (Tables S2) and 215 at 48 hpi (Tables S3). The AS events that occurred in the #241-specific DASGs may involve general response mechanisms that are specific to resistible cultivars. Thus, these #241-specific DASGs were absent in Nip and were retained for further analysis to explore the linkage between AS events and durable rice blast resistance.

Investigation of Infection-Specific AS Isoforms Generated by #241-Specific DASGs
During the analysis of the expression level of #241-specific DASGs-generated transcripts/isoforms at each stage, we introduced the concept of the transcript usage index (TUI), which indicates the expression level of individual transcripts compared to the

Predicted the Protein Domains Alterations Caused by AS Occurred in #241-Specific DASGs
Considering that AS may change the protein domain [15], we investigate the domain architecture of transcripts/isoforms produced by #241-specific DASGs. There were 68 #241-specific DASGs that underwent domain alterations due to AS, which corresponds to 312 transcripts/isoforms (Table S6). Among them, 56 #241-specific DASGs exhibited domain loss or gain. For example, three annotated transcripts (Os06t0227200-01, Os06t0227200-02, and Os06t0227200-03) of Os06g0227200, encode proteins with Dev_Cell_Death (PF10539) domains, which are absent in the novel isoform MSTRG.20959.2. A similar phenomenon was also detected in the bHLH TF gene (Os01g0952800) and the Glutathione S-transferase gene (Os12g0210300). Twelve #241-specific DASGs underwent domain changes, such as a gene encoding PPR repeat-containing protein (Os02g0829850). In addition, one gene encoding an EF-hand motif-containing protein, Os02g0802400, was found to have undergone domain gain.

Expression Analysis of Rice Spliceosome Involving in Regulation of AS Events in #241 during Interaction with M. oryzae
In Eukaryotes, mature mRNAs are generated from precursor mRNAs (pre-mRNAs) through RNA splicing. This process is performed by the spliceosome, a highly dynamic ribonucleoprotein complex that contains small nuclear ribonucleo-proteins (snRNPs) and non-snRNP proteins [28]. To study spliceosomes involved in the regulation of AS, we retrieved annotated rice spliceosome protein sequences from the PlantGDB database (http://www.plantgdb.org/, accessed on 11 April 2022). A BLASTP search was then used to check for putative spliceosomes missing from the database. In total, 281 rice spliceosome coding genes were obtained (Table S7). Among them, 60, 117, 42, and 62 genes were assigned into subfamilies of the small nuclear ribonucleoprotein (snRNP) subfamily, the splicing factor subfamily, the splicing regulation subfamily, and the novel spliceosome proteins subfamily, respectively (Table S7). Thirteen and four spliceosome genes were significantly differentially expressed at 24 and 48 hpi ( Figure S3a,b). Based on their expression analysis, eight spliceosome genes were significantly up-regulated at 24 hpi ( Figure S3c). Os02g0150100, a gene that encodes one U4/U6, was up-regulated by 4.66-fold, followed by Os03g0569900 and Os02g0767100, which encode one HSP73/HSPA8 and one CBP20, respectively. At 48 hpi, four up-regulated spliceosome genes were identified, which included Os03g0110400 (encoding SF2/ASF protein) and Os04g0504800 (encoding NPW38 protein) ( Figure S3b).

Discussion
Alternative splicing (AS) occurs widely in eukaryotes and shapes the diverse proteome without gene expansion. At present, AS events are reported to be involved in multiple biological processes in rice, such as development and the responses to abiotic or biotic stresses [29][30][31]. However, it is still unclear about the genome-wide AS profile of rice in response to M. oryzae. In this study, we designed a transcriptome assembly pipeline to combine public RNA-seq data and explored the relationship between AS and durable rice resistance through the investigation and comparison of differentially alternatively spliced genes (DASGs) in Pi21-RNAi transgenic line (#241) and Nip.
Alternative donor (A3SS) splicing is the most common AS type in both the susceptible and resistible rice cultivars during interaction with M. oryzae Guy11, which is consistent with the results of previous research [32]. However, [26] found that intron retention (RI) was the most common AS type in M. oryzae during its infection procession. Thus, the different AS preferences of rice and M. oryzae were discovered in this study. Totals of 187 and 191 #241-specific DASGs were found at 24 and 48 hpi, which is more than were found in Nip. This suggests that AS may play an important role in durable rice resistance at an early interaction stage.
In detail, #241-specific DASGs at 24 hpi encoding protein kinase, leucine-rich repeat containing proteins, bHLH TFs and WRKY TFs, and ATPases were found involved in responses to plant hormone stimulation, such as jasmonic acid (JA) and ethylene (ET), which play an essential role in transducing the activation of plant defense systems [33]. We thus inferred that these DASGs might be related to the regulation of response to plant hormones.
According to [30], TF genes of WRKY and bHLH were also found to undergo frequent AS events in response to salt stress and pathogen infection. For example, [15] proposed that two WRKY TFs of the japonica rice cultivar "Zhonghua 17" (ZH17), OsWRKY62 (Os09g0417800) and OsWRKY76 (Os02g0181300), generate short alternative variants to enhance rice resistance against M. oryzae. However, we could not find OsWRKY62 and OsWRKY76 in the list of #241-specific DASGs, which may originate from different genetic backgrounds of the cultivars used. Similar regulation was also reported in the LAMMER kinase gene OsDR11 (Os12g0460800), which produced two alternative variants, with contrasting functions, in response to Xoo infection [34]. We also found that this gene underwent an A3SS AS event and a retained intron in response to M. oryzae Guy11 infection (Table S2). Alternative splicing of OsDR11 (Os12g0460800) may play an important role in rice defense against pathogens. Furthermore, Os01g0952800, encoding bHLH TF, generates three alternative variants with down-regulated expression levels, which may be involved in ethylene-activated signaling pathways. Notably, the wheat bHLH TF gene TabHLH060 was also inhibited by ethylene and enhanced the susceptibility of transgenic Arabidopsis thaliana [35]. Thus, we suggested that Os01g0952800 was inhibited by ethylene, which might be related to its alternative splicing events. Detailed annotation of Os05g0425700 showed that a signature F-box motif at the N-terminus and leucine rich repeats (LRR_6) at the C-terminus, which represents that Os05g0425700 encodes an F-box protein. In Arabidopsis, F-box protein Coronatine insensitive 1 (COI1) is a pivotal factor in the JA signal response [36]. However, how AS affects the F-box protein-regulating JA signal response is still unknown and this study provides insight into this aspect.
Among #241-specific DASGs at 48 hpi, Os08g0295100, a ubiquitin-like protein encoding gene, underwent an SE AS event in response to M. oryza infection and may be involved in response to endoplasmic reticulum (ER) stress. Three classes of ubiquitination pathway-associated enzymes were summarized by [37]. A BLAST search against the Uniprot database revealed that Os08g0295100 encodes an E3 ubiquitin ligase. Under diverse stress, the demand for plant protein exceeds the plant system capacity, which results in the accumulation of misfolded or unfolded proteins and sets off ER stress [38]. In response to ER stress, unfolded protein response (UPR) was activated to lighten the load of unfolded proteins in ER, during which E3 ubiquitin ligase is required for marking unfolded proteins with adapters so that the ER-associated degradation (ERAD) system removes these proteins. In this study, we confirmed that multiple isoforms were generated by Os08g0295100, and primary isoforms switched from long-spliced variants (MSTRG.25458.2: 1020 bp) to short-spliced variants (MSTRG.25458.1: 630 bp), which may affect ubiquitination of unfolded target proteins and durable rice resistance. Furthermore, the PPI network of Os08g0295100 predicted that Os08g0295100 might interact with a SNARE protein (Uniprot ID: Q9LGF8). In Arabidopsis thaliana, SNARE protein SYP61 and ubiquitin ligase ATL31 cooperatively regulated the response to carbon/nitrogen conditions, during which SYP61 interacted with ATL31 and was ubiquitinated [39]. Thus, we suggested a similar cooperation with Os08g0295100 and Q9LGF8 may also play a role in response to rice blasts. However, the relationship between AS events of Os08g0295100 and its function still needs further research.
Through investigation of infection-specific AS isoforms generated by #241-specific DASGs, the expression levels of 19 and 26 transcripts/isoforms were specifically dominant (TUI value > 0.5) at 24 and 48 hpi in the #241-M. oryzae interaction. Os06t0227200-03, the infection-specific primary transcript generated by Os06g0227200 at 24 hpi, encodes protein with the Dev_Cell_Death (DCD) domain, which was originally detected in the soybean NRPgene sequence and mediates signaling in plant development and programmed cell death (PCD) [40]. Interestingly, an A3SS AS event results in a primary transcript/isoform switch from MSTRG.20959.2 to Os06t0227200-03 during interactions with rice blasts (Figure S2a,c), which results in DCD domain gain (Table S6). Hence, it is reasonable to infer that M. oryzae-induced AS events enhance the Os06t0227200-03 production, which activates PCD to limit the spread of M. oryzae invasive hyphae. Moreover, Os07g0495900 encodes AAA-type ATPase and underwent an A5SS AS event at 24 hpi, which resulted in primary transcript/isoform switching from the novel isoform MSTRG.23603.1 to the annotated transcript Os07t0495900-02. Os07t0495900-02 was also associated with a response to JA (Figure 3a). The AAA-type ATPase from tobacco (Nicotiana tabacum), NtAAA1, was induced by JA and ethylene during the hypersensitive response (HR) resulting from infection by tobacco mosaic virus (TMV) and Pseudomonas syringae [41]. Thus, we suspect that JA may stimulate AS events, as shown by the induction of Os07t0495900-02.
Taken together, TFs of WRKY and bHLH, F-box protein with leucine rich repeats, and AAA-type ATPase underwent AS events and may be involved in the response to the stimulation of plant hormones, especially JA and ET. More importantly, we confirmed that skipped exon AS events in the E3 ubiquitin ligase coding gene Os08g0295100 enhanced the generation of its short isoform MSTRG.25458.1 during interaction with M. oryzae. In silico analysis revealed that Os08g0295100 may interact with the SNARE protein Q9LGF8, which might cooperatively regulate rice defense. This study will contribute to a greater understanding of the AS landscape in rice blast-susceptible and resistible Oryza sativa, which helps to explore AS events associated with durable rice resistance.

Data Collection
Transcriptome data of Pi21-RNAi transgenic rice line (#241) and "Nipponbare" at 24 and 48 hpi of interacting with M. oryzae Guy11 were retrieved from the Ensembl website with accession of SRA492222, which includes millions of 50-bp single-ended RNA-seq reads. Detailed information of the collected data is provided in Table S1. Oryza sativa reference genome IRGSP-1.0 was downloaded from the Ensembl website (http://plants. ensembl.org/Oryza_sativa/Info/Index, accessed on 1 April 2022). FastQC v0.11.8 and Trimmomatic v.38 were used to assess read quality and to remove poor-quality reads or reads that consisted of adapters only.

Parse of Gffcompare Class Code to Distinguish Annotated Genes and Novel Genes/Isoforms
Clean DNA sequencing reads were mapped to the annotated rice genome using HiSAT2 with the parameters: no-mixed and no-discordant. The aligned reads for each sample were used to assemble the transcripts using StringTie v 2.1.1. The assemblies produced by StringTie were merged with the reference annotation file into one GTF file using the merge command in StringTie. The Transdecoder v5.5.0 was used to check the coding potential of the assembled transcripts among the StringTie-merged GTF files, and the assembled transcripts without protein-coding potential were filtered out. The remaining assembled transcripts were compared with the annotated transcripts for the classification of novel genes or novel transcripts using gffcompare v0.11.5, which was based on the gffcompare class code parse method proposed by [22].

Quantification of Transcript Abundances and Identification of Differentially Expressed Genes (DEGs)
The determination of the relative abundance of the annotated and assembled transcripts at different host-pathogen interaction stages was performed using StringTie, which outputs FPKM (fragments per kilobase of transcript per million mapped reads) values for each transcript. The expression profiles of each transcript were then converted into gene abundance using TBtools [42], with StringTie-merged GTF files as input. The differential gene expression analysis was performed using DEGseq. Genes with a combination of an adjusted p-value < 0.05 and |fold change| ≥ 1 were considered to be differentially expressed.

Identification of Genes That Undergo AS Events
Among the AS detection tools available, SUPPA2 can give an effective and accurate prediction of AS events, especially at low sequencing depth and with short reads. In this study, we used SUPPA v2.3 to generate the AS profiles of genes in the StringTie-merged GTF file, including annotated genes and novel genes based on the transcript assembly results. SUPPA first generates annotated files for seven types of AS events (A5SS, A3SS, SE, RI, MXE, AFE, and ALE). SUPPA2 then quantifies the percentage spliced index (PSI; ψ) for each sample, which indicates the inclusion of sequences into transcripts based on the normalized transcript abundance values (FPKM) and thus represents the AS event inclusion level. Transcripts with FPKM values < 1.0 were eliminated. The difference of splicing change (∆PSI; dpsi value) from each AS event across multiple interaction stages and their significance were calculated, which provided high-confidence AS events between the two conditions.

Gene Ontology Enrichment Analysis
The total proteins, encoded by the annotated and StringTie-merged transcripts, were uploaded on the eggNOG site (http://eggnog-mapper.embl.de/, accessed on 5 April 2022) and Gene Ontology (GO) annotation was performed. GO enrichment analyses were performed using TBtools. The GO terms with FDR < 0.05 and corresponding GO levels > 4 were retained and defined as enriched GO terms.

Visualization of the AS Profiles of the DASGs and Protein-Protein Interaction (PPI) Network
The gene models and corresponding Sashimi plots were visualized via the Integrative Genomics Viewer (IGV) (Robinson et al., 2017). Heatmaps were generated with the heatmap.2 function in R. The interacting proteins encoded by DASGs that were unique to the Pi21-RNAi transgenic rice line were identified by the STRING database (https://cn.string-db.org/, accessed on 9 April 2022) using the default parameters.

PCR Amplification
Conidia of Magnaporthe oryzae strain Guy11 were used for leaf inoculation of Tetep plants, another resistible indica rice cultivar [43]. Conidial suspensions were adjusted to 5 × 105 spores/mL, and 5 mls of the suspensions were sprayed onto the leaves of young plants. Uninoculated leaves of "Tetep" plants were used as the control samples, and leaves of "Tetep" plants at 24 and 48 h post-inoculatio (hpi) were used as the treatment samples. PCR amplification followed the method of [44]. Total RNA was isolated from rice leaves using the Qiagen RNAeasy Mini Kit (Qiagen Inc., Valencia, CA, USA) according to the manufacturer's instructions. The purified RNA was converted to cDNA by a Takara PrimeScriptTM RT reagent kit. PCR amplification primers were designed via primer version 4.0 (https://bioinfo.ut.ee/primer3-0.4.0/, accessed on 7 April 2022) and are also provided in Figure S1. PCR amplifications were performed on an Eppendorf G-storm GS2 Mastercycler. The amplification conditions were a 4 min denaturation step at 94 • C, followed by 32 cycles of 30 s at 94 • C, 30 s at 55 • C, and 1 min at 72 • C, with a final extension step of 10 min at 72 • C.
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/agronomy12102414/s1, Table S1. Data accession and mapping rate of data used in this study; Table S2. #241-specific in Guy11_241_24h vs. 0h comparison; Table S3. #241-specific DASGs in Guy11_241_48h vs. 0 h comparison; Table S4. GO enrichment analysis for #241-specific DASGs might be associated with rice defense; Table S5. Annotation of primary transcripts/isoforms generated by #241-specific DASGs; Table S6. #241-specific DASGs underwent domain alterations; Table S7. Overview of rice splicesome used in this study; Figure S1. Multiple sequence alignment of transcript/isoform generated by Os08g0295100. P1 and P2 represent forward and reverse primers used for PCR amplification. L1 and L2 fragments represent exons absent in MSTRG.25458.1; Figure S2. Sashimi plot and protein-protein network analysis for Os06g0227200. (a) Sashimi diagram showing the exon-intron structure of representative annotated transcripts and novel isoforms. The red dashed box encloses the location of AS events. (b) The PPI network of the protein encoded by Os06g0227200. Green, red, and blue lines indicate interactions predicted from gene neighborhood, fusions, and co-occurrence. Light green, black, and dark blue lines indicate additional interactions inferred from text mining, co-expression, and protein homology. (c) Expression profile of annotated transcript and novel isoforms produced by Os06g0227200; Figure S3. Comparison of rice spliceosome component genes that are significantly differentially expressed at (a) 24

Conflicts of Interest:
The authors declare no conflict of interest.