Local and systemic transcriptome and spliceome reprogramming induced by the root-knot nematode Meloidogyne incognita in tomato

Abstract Root-knot nematodes (Meloidogyne spp.) are widely spread root parasites that infect thousands of vascular plant species. These highly polyphagous nematodes engage in sophisticated interactions with host plants that results in the formation of knot-like structures known as galls whose ontogeny remains largely unknown. Here, we determined transcriptome changes and alternative splicing variants induced by Megalaima incognita in galls and neighboring root cells at two distinct infective stages. M. incognita induced substantial transcriptome changes in tomato roots both locally in galls and systemically in neighboring cells. A considerable parallel regulation of gene expression in galls and neighboring cells were detected, indicative of effective intercellular communications exemplified by suppression of basal defense responses particularly during the early stage of infection. The transcriptome analysis also revealed that M. incognita exerts a tight control over the cell cycle process as a whole that results in an increase of ploidy levels in the feeding sites and accelerated mitotic activity of the gall cells. Alternative splicing analysis indicated that M. incognita significantly modulates pre-mRNA splicing as a total of 9064 differentially spliced events from 2898 genes were identified where intron retention and exon skipping events were largely suppressed. Furthermore, a number of differentially spliced events were functionally validated using transgenic hairy root system and found to impact gall formation and nematode egg mass production. Together, our data provide unprecedented insights into the transcriptome and spliceome reprogramming induced by M. incognita in tomato with respect to gall ontogeny and nematode parasitism.


Introduction
Root-knot nematodes (Meloidogyne spp.) are obligatory root parasites that infect almost all vascular plant species including many economically important crops, vegetables, and ornamental plants, causing significant economic losses [1].With over 100 species of Meloidogyne species described so far [2], Megalaima incognita is the most globally widespread and damaging species.These highly polyphagous nematodes engage in sophisticated interactions with host plants that results in the formation of a permanent feeding site, known as giant-cells, to support nematode development and maturation.The infective second-stage juveniles (J2) penetrate plant roots at the root tip region, migrate intercellularly, and select few vascular cylinder cells as potential feeding cells.If these cells are compatible, the infective J2s become sedentary and guide these cells to enlarge in size via recurrent mitosis without cytokinesis to form multinucleated and metabolically hyper-active giant-cells [3,4].Giant-cell formation stimulates neighboring cells to divide asymmetrically in a disordered way, resulting in the formation of knot-like swellings termed galls at the site of infection.The infective J2s feed from the giant-cells and molts to third-stage juveniles (J3), fourth-stage juveniles (J4), and finally adult females laying hundreds of eggs inside the egg [1,3].
The identity of nematode-induced galls and giant-cells is believed to be established through specific transcriptional reprogramming, which is mediated via a broad set of transcription factors in addition to other transcriptional and posttranscriptional regulators.M. incognita-induced transcriptome programming in host plants has been documented in several studies [5][6][7][8][9][10][11][12][13][14], providing interesting information about the importance of phytohormone signaling, suppression of defense responses, cytoskeletal organization, solute transport, cell wall modifications, regulation of metabolic pathways, for successful nematode parasitism of host plants.However, most of these transcriptome studies were conducted using whole root tissues, and hence lack the fine resolution required for deepening our understanding of cellular responses to M. incognita infection and the molecular mechanisms underlying gall formation and development.
Post-transcriptional regulation of mRNA by alternative splicing (AS) has recently emerged as an important regulatory mechanism determining gene function and specificity in various developmental and stress contexts [15][16][17][18][19]. Exon skipping and intron retention, the two main splicing events, impact transcriptome diversity and function.For example, the generation of multiple mature mRNA transcripts from a single pre-mRNA through combinatorial exclusion or inclusion of exons has the potential to produce protein isoforms with varied structural and functional specificities.Intron-retention events frequently produce non-functional transcripts that could play a role in fine-tuning the transcript level of certain disease resistance genes to reach appropriate activity and maximize plant fitness [19][20][21].In some cases, intron-retained transcripts can be translated, producing truncated proteins that may have dominant-negative functions [18,20,22].
Despite recent experimental data indicate that AS plays fundamental role in various pathosystems [17,19], information about the frequency and importance of AS events during plantnematode interactions has been reported in only a limited number of studies [23][24][25][26].For example, soybean cyst nematode (SCN, Heterodera glycines) has been shown to induce 1979 splicing events in 1576 unique genes in soybean roots upon infection with exon skipping being the most common differentially spliced events [26].Importantly, effectors from cyst and root-knot nematodes have been shown to interacts with host spliceosomal proteins in order to reprogram plant mRNA splicing machinery and facilitate infection.The 30D08 effector from the beet cyst nematode Heterodera schachtii and the MiEFF18 effector from M. incognita were reported to physically associate with the auxiliary spliceosomal SMU2 and SmD1 from Arabidopsis, respectively, thereby modulating the splicing patterns of plant genes to facilitate nematode infection [23,24].Despite these compelling evidences, the magnitude, patterns, and functions of differentially spliced events occurring during the formation and development of M. incognita-induced galls remain elusive.
In this study, we determined transcriptome changes and alternative splicing variants, induced by M. incognita in galls and neighboring root cells at two distinct infection stages.Our data indicate that M. incognita induces robust gene expression changes in tomato roots both locally in galls and systemically in neighboring cells.This transcriptome reprograming occurred in an infection stage-dependent manner and involved parallel regulation in galls and neighboring cells.Our analysis also revealed that M. incognita greatly modulates pre-mRNA splicing particularly during the early stage of infection as a total of 9064 uniquely differentially spliced events from 2898 genes were identified.The functional importance of certain splicing variants was examined using transgenic hairy root system and found to alter nematode parasitism by impacting gall formation and nematode egg mass production.Overall, our results shed light on how M. incognita alters host gene expression and pre-mRNA splicing to regulate various biological processes and molecular functions required for gall organogenesis and nematode parasitism.

M. incognita induces substantial transcriptome changes in tomato roots both locally in galls and systemically in neighboring cells
To determine local and systemic transcriptome changes induced by M. incognita in tomato roots, we inoculated 10-day-old tomato seedings with second-stages juveniles.At 4-and 11-day post inoculation (dpi), galls, neighboring roots tissues, and the corresponding non-infected whole roots we collected in three biological samples for RNA extraction and RNA-seq library preparation (Supplementary Fig. S1).The 4-and 11-day time points were selected to represent young and fully developed galls.Using a false discovery rate less than 5%, we identified 3895 and 6944 differentially expressed genes (DEGs) in galls and adjacent root tissues, respectively, in comparison with non-infected roots at 4 dpi (Supplementary Data Table S1 and S2).Similarly, at 11 dpi 3582 and 1063 DEGs were identified in galls and adjacent root tissues, respectively, in comparison with the corresponding noninfected roots (Supplementary Data Table S3 and S4).These data indicate that M. incognita induces robust gene expression changes in tomato roots both locally in galls and systemically in neighboring cells.

M. incognita induces transcriptome reprograming in an infection stage-dependent manner
We examined whether M. incognita induces distinct transcriptome reprograming at the 4 and 11 dpi.Venn diagrams revealed that 801 and 1443 genes were uniquely differentially expressed in galls at the 4 and 11 dpi time points, respectively (Fig. 1A).Likewise, 4062 and 74 genes were exclusively differentially expressed in the neighboring root tissues at the 4 and 11 dpi time points, respectively (Fig. 1A).This finding implies that M. incognita induces transcriptome reprograming in galls and neighboring cells in an infection stage-dependent fashion.The 4-dpi gall specific DEGs (801) were enriched in genes involved in small molecule biosynthetic process, whereas the 11-dpi gall-specific DEGs (1443) were enriched in genes implicated in mitotic cell cycle process and microtubule-based process.The DEGs unique to the neighboring cells at 4 dpi (4062) were enriched in genes implicated in cytoplasmic translation, cellular catabolic process, vesicle-mediated transport, phosphorus metabolic process, protein folding, targeting, and modification, mRNA splicing, organelle organization, cell communication, and cellular response to stimulus (Fig. 1B).The DEGs unique to the neighboring cells at 11 dpi (74) did not show enrichment in a particular biological process (Fig. 1B).

M. incognita induces parallel transcriptome reprograming in galls and neighboring cells
Cross comparison of the identified DEGs revealed that 2470 (2066 + 404) and 674 (517 + 157) genes were similarly up-or downregulated in galls and adjacent root tissues at 4 and 11 dpi, respectively (Fig. 2 A and B).Only three genes at these two time points showed opposite regulation in galls and neighboring root tissues (Fig. 2A and B).The finding that 63.5% of the DEGs identified in galls at 4 dpi were similarly regulated in neighboring root tissues implies that this parallel regulation is not random but is biologically meaningful.Interestingly, gene ontology (GO) analysis of the 2066 downregulated genes both in galls and in neighboring cells at 4 dpi were enriched in genes involved in defense response to other organisms, as well as those implicated in cell wall organization or biogenesis and transmembrane transport (Fig. 2C).Protein class analysis revealed that water and metal ion transmembrane transporters and peroxidases were highly overrepresented among these downregulated genes.The similarly upregulated genes (404) in galls and in neighboring root tissues at 4 dpi were enriched in genes involved in wounding response, regulation of cell communication and regulation of gene expression (Fig. 2C).Remarkably, protein class analysis of these 404 upregulated genes revealed enrichment for transcription factors and metabolite interconversion enzymes, particularly those with hydrolase, oxidoreductase, phospholipase, and oxygenase functions.The 674 genes that were similarly regulated in galls and adjacent roots tissues at 11 dpi included 517 downregulated genes and 157 upregulated genes.Response to cytokinin and response to auxin were the significantly enriched GO terms among the 517 downregulated genes (Fig. 2C).The 157 similarly upregulated genes in galls and neighboring regions were involved in a wide range of biological processes but no such enrichment for GO

Transcription factors regulatory networks controlling gene expression in galls and neighboring regions
The identification of substantial numbers of commonly and uniquely regulated genes between galls and adjacent roots tissues at both time points prompted us to identify transcription factors with potential roles in this regulation.A total of 261 and 319 transcription factor-encoding genes were identified as differentially expressed in galls and adjacent root tissues, respectively at 4 dpi (Supplementary Data Table S5).Similarly, at 11 dpi, 212 and 90 transcription factors were determined as differentially expressed in galls and adjacent root tissues, respectively (Supplementary Data Table S5).After eliminating duplicates, a total of 487 unique transcription factors were identified as differentially expressed in galls and/or adjacent root cells at both time points.As shown in Fig. 3, 236 of these transcription factors belongs to the bHLH, bZIP, ERF, MYB, and WRKY families, implying key functions of these family members in establishing the regulatory networks controlling gene expression in galls and neighboring regions.
Cross-comparison of the transcription factor gene lists revealed that 162 and 62 transcription factors (highlighted in back dots in Fig. 3) were similarly regulated in galls and adjacent root tissues at the 4 and 11 dpi time points, respectively.These factors belong to various families, including AP2, bHLH, bZIP, GATA, ERF, TCP, TGA, and Trihelix, for example, suggesting roles of these factors in coordinating gene expression in galls and adjacent root cells.Cross comparison also revealed that 256 and 178 transcription factors (highlighted in red dots in Fig. 3) were uniquely regulated in galls or neighboring root tissues at 4 and 11 dpi, respectively.Notably, all members of the E2F, GTE, and HD transcription factor families exhibited unique expression patterns in galls or neighboring root tissues (Fig. 3), suggesting key regulatory functions of these factors in mediating gene expression specificity.Examining the expression patterns of the 487 identified transcription factors revealed that two-thirds of these transcription factors were downregulated (Fig. 3 and Supplementary Data Table S5).Remarkably, all members belonging to the TCP, TGA, HD-like, and NF-YA families were downregulated (Fig. 3).Similarly, the majority of the bZIP, GATA, and GRAS family members were also downregulated (Fig. 3).Together, these data suggest that downregulation of key regulatory factors is critical for gall formation and development.

Mitotic cell cycle processes are highly activated in fully developed galls
The enrichment of mitotic cell cycle-related genes among the 11dpi gall-specific DEGs is remarkable considering the key role of cell cycle in gall ontogeny and nematode parasitism.Therefore, we mapped all differentially expressed cell cycle-related genes, identified in the 11-dpi galls in comparison with the corresponding non-infected roots, across various cell cycle phases.As shown in fig.4, these genes were distributed across all cell cycle phases including interface phase (G1, S, and G2) and mitotic (M) phase (mitosis and cytokinesis).

Local and systemic suppression of defense response in galls and neighboring root cells
GO term analysis implied local and systemic suppression of defense responses in galls and neighboring root cells during the early stage of infection.To provide more details for this finding, we compared the levels of downregulation of a set of genes known to be involved in basal defense response between galls and neighboring root cells at both time points.The expression of 11 defensin-like genes were simultaneously downregulated in the 4-dpi samples.However, the level of downregulation in neighboring root cells tended to be at a lower magnitude (Fig. 5).
Additionally, genes encoding pathogenesis-related (PR) proteins, respiratory burst oxidases, and lipoxygenase exhibited the same patterns of downregulation (Fig. 5), supporting the notion of systemic suppression of defense responses upon nematode infection.Nevertheless, at the 11-dpi time point, the majority of these genes showed similar downregulation in galls but not in the neighboring root cells, suggesting a transition of defense suppression from local/systemic state at the early stage of infection into local state at later stage on infection.
It may be important to mention that a number of defenserelated genes, including PR genes, defensin-like genes, receptor  S1  and S2).The upregulation of these defense-related genes may ref lect that tomato plants ineffectively respond to nematode infection or ref lect a mechanism to mitigate damage.

M. incognita profoundly modulates pre-mRNA splicing during infection
We next analyzed our RNA-seq data to identify differentially spliced events using the JunctionSeq package [27].Using a FDR cut-off of 0.05, we identified 2272 and 6580 differentially spliced events in galls and neighboring root tissues, respectively, as compared with non-infected roots at 4 dpi (Supplementary Data Tables S7 and S8, and Fig. 6A and B).Similarly, 933, and 1603 differentially spliced events were identified in galls and adjacent root tissues, respectively, compared with non-infected roots at 11 dpi (Supplementary Data Tables S9   and S10, and Fig. 6C and D).A total of 9064 unique differentially spliced events from 2898 genes were identified after eliminating duplicates.These results indicate that M. incognita profoundly modulates pre-mRNA splicing particularly during the early stage of infection.

M. incognita induces suppression of intron retention events and exon usage both in galls and neighboring cells
The differentially spliced events were further grouped into two main categories: intron retention (IR) and differentially used exons (DUEs).Of the 2272 differentially spliced events identified in the 4-dpi galls, 740 events (402 IR and 338 DUE) were significantly more abundant in galls, and 1532 events (962 IR and 570 DUE) were significantly more abundant in the corresponding non-infected root tissues (Fig. 6A).Of the 6580 differentially spliced events identified in the neighboring regions at 4 dpi, 1134 events (524 IR and 610 DUE) were significantly more abundant in neighboring  S6.regions, and 5446 events (3319 IR and 2127 DUE) were significantly more abundant in the corresponding non-infected root tissues (Fig. 6B).These data imply that suppression of IR events and increasing exon skipping occur in galls and to a much higher extent in the neighboring cells.
This trend was also observed in galls and neighboring cells at 11 dpi despite the fact that the numbers of differentially spliced events were substantially reduced.The 933 differentially spliced events identified in galls at 11 dpi included 337 events (180 IR and 157 DUE), which accumulated more significantly in galls, and 596 events (377 IR and 219 DUE), which accumulated more significantly in the corresponding non-infected root tissues (Fig. 6C).The 1603 differentially spliced events identified in the neighboring regions at 11 dpi included 965 events (562 IR and 403 DUE), which accumulated more significantly in neighboring regions, and 638 events (307 IR and 331 DUE), which accumulated more significantly in the corresponding non-infected root tissues (Fig. 6D).
To validate IR and DUE data, we used reverse-transcription quantitative PCR (RT-qPCR) to measure the abundance of 6 IR events (Solyc08g067360, Solyc07g009070, Solyc10g084780, Solyc01g080150, Solyc04g009310, and Solyc01g095140) and 6 exon skipping events (Solyc07g054410, Solyc08g082130, Solyc01g057420, Solyc01g008090, Solyc10g048030, and Solyc09g011500) in galls or neighboring cells at 4 or 11 dpi relative to non-infected root tissues.With the exception of the Solyc01g008090 exon skipping event, which did not show significant change in abundance, the expression data of the remaining 11 events confirmed the splicing patterns obtained with JunctionSeq (Fig. 6E).

Genes involved in oxidation reduction are subjected to alternative splicing in gall and neighboring cells
We next performed GO term enrichment analysis to identify biological processes and molecular functions that are overrepresented among the differentially spliced genes (DSGs) showing increased or decreased IR/exon usage events in galls and neighboring cells at 4 and 11 dpi.At 4 dpi, genes showing increased IR and/or exon usage in galls were enriched in functions related to catalytic activity and wounding response, whereas those showing reduced IR and/or exon usage were enriched in genes related to oxidation-reduction process and transmembrane transporter activity.Similarly, genes showing increased IR and/or exon usage in the neighboring cells at 4 dpi were enriched in functions related to catalytic activity, hydrolase activity, and wounding response.However, genes with reduced abundance of IR and/or exon usage events in the neighboring region at 4 dpi were enriched in biological processes and molecular functions related to cellular lipid metabolic processes, protein and small molecule metabolic processes, macromolecule biosynthetic process, transport, oxidoreductase activity, and RNA binding.At 11 dpi, DSGs with increased IR and/or exon usage in galls showed an enrichment for oxidoreductase activity, while those exhibiting reduced IR and/or exon usage were enriched in oxidoreductase activity and transmembrane transporter activity.DSGs with increased IR and/or exon usage events in the neighboring region at 11 dpi were enriched in functions related to wounding responses and oxidoreductase activity, whereas those genes showing reduced IR and/or exon usage events were enriched in genes related to catalytic activity and primary metabolic process (Supplementary Fig. S2).These data clearly show that a significant number of genes involved in oxidation-reduction are subjected to alternative splicing both in gall and neighboring cells at early and late stages of nematode infection.

Differentially used exons impact protein domain arrangement
Differential usage of exons in specific splicing variants have the potential to impact protein domain organization, and hence protein function [18].Therefore, we examined if DUEs identified in galls and neighboring cells would affect protein domain structure by scanning the amino acid sequences of significantly DUEs against functionally annotated protein domains in the Pfam database.Numerous functional domains, including cytochrome P450, glycoside hydrolase, peroxidase, peptidase, helicase, Sadenosyl-l-methionine-dependent methyltransferase, and small GTPase were identified among DUEs (Fig. 7), suggesting that DUEs have the potential to alter the activity of the encoding enzymes.Domains functioning in immunity and defense responses were also identified, including protein kinase, lipoxygenase, MLO, and Tify (Fig. 7).

Differentially used exons impact gall formation and egg mass production
We examined the functional impact of exon splicing events on plant response to nematode infection by overexpressing splicing variants lacking DUEs using transgenic hairy root system.The overexpression constructs included the coding sequences of Solyc03g094160 lacking exon 4, Solyc12g013650 lacking exon 2, Solyc01g079260.4.1 lacking exon 3, and Solyc11g032100.2lacking exon 2. In all cases, full-length coding sequences of these genes were also overexpressed for comparison.In addition, transgenic hairy roots expressing the empty vector were generated and included as control.Transgenic hairy roots were identified using GFP screens, and the composite plants were then inoculated with M. incognita second-stage juveniles.The functionality of the overexpression constructs was confirmed by quantifying mRNA abundance of the overexpressed genes in three biological samples randomly collected from different GFP-positive hairy roots.Approximately, 6 weeks after inoculation number of galls per root system as well as number of eggs per gram root tissues were determined.Overexpression of both full-length and splicing variants of Solyc03g094160 negatively impacted the number of galls as compared with control plants expressing the empty vector (Fig. 8A).However, the splicing variant produced significantly lower number of galls as compared with the full-length variant (Fig. 8A).Overexpression of the full-length variant resulted in 54.4% reduction in number of galls, whereas the splicing variant showed 73% reduction in comparison with control plants (Fig. 8A).
Both variants showed similar impact on egg number (Fig. 8C).
In comparison with control plants expressing the empty vector, both full-length and the splicing variant of Solyc12g013650 sig-nificantly reduced the number of galls but the full-length variant produced significantly less galls compared with the splicing variant (Fig. 8A).Both variants showed similar non-significant reduction in egg numbers relative to control plants expressing the empty vector (Fig. 8C).Transgenic hairy roots overexpressing the full-length and splicing variants of Solyc01g079260 exhibited distinct impacts on gall and egg numbers when compared with control plants (Fig. 8B and D).While the full-length variant produced significantly less gall, the splicing variant produced gall number similar to that of control plant (Fig. 8B).Only the splicing variant significantly reduced the number of eggs as compared with control plants (Fig. 8D).Both full-length and splicing variants of Solyc11g032100.2produced similar number of galls per root system, but the number of eggs determined on the transgenic roots overexpressing the splicing variant were twice the number determined on the transgenic roots overexpressing the full-length variant (Fig. 8B and D).Taken together, these data support a role of DUEs in gall formation and/or egg mass production.

Discussion
In this study, we investigated local and systemic transcriptome and spliceome reprogramming in the M. incognita-induced galls and neighboring root cells, and revealed several key features of this reprogramming.M. incognita induced substantial transcriptome changes in tomato roots both locally in galls and systemically in neighboring cells at early and late stages of infection.Previous studies documenting transcriptome changes using whole roots as experimental materials in response to M. incognita infection suggest a systemic induction of gene expression changes [5][6][7][8][9][10][11][12][13][14].However, the characteristics and magnitude of this induction were not previously reported.Our finding that 63.5% of the DEGs identified in galls at 4 dpi were similarly regulated in neighboring root tissues indicates that this systemic regulation is of biological significance particularly at the early stage of infection.However, at the later stages of infection, the parallel regulation of gene expression in gall and neighboring root cells seems to be less profound as only 19% of the DEGs identified in galls at 11 dpi were similarly regulated in neighboring cells.Cell-to-cell communication during plant development are generally mediated through various signaling mechanisms involving reactive oxygen species (ROS), small peptides, phytohormone signaling and the movement of transcription factors and small non-coding RNAs [28,29].The involvement of these signaling mechanisms in establishing plantnematode interactions [30][31][32][33] supports such a role in effective intercellular communication between galls and neighboring root cells considering a significant number of genes coding for these signaling components are identified as differentially expressed.
Coordinated intercellular communications between plant cells is required for proper plant responses to infection relating to defense responses and resource allocation [34].In this regard, our analysis provided new insights into the coordinated cell-tocell communications between the nematode-induced galls and neighboring root cells.We observed downregulation of a substantial number of basal defense response genes both in galls and neighboring root cells at the 4-dpi time point (Fig. 5).However, at 11 dpi, the downregulation of the large majority of these genes was restricted to gall cells, suggesting a transition of defense suppression from local/systemic state into local state.This hypothesized transition is likely the result of a switch in plant communication system prioritizing the production and transport of nutrients to the nematode feeding cells.Our analysis also provided evidence involving jasmonic acid signaling in the coordinated cell-to-cell communications between galls and neighboring root Figure 8. Impact of exon skipping events on plant response to nematode infection.The differentially spliced variants of Solyc03g094160 (lacking exon 4), Solyc12g013650 (lacking exon 2), Solyc01g079260 (lacking exon 3), and Solyc11g032100 (lacking exon 2) were overexpressed in the susceptible tomato cultivar Heinz 1706 using transgenic hairy root system.Composite plants with GFP positive roots were inoculated with approximately 500 second-stage juveniles of Megalaima incognita.Six weeks after inoculation, the number of galls per plant (A and B) and number of eggs per gram of root tissues (C and D) were determined.Full-length coding sequences of these four genes were also overexpressed for comparison.Data are presented as mean ± SE.Each dot represents one data point.Statistically significant differences between overexpression and control plants expressing the empty vector were calculated using ANOVA ( * P < 0.05).cells.Notably, auxin and cytokinin signaling genes were enriched among the similarly downregulated genes in galls and neighboring cells at 4 dpi, whereas genes involved in jasmonic acidmediated signaling pathway were enriched among the similarly upregulated genes, providing suggestive evidence for suppression of auxin and cytokinin signaling in these tissues.This suggestion is supported by our finding that several auxin response factors (ARFs) and two-component response regulators (type-B ARRs), which are involved in auxin-and cytokinin-mediated signaling [35,36], respectively, were downregulated in galls and neighboring cells at 4 dpi.In contrast, genes encoding auxin/indole-3acetic acid proteins (Aux/IAAs) and cytokinin oxidases, which inhibit auxin and cytokine signaling, respectively [35,37,38], were upregulated.Unlike auxin and cytokinin, the activation of numerous genes involved in jasmonic acid biosynthesis and signaling suggests a role of jasmonic acid and its metabolite methyl jasmonate, produced upon nematode infection, in systemic signal transmission and cell-to-cell communication.The function of JA and MeJA in short-and long-distance signal transmission has been well established in various plant species, including tomato [39,40].In this regard, the systemic reprograming of cells adjacent to the infection sites at early infection stage seems to be critical for creating favorable environments for successful nematode infection by facilitating nutrient f low to the giant-cells.The enrichment of genes involved in translation, cellular catabolic process, vesicle-mediated transport, protein folding, targeting, and degradation, organelle organization, and cell communication supports this view.
Plant transcription factors play key regulatory functions controlling various developmental and physiological processes that determine cell identity [41].Here, we cataloged about 800 differentially expressed transcription factors belonging to various families, expanding the limited list of regulatory factors functioning during M. incognita -tomato interactions.Member of the bHLH, bZIP, GATA, ERF, MYB, WRKY, and ZF families constituted a significant portion of these factors.While members of these families have been shown to regulate plant-nematode interactions [42][43][44][45][46][47][48].Combinatorial interactions between TFs belonging to the same or different families may determine the specificity of the regulatory networks required for gall formation and development.Consistent with this notion, various members of the bHLH, MYB, and bZIP families have been found to form homo-and heterodimers with proteins from the same or other families that impact their regulatory function and specificity [42,[49][50][51].Our analysis also suggests that members of the transcription factor families E2F, GTE, and HD may contribute to gene expression specificity because they exhibited unique expression patterns in galls or neighboring root tissues.Our analysis also revealed that all differentially expressed transcription factors belonging to the TGA family were downregulated.This may be related to the function of these factors in regulating defense-related genes through physical association with NPR1 (NON-EXPRESSOR OF PR1), a key regulator of immunity [52].These factors also function in plant immunity in an NPR1-independent manner [53,54].Notably, several Trihelixand Scarecrow-like transcription factors exhibited distinct and overlapping expression patterns in galls and neighboring cells, suggesting that spatial and temporal activities of these factors may contribute to nematode infection success.The potential functional importance of these factors is supported by the findings that members of these families were targeted by M. incognita effectors [55,56].
Unlike the TGA family members, we found that all differentially expressed transcription factors of the E2F family were upregulated, indicative of cell cycle activation during nematode parasitism, considering the fundamental functions of these factors in regulating genes involved in cell cycle progression and DNA synthesis [57].This indication was further supported by the enrichment of genes involved in mitotic cell cycle among the 11-dpi gall specific DEGs in which 66 out of 72 genes were upregulated.This is consistent with previous findings showing activation of several core cell cycle genes in nematode inducedfeeding sites and galls [58][59][60][61][62][63][64].The functional distribution of these genes across various cell cycle stages and checkpoints (Fig. 4) shows how M. incognita exerts tight control over the cell cycle process as a whole during infection that results in an increase of ploidy levels in giant-cells and accelerated mitotic activity of surrounding cells.Interestingly, among the cell cycle-related genes we identified several involved in cell cycle checkpoints including, for example, cyclin-dependent kinases, AURORA kinases, and the mitotic checkpoint proteins BUBR1, BUB1 and BUB3.The abundance and differential regulation of these genes may play roles in switching ordinary mitotic cell cycle in giant-cells to aberrant cell cycle involving acytokinetic mitosis and endoreduplication to increase nuclear ploidy levels required for giant-cell ontogeny and enlargement [65,66].
Our analysis of AS events in galls and adjacent root tissues revealed key features of spliceome reprograming during M. incognita infection.First, spliceome reprograming induced by M. incognita is profound accounting for 9064 events from 2898 genes.Second, M. incognita modulates pre-mRNA splicing at 4 dpi to a much higher degree compared to the 11-dpi time point (8852 events vs 2536).This is consistent with the concept that AS is dynamic in actively dividing cells and differentiating tissues [18,67,68].However, we found that AS events in adjacent root tissues were more prevalent than in galls (8183 events vs 3205).This indicates that cells adjacent to M. incognita-induced galls undergo massive posttranscriptional programming involving notable AS, highlighting the essential roles of these cells in establishing the compatibility of the interactions.Third, we found that 52% (1179 events) of the AS events detected in galls at 4 dpi were also identified in the adjacent root cells but this overlap was diminished to only 16% at the 11-dpi time point.This finding indicates that, similar the transcriptome programming, AS in galls and adjacent root tissues are highly coordinated specifically at the early stages of infection despite the mechanism underlying this coordination remain to be determined.Fourth, M. incognita-induced AS appears to be linked to gene expression levels particularly in galls.For example, 58.7% and 54% of the DSGs identified in galls at 4 and 11 dpi, respectively, were differentially expressed.In adjacent root cell, this percentage dropped to 48.4% and 14.5% for the 4 and 11 dpi, respectively.A similar association has been previously reported [18].The substantial overlaps between DSGs and DEGs particularly in galls indicate that transcriptome and spliceome reprograming driven by M. incognita infection is mechanistically linked.This can be explained by the functional association between the RNA polymerase II (pol II) and splicing factors where pol II elongation rates regulate AS, and splicing factors can reciprocally alter pol II elongation [69].This association can also be mediated through a set of co-regulated transcription factors and splicing factors or through regulatory factors with dual functions in gene transcription and splicing.
The splicing patterns revealed suppression of IR and exon usage events in galls and to a much higher degree in the neighboring cells at 4 dpi.However, at 11 dpi, this pattern was reversed in the neighboring cells with IR and exon usage events being more abundant as compared with the corresponding non-infected root tissues.These data suggest that suppression of IRs facilitates M. incognita infection of tomato.This can be explained by the fact that IR frequently produces premature termination codons [22], and suppression of these events in galls and neighboring cells would reduce the accumulation of non-functional transcripts, thereby enhancing the translation efficiency of functional transcripts, and hence plant susceptibility.In accord with this view, it has been shown that IR events accumulate to a much higher magnitude in the resistant interactions compared with susceptible interactions [70,71].Our analysis also suggests that DUEs, which lead to exon skipping events, may contribute to the establishment of plant-nematode interactions.Protein domain analysis revealed that differential usage of exons can alter protein domain arrangement, and subsequently the activity of various enzymes, including protein kinase, lipoxygenase, and peroxidase, for instance.The enrichment of the peroxidase domain among the DUEs is remarkable and suggests that AS of peroxidase-encoding genes plays a central role in the regulation of plant cellular responses to oxidative stress during infection.The importance of exon usage for the establishment of plant-nematode interactions was further validated by dur data showing that overexpression of 4 splicing variants lacking DUEs negatively impacted gall formation and/or egg mass production.These splicing variants may produce protein isoforms with distinct structures, biological functions, cellular localizations, and protein interaction networks that interfere with successful nematode parasitism.
In conclusion, our study has revealed extensive transcriptome and spliceome reprogramming that occurs during the compatible interaction between M. incognita and tomato plants, and has pointed to possible mechanisms of action that facilitate this interface.Our analyses also indicate that widespread modulations of gene expression and AS events occurring systemically in cells adjacent to the M. incognita-induced galls represent important features of this interaction.Together, our data documented in the current study provide the foundation for future studies to investigate the importance of local and systemic changes in gene expression and AS with respect to gall ontogeny and nematode parasitism.

Plant material and growth conditions
Tomato (Solanum lycopersicum) cultivar Heinz 1706 was used in all experiments described in this study.Plants were grown in controlled growth chambers at 26 • C under 16-h light/8-h dark conditions.

Nematode inoculation, tissue collection, and RNA-seq library construction
Seeds were surface-sterilized and germinated in sterile topsoil (Promix BK25 Peat/Bark Based).Then, 2-week-old tomato seedlings were transferred to pots containing a sterilized mixture of sand and topsoil at a 3:1 ratio.Each plant was inoculated with about 300 second-stage juveniles of M. incognita or a mock 0.1% agarose solution (control).Inoculated and control plants were organized in a randomized complete block design in a growth chamber at 26 • C under 16-h light/8-h dark conditions and 65% humidity.Galls, neighboring regions, and non-infected whole roots were collected at 4-and 11-day post inoculation (dpi) in three biological replicates, resulting in a total of 18 samples.Approximately, 250 galls and neighboring regions were collected for each replicate.Isolation of DNA-free RNA was performed using Direct-zol™ RNA Miniprep Plus Kit (Zymo Research) following the manufacturer's protocol.rRNA-depleted RNA-seq libraries were prepared using rRNA-depleted RNA and a NEBNext Ultra Directional RNA Library Prep Kit for Illumina (New England Biolabs, Beverly, MA, USA) according to the manufacturer's protocol.The libraries were sequenced using the NovaSeq 6000 platform (Illumina Corp., San Diego, CA, USA) with 150-bp paired-end reads.

Alternative splicing analysis
Differentially used exons (DUE) and intron retention (IR) events were determined using the R package JUNCTIONSEQ [27], which utilizes QORTS quality-control/data-processing software [76] for read quality control counts.Reads mapped to each nonoverlapping exons and splice junctions were counted and used to determine DUE and IR events using DESEQ2 [75] with a false discovery rate (FDR) cut off of 5% as recently described [18].

Gene ontology enrichment analysis
Gene Ontology (GO) term enrichment analysis of DEGs and DSGs was performed using PANTHER version 18.0 with Fisher's exact test and Bonferroni multi-test correction with a cut off P-value of 0.01 for significance.

Generation of transgenic tomato hairy roots and nematode infection assays
The full-length coding sequences as well as differentially spliced variants of Solyc01g079260 (lacking exon 3), Solyc11g032100 (lacking exon 2), Solyc12g013650 (lacking exon 2), and Solyc03g094160 (lacking exon 4) were synthesized and cloned in the pG2RNAi2 vector under the control of a strong ubiquitin promoter.Primer sequences used for cloning are provided in Supplementary Table S11.The pG2RNAi2 binary vector contains a super folded green f luorescent protein (sGFP) for the identification of transgenic hairy roots.All constructs were confirmed by sequencing, and then were transformed into Agrobacterium -rhizogenes strain K599 and used for generating transgenic tomato hairy roots in the susceptible Heinz 1706 background.The roots of 1-week-old seedling were removed and the seedlings were vacuum-infiltrated with agrobacterium suspension (OD600 of 1.0) four times each for 90 seconds.The infiltrated seedlings were maintained under high humidity conditions (65%) in a growth chamber with 16h light/8-h dark conditions at 26 • C for 3 weeks.The generated hairy roots were then screened using dissecting microscope containing GFP filter to identify transgenic hairy roots.Composite plants with GFP-positive roots were planted in pots containing topsoil (Premier Pro-Mix BK25 Peat/Bark Based Growing Soil), and inoculated with approximately 500 freshly hatched second-stage juveniles of M. incognita.Six weeks after inoculation, the number of galls per plant and number of eggs per gram of root tissues were determined.The infection experiments were organized in a randomized complete block design with between 15 and 20 replicated plants per construct.Statistically significant differences between overexpression and control plants expressing the empty vectors were calculated using analysis of variance with P < 0.05.

Reverse-transcription quantitative PCR assays
To confirm the functionality of overexpression constructs, total RNA was extracted from transgenic hairy roots overexpressing the full-length coding sequences and differentially spliced variants of Solyc01g079260, Solyc11g032100, Solyc12g013650, and Solyc03g094160 as previously described [77] and subjected to RT-qPCR quantification using gene specific primers (Supplementary Data Table S11).Total RNA was also extracted from transgenic hairy roots expressing the empty binary vector and used as control.To quantify the expression levels of various splicing events, total RNA was isolated from galls and neighboring root regions at 4 and 11 dpi.Total RNA was also isolated from noninfected root tissues at each both time points and used as control.Primers were designed to specifically amplify the directed intron or exon in each splicing event (Supplementary Data Table S11).RT-qPCR quantifications were conducted using Verso SYBR green 1-step RT-qPCR (Thermo Fisher Scientific), with 20 ng of RNA per reaction, following the manufacturer's instruction.Gene expression values were normalized using actin (Solyc11g005330.2.1) as a reference gene.Relative fold change values were calculated using three biological samples.Each reaction was performed with three technical replicates.Statistically significant differences between control and overexpression samples were determined by t tests (P < 0.05).

Figure 1 .
Figure 1.Megalaima incognita induces transcriptome reprograming in an infection stage-dependent manner.A: Venn diagram showing overlapping and uniquely differentially expressed genes (DEGs) between galls and neighboring root cells (NR) at 4 and 11 dpi.B: Heat map representation of GO terms significantly enriched among DEGs unique to galls and neighboring root cells at 4 and 11 dpi.GO term enrichment analysis was performed using PANTHER version 18.0 with Fisher's exact test and Bonferroni multi-test adjustment with a cutoff P-value of 0.01 for significance.

Figure 3 .
Figure 3. Classification of differentially expressed transcription factors identified in galls and adjacent roots cells at 4 and 11 dpi.Shown are fold change values of differentially expressed transcription factors belonging to various families/subfamilies.Each transcription factor is represented by a dot.Black dots denote transcription factors that were similarly regulated in galls and adjacent roots cells at the 4 or 11 dpi time points.Red dots denote transcription factors that were specifically regulated in galls or adjacent roots cells at 4 or 11 dpi.

Figure 4 .
Figure 4. Schematic representation showing distribution of differentially expressed cell cycle-related genes identified in the 11-dpi galls across various cell cycle phases.Upregulated genes are indicated in black, whereas downregulated genes are indicated in red.Gene expression values and annotations are provided in Supplementary Data TableS6.

Figure 5 .
Figure 5. Local and systemic suppression of defense response in galls and neighboring root cells.Fold change plot of downregulated defense-related genes in galls and neighboring root cells at 4 and 11 dpi.Bubble size represents the level of downregulation for each gene.

Figure 6 .
Figure 6.Modulation of pre-mRNA splicing during Megalaima incognita infection.A-D: MA plots comparing fold change values against the mean of the normalized counts of differentially spliced events determined in galls (A and C) and neighboring region (NR) (B and D) at 4 and 11 dpi.Red dots represent significantly differentially spliced events.DUE: differentially used exons; IR: intron retention.E: RT-qPCR quantification of the expression levels of differentially spliced events in galls and NR at 4 and 11 dpi.Shown are the fold change values of 6 exon skipping (ES) and 6 IR events in galls or NR at 4 and 11 dpi relative to the corresponding non-infected control root tissues.Data are the mean ± SE of three biological replicates.

Figure 7 .
Figure 7. Impact of differentially used exons on protein domain organization.The bubble plot illustrates the most common protein domains identified in the differentially used exons in galls and neighboring root cells (NR) at 4 and 11 dpi.The size of each bubble represents the number of functional protein domains allocated to each group.