Skip to main content

Comprehensive deciphering the alternative splicing patterns involved in leaf morphogenesis of Liriodendron chinense

Abstract

Alternative splicing (AS), a pivotal post-transcriptional regulatory mechanism, profoundly amplifies diversity and complexity of transcriptome and proteome. Liriodendron chinense (Hemsl.) Sarg., an excellent ornamental tree species renowned for its distinctive leaf shape, which resembles the mandarin jacket. Despite the documented potential genes related to leaf development of L. chinense, the underlying post-transcriptional regulatory mechanisms remain veiled. Here, we conducted a comprehensive analysis of the transcriptome to clarify the genome-wide landscape of the AS pattern and the spectrum of spliced isoforms during leaf developmental stages in L. chinense. Our investigation unveiled 50,259 AS events, involving 10,685 genes (32.9%), with intron retention as the most prevalent events. Notably, the initial stage of leaf development witnessed the detection of 804 differentially AS events affiliated with 548 genes. Although both differentially alternative splicing genes (DASGs) and differentially expressed genes (DEGs) were enriched into morphogenetic related pathways during the transition from fishhook (P2) to lobed (P7) leaves, there was only a modest degree of overlap between DASGs and DEGs. Furthermore, we conducted a comprehensively AS analysis on homologous genes involved in leaf morphogenesis, and most of which are subject to post-transcriptional regulation of AS. Among them, the AINTEGUMENTA-LIKE transcript factor LcAIL5 was characterization in detailed, which experiences skipping exon (SE), and two transcripts displayed disparate expression patterns across multiple stages. Overall, these findings yield a comprehensive understanding of leaf development regulation via AS, offering a novel perspective for further deciphering the mechanism of plant leaf morphogenesis.

Peer Review reports

Introduction

In eukaryotes, precursor RNA (pre-mRNA) splicing, which removed the introns from the nascent mRNAs and ligated the exons in the right order, is an essential step in translatable mRNA synthesis and gene expression [1, 2]. Besides constitutive splicing in per-mRNAs, spliceosome flexibly selected splice site to generate structurally variable transcripts from individual pre-mRNAs, known as alternative splicing (AS) [3]. As an important post-transcriptional regulation mechanism, AS can greatly enhance genome coding capacity, transcriptome plasticity and proteome diversity [4]. Genome-wide analysis of AS have been shown that extensively influence diverse biological processes in plants. These include responses to various environmental stresses and physiological signals, along with the dynamic regulation of plant development such as fruit development and early plant morphogenesis [5,6,7,8,9,10]. However, some alternative transcripts might not necessarily perform biological functions. AS often impacts mRNA stability by introducing premature termination codons (PTCs). Transcripts harboring PTCs are targeted by the nonsense-mediated decay (NMD) pathway and rapidly degraded to maintain mRNA stability [11, 12].

The occurrence of alternative splicing is contingent upon the spliceosome, a complex comprising five small nuclear ribonucleoprotein (snRNPs) particles (U1, U2, U4, U5, and U6) and an array of non-snRNP proteins, including the U2 auxiliary factor (U2AF) [11, 13]. These RNA-protein complexes effectively recognize the 5’ and 3’ splicing sites and define exons and introns in a precise and ordered assembly [14]. Furthermore, trans-acting splice factors including serine/arginine-rich (SR) proteins and heterogeneous nuclear ribonucleoproteins (hnRNPs) also regulate the occurrence of AS through binding to cis-regulatory elements within exons (ESRs) and introns (ISRs), respectively [15, 16]. Categorically based on splice site selection, five distinct alternatively spliced subtypes have been identified in plants: retained introns (RI), skipped exons (SE), alternative 5’ splice sites (A5SS), alternative 3’ splice sites (A3SS), and mutually exclusive exons (MEX) [17]. Nearly 95% of human genes experience alternative splicing, with SE standing out as the primarily form [18], whereas in plants, more than 60-85% of genes are regulated by AS, with RI being the most prevalent form [19,20,21]. Moreover, the distribution of various AS types may exhibit significant variability between dicots and monocots [22]. While the fundamental splicing mechanism remains comparatively conserved, the intricacies of intron processing and the principles of alternative splicing appear to have diverged among plant species. Consequently, a more comprehensive investigation is imperative in the domain of plants to unveil the multifaceted and diversified mechanisms governing alternative splicing.

Leaf morphology can be diverse to facilitate photosynthesis and gas exchange in plants [23]. Plant leaf development is a complex biological process in which plant cells coordinate cell growth, cell division and cell differentiation at the three axis of proximal-distal, medial-lateral and abaxial-adaxial [24]. Previous studies have shown that the regulation of transcription factors is widely involved in plant leaf development. For instance, the Class I KNOTTED1‐like homeobox proteins (KNOXI) and the ASYMMETRIC LEAVES1(AS1) transcription factors antagonistically regulate the initiation of leaf primordium [23, 25]. Another noteworthy factor, CUP-SHAPED COTYLEDON2 (CUC2) and PIN‐FORMED1(PIN1) feedback regulate the auxin maximum value at the serrated tip of the leaf margin to promote the intricate sculpting of leaf margin morphology [26]. Moreover, AINTEGUMENTA (ANT) and GA20ox1 regulated leaf size by influencing the contents of auxin and gibberellin (GA), thus reflecting the interconnectedness of growth mechanisms [27, 28]. It is apparent that most studies delved into transcriptional regulation, with limited researches focusing on the impact of post-transcriptional mechanisms on leaf morphogenesis and development. Yet, AS is being progressively considered an important machinery for modulating plant development and response to various environment stresses. Examination of global AS patterns in wild strawberry (Fragaria vesca L.) throughout different fruit development stages reveals that shifts in AS modes could potentially contribute to rapid changes in gene expression during fruit set [7]. GS3 is regulated by AS to fine-tune grain size in rice (Oryza sativa L.) [29]. Furthermore, large-scale studies in plants have concentrated on dissecting AS patterns and dynamic expression changes within their isoforms across different tissues or developmental phases [30, 31]. In general, the expression level of the constitutive transcripts is significantly higher than variable isoforms, however, certain conditions can elicit shifts in this expression pattern, leading to different dominant isoforms [32]. These events have been defined as differential transcript usage (DTU), a phenomenon that holds promise for further offering insights into the occurrence of AS or for pinpointing potential candidates for functional analysis [33].

As an excellent ornamental tree, Liriodendron chinense (Hemsl.) Sarg. is widely used in garden cultivation thanks to its unique leaf shape, gorgeous flowers and straight trunk. The distinctive ‘mandarin jacket’ leaf shape of Liriodendron is widely recognized, rendering it a valuable subject for unraveling mechanisms of leaf morphogenesis. Insights into the dynamic leaf development of L. chinense at the transcriptional level have suggested the potential involvement of various transcription factors and genes related to phytohormone metabolism in the process [34]. On this basis, preliminary functional analyses have been conducted on core genes such as LcKNOX, LcCUC2-like [35, 36]. However, previous investigations have predominantly focused on the transcriptional level, leaving a gap in the comprehensive exploration of its post-transcriptional regulatory mechanism. While Tu et al have provided an overview of genome-wide alternative splicing events in L. chinense [37], no attention has been paid to the specific biological process of leaf morphogenesis. It is noteworthy that the wealth of genome resources and transcriptome data, derived from multiple samples, provide the possibility for a profound and tailored analysis of alternative splicing at the post-transcription level [34, 37, 38]. To uncover genome-wide AS events during leaf morphogenesis of L. chinense, we embarked on multiple available transcriptomes were analyzed in this study [34, 37]. Firstly, we conducted Iso-seq to produce an enhanced genome annotation file. Subsequently, we reanalyzed the RNA-seq data using this improved annotation to retrieve genome-wide AS events. Then, differential splicing events and genes between different developmental stages of L. chinense leaves were identified using RNA-seq data. Finally, we meticulously annotated and characterized the transcripts involved in leaf morphogenesis. Taken together, through secondary in-depth analysis of the available transcriptome data, we revealed the alternative splicing landscape throughout the leaf morphogenesis of L. chinense and characterized key candidate genes regulated by AS. Our study provides a new perspective for further exploring the mechanism of plant leaf morphogenesis.

Materials and methods

Data sources and plant materials

Published PacBio full-length transcripts as well as Illumina data of L. Chinense were analyzed to unveil alternative splicing events. Total RNA from seven L. chinense tissues (shoot apices, leaves, bracts, sepals, petals, stamens and pistils, 21 samples) were mixed for SMRT bell libraries construction, while all samples also were used respectively for Illumina sequencing (PRJNA559687) [37]. Besides, Illumina RNA-seq data from different young leaves stages include fishhook shaped stage (P2) and deeper lobe stage (P7) were used for differential analysis to find alternative splicing genes regulating lobed leaf development (SRR8101043, SRR8101042, SRR8101041 and SRR8101040) [34]. In addition, plant materials used in the experiment were collected from the forest farm affiliated to Nanjing Forestry University, Jiangsu Province, China (119°13′20″E, 32°7′8″N). Various tissues were removed from a mature L. chinense tree that was originally from Lushan, Jiangxi Province (116°0′E, 29°32′N) and young leaves in distinct developmental stages were removed from an adulted tree originated from Wuyi Mountain, Fujian Province (117°0′E, 27°N).

Reads mapping and transcript assembly

An optimized pipeline was constructed for integrating transcriptome data resources (Fig. 1). Firstly, Illumina clean reads were mapped to the L. Chinense reference genome [38] using Tophat2 (v.2.1.1, http://ccb.jhu.edu/software/tophat/index.shtml) [39]. Based on alignment, the unique mapped reads were assembled by StringTie (v.2.2.1, https://ccb.jhu.edu/software/stringtie/) [40] with the default parameters. Then, the transcripts obtained from all samples were merged with StringTie --merge. Meanwhile, Minimap2 (v.2.24, https://github.com/lh3/minimap2) [41] was used to map high-quality full-length isoforms to the reference genome [42]. Mapped isoforms were further collapsed into non-redundant transcript groups using cDNA_Cupcake (v.29.0.0, https://github.com/Magdoll/cDNA_Cupcake).

Fig. 1
figure 1

The workflow for acquiring the improved annotation file and AS analysis

To obtain an improved annotation file, transcripts from both short-read RNA-seq data and long-read Iso-seq data were integrated to consolidate transcript information and remove redundancy using cuffmerge in the Cufflinks suit (v.2.2.1, https://github.com/cole-trapnell-lab/cufflinks). Transcripts with class_code “=” (exact match of intron chain) or “j” (multi-exon with at least one junction match) were considered as isoforms of known genes [7]. Based on the above, the improved transcriptomes were used as reference for downstream analysis.

Identification of leaf development stage differentially expressed genes (DEGs)

The fragments per kilobase per million mapped fragments (FPKM) values were calculated to estimate gene expression levels. The expression matrix of different stages was generated by cuffdiff in the Cufflinks suit (v.2.2.1, https://github.com/cole-trapnell-lab/cufflinks). DEseq2 (v.1.34.0, https://bioconductor.org/packages/DESeq2/) was used to calculate differential expression levels between two stages for leaf development [43]. Genes with p-value < 0.05 and |log2fold change| > 2 were further considered differentially expressed.

Identification of AS events and discovery of differential alternative splicing genes (DASGs) between leaf development stages

To classify the AS events, SUPPA (v.2.3, https://github.com/comprna/SUPPA) [44] and AStalavista (v.4.0) [45] were employed using the annotation file assembled from Iso-seq and Illumina RNA-seq data. Five main types of AS events were extracted and counted, namely RI (structure: 0,1^2-), SE (structure: '0,1-2^), A3 (structure: 1-,2-), A5 (structure: 1^,2^), MXE (structure: 1-2^,3-4^). AS location was visualized based on the assembled transcripts using R package RIdeogram [46]. Base-composition bias was visualized using R package ggplot2 (v.3.4.1, https://ggplot2.tidyverse.org/) [47] and ggseqlogo [48].

Moreover, aligned data from shoot apices, leaves, P2 and P4 were further used to identify differentially alternative splicing events between leaf development stages by SUPPA with default parameters. The ‘psiPerEvent’ was performed to extract the proportion spliced-in (PSI) value of local AS events in each sample. Then, the events that are differentially spliced between a pair of conditions were calculated with ‘diffSplice’. Differential events were filtered with P-value < 0.05. The events examined include retained introns (RI), retained partial exons in 3’ splice site (A3SS) or 5’ splice site (A5SS), skipped exons (SE) and two mutually to exclusive exons (MXE). Finally, the sashimi plot was plotted based on the built-in rmats2sashimiplot function of Rmats.

Gene Ontology (GO) enrichment analysis and visualization

The predicted coding region of transcripts were determined using TransDecoder (v.5.5.0, https://github.com/TransDecoder). Then, the protein from the longest ORF were functionally annotated using eggnog (v.5.0, http://eggnog5.embl.de/#/app-/home) [49]. Gene Ontology (GO) enrichment analysis was performed to explore the potential functions of DASGs and DEGs using the R package clusterProfiler (v.4.2.2, https://bioconductor.org/packages/release/bioc/html/cluster-Profiler.html) [50]. GO terms with P-value ≤ 0.01 were considered as significantly enriched categories. The R package ggplot2 (v.3.4.1, https://ggplot2.tidyverse.org/) was used for visualization.

RT-qPCR of AS transcripts

Total RNA was extracted from samples with the RNA-prep pure kit (Tiangen, Beijing, China) according to the instructions. Then, cDNA was synthesized from 500ng total RNA using Evo M-MLV RT Premix AG11706 (Accurate Biotechnology (Hunan) Co.,Ltd) in a 10 μL reaction volume. Using Oligo 7.0 software, primers for RT- qPCR were designed adhering to the instructions.

Phylogenetic analysis

Protein sequences of AIL5 isoforms from various plant species were obtained by conducting BLASTP search. These sequences were derived from NCBI database (https://www.ncbi.nlm.nih.gov/), FernBase (https://fernbase.org/) and multiple gymnosperms genomes [51,52,53,54]. Multiple sequences were aligned using Mafft (v 7.310) with default parameters [55]. Then, the unrooted phylogenetic tree was constructed using the maximum-likelihood (ML) method in IQ-TREE (v 1.6.12) [56]. The phylogenetic tree was visualized using iTOL (https://itol.embl.de/) online tool [57]. Gene structure was displayed using GSDS (v2.0, https://gsds.gao-lab.org/) [58].

Results

Overview of Liriodendron Chinense transcriptome assembly

To comprehensively investigate AS events in Liriodendron Chinense, we employed a dataset encompassing 25 RNA-seq samples. This dataset incorporated seven tissues (shoot apices, leaves, bracts, sepals, petals, stamens and pistils) as well as two leaf development stages (P2 and P7). Following alignment to the genome and the subsequent filtration of discordant and multiple alignments, an average of about 32.49 (82.31%) million unique mapped reads was utilized for subsequent transcriptional assembly (Supplementary Table S1). For each tissue, we independently constructed transcript assemblies according to a set of criteria. Then, the transcriptomes of these tissues were merged. Additionally, the annotation file was further improved by integrating the high-quality full-length isoforms and the reference genome. It should be noted that some transcripts were filtered out based on the “class code” and were not assembled within the reference annotation due to simple consideration for the main types of transcripts that match the reference genome. An updated annotation file was generated after filtering and removing redundancy as the reference to identify AS events (Supplementary File S1). From the annotation file, we identified a total of 90,376 transcripts associated with 32,445 genes. Notably, approximately 47% of these genes (15,306) featured multiple isoforms, collectively contributing 73,238 transcripts. Additionally, around 68.6% of transcripts (50,217) harbored more than one intron yet not exceeding ten introns.

Alternative splicing landscape in Liriodendron Chinense

To uncover genome-wide alternative splicing AS events, we analyzed the annotation file through both SUPPA and AStalavista respectively and focused on the five major types of alternative splicing events (Fig. 2a). A total of 50,259 AS events were identified simultaneously involving 10,685 genes (32.9 %), including retained intron (IR, 14,341 events, 28.5 %), skipped exon (ES, 11,588 events, 23.1 %), alternative 3’ splice site (A3SS, 13,569 events, 27 %), alternative 5’ splice site (A5SS, 9,659 events, 19.2 %) and mutually exclusive exons (MEX, 1,102 events, 2.2 %) (Fig. 2b, Supplementary File S2). Among them, IR were the most abundant splicing events, followed by A3SS, while MEX were relatively less frequent. Notably, the occurrence rate of IR events in L. chinense is considerably lower than what is observed in other plants [59].

Fig. 2
figure 2

Genome-wide alternative splicing (AS) landscape of L. Chinense. a The five representative forms of AS, namely alternative 3’ splice sites (A3SS), alternative 5’ splice sites (A5SS), retained introns (RI), skipping exons (SE) and mutually exclusive exons (MEX). b The number of each AS type identified using different software. c and d Analysis of splicing site bias, including the donor site (c) and acceptor site (d). e Distribution of AS events on chromosomes. Different colors represent different types of AS

The process of AS depends on the effectively recognition of the 5’ and 3’ splicing sites by splicing factors. Therefore, we analyzed the nucleotide frequency of transcripts using .gtf file to unravel the distribution of nucleotide at splicing sites. The results showed that obvious nucleotide bias was found at the donor/acceptor site. Remarkably, the core sequence elements ‘GT’ (donor site) (Fig. 2c) and ‘AG’ (acceptor site) (Fig. 2d) were identified as highly conserved, constituting 97.45% of all splice sites. This finding is consistent with the previously reported “GT—AG” rule of remove introns [12]. Besides, the U2 type splice site “GC—AG” constituted 1.44% of all splice sites, and the U12 type splice site “AT—AC” accounted for 0.17% (Supplementary Table S2). Despite the low distribution frequency, these splicing site types without canonical borders are also relatively conserved during pre-mRNA splicing across species [60, 61].

Additionally, we drew chromosome distribution maps to visually represent the occurrence of five distinct types of alternative splicing events across the genome, with the exception of those splicing events distributed within contigs (Fig. 2e). The distribution of these AS events exhibited an uneven pattern across all 19 chromosomes, while following a consistent trend across the five splicing types, that is, more distribution at both ends of the chromosome and few in the middle. Intriguingly, the distribution of AS events did not exhibit a strict positive correlation with chromosome length. For instance, Chr 2 harbored the highest proportion of AS events at 8.28%, whereas Chr 19 displayed the lowest percentage at 2.73%. Notably, even the shortest Chr 9 contributed 4.39% to the overall tally of AS events.

Characteristics of Alternative splicing events during leaf development

Subsequently, we focused on the AS events during the leaf development, encompassing the five major types, A3SS, A5SS, RI, SE and MXE, aiming to explore the role of alternative splicing in leaf morphogenesis (Fig. 3a). Through the analysis of Illumina RNA-seq data spanning four stages of leaf development, a total of 16,114, 14,973, 14,649 and 13,543 alternative splicing events were identified in SA, P2, P7 and LE samples, respectively (Fig. 3b, c). Notably, from the early stage of leaf morphogenesis to mature leaves, this temporal analysis revealed a gradual reduction in the number of splicing events, suggesting a propensity for alternative splicing to primarily occur during the active growth stages.

Fig. 3
figure 3

AS analysis at different stages of leaf development. a Four development stages of L. Chinense leaves, including shoot apices (SA), fishhook shaped stage (P2), deeper lobe stage (P7) and leaves (LE). b and c The number of each type of AS events (b) and genes (c) at different stages. d Venn diagrams of the same AS type at different development stages

Delving into the various splicing types, we observed that SE was the most abundant AS events, ranging from 28.7% to 30.4% across the four stages and involving the most AS genes (Fig. 3b, c). In contrast, RI events were less frequently detected, comprising only 13.8% to 16.8% of events and surpassing only MXE events in frequency. This distribution may be attributed to limitations in genome assembly quality and the short-read lengths of Illumina RNA-seq data. Furthermore, statistical analysis revealed that there were 1,620, 1,684, 671, 2,068, and 236 genes consistently undergoing A3SS, A5SS, RI, SE, and MXE events, respectively, across all four developmental stages (Fig. 3d), indicating the conservation of some AS events across different leaf development stages. However, the Venn diagram reveals that still numerous AS events are distinct at different developmental stages, indicating stage-specific regulation of AS (Fig. 3d).

Analysis of DASGs and DEGs during leaf development and corresponding pathways enrichment

To gain a deeper understanding into the dynamic shifts in splicing that unfold during leaf development, we performed a comparative analysis of the relative expression of alternative splicing events across four developmental stages. Subsequently, the differentially alternative splicing genes (DASGs) and events were identified, sorted and integrated. (Supplementary Table S3). Correspondingly, the same comparative groups were applied to discern significantly differentially expressed genes (DEGs) (Supplementary Table S4). Based on statistical findings, a total of 804 DAS events, sourced from 548 gene loci, were delineated through the calculation of differences in mean PSI values between P2 and P7 stages (Fig. 4a). Notably, during the transition from young lobed (P7) to mature leaves (LE), 3,734 DAS events were detected with the highest frequency, from 2,243 gene loci (Fig. 4b). Besides, a comparative analysis between shoot apices (SA) and leaves (LE) led to the identification of 1,111 DAS events from 744 gene loci (Fig. 4c). However, only 87 genes were identified as shared DASGs in all three aforementioned comparisons (Fig. 4d). By contrast, 347 genes exhibited DEG commonality across the array of the comparisons (Fig. 4e). The formation of leaf organs is an orderly and complex process, in which different isoforms perform their functions at different stages.

Fig. 4
figure 4

Differential alternative splicing genes (DASGs) and differentially expressed genes (DEGs) during leaf development, and corresponding pathways GO enrichment. a-c Overlaps between DASGs and DEGs during different developmental stages, including P2 vs P7 (a), P7 vs LE (b) and SA vs LE (c). d and e Venn diagrams of DASGs (d) and DEGs (e). f The top 20 most significantly enriched GO terms in biological process of DASGs and DEGs. The red font marks the pathways related to leaf morphogenesis

As an important post-transcriptional regulatory mechanism, alternative splicing affects gene transcription and expression dynamics. To explore the potential relationship between alternative splicing and gene expression, we engaged in a comparative analysis of DASGs and DEGs. Intriguingly, only 13 genes (2.37%) with DAS events were affiliated with differentially expressed during the development of tender leaves (Fig. 4a). There are relatively more overlaps, 454 genes (20.24%), between differential splicing and differential expression during the transformation process of mature leaves (Fig. 4b). In the three comparisons, despite not exhibiting the largest overlap between differentially spliced and differentially expressed genes between leaf buds and leaves, they demonstrate the highest proportion (27.55%) in differentially spliced genes (Fig. 4c). The overall comparison revealed a limited number of genes that establish an intersection between DASGs and DEGs.

To investigate the influence of development-induced alternative splicing events on biological processes, we performed GO functional enrichment analysis on both DASGs and DEGs. During the transition from fishhook (P2) to lobed (P7) leaves, enrichment analysis revealed that genes with DAS events are primarily involved in developmental and morphogenetic pathways, such as regulation of cell size, formation of plant organ boundary and anatomical boundary (Fig. 4f). Notably, these biological processes coincided with the development stages of leaves. Comparatively, differentially expressed genes are also mainly enriched in development and morphogenetic pathways, especially leaf and shoot morphogenesis (Fig. 4f). Nonetheless, it is noteworthy that the cohort of DEGs captured by these pathways exhibits significant divergence from the subset of DASGs captured, with only one intersection Lchi00862, a homeobox transcription factor (Supplementary Table S5). Moreover, the analysis of differential genes during the transition to mature leaves showed that both DASGs and DEGs exhibited predominant involvement in responding to light stimulation and photosynthesis pathway (Fig. 4f). During this stage, leaves generate carbohydrates and ATPs for plant growth via photosynthesis and metabolic activities. Previous studies paid more attention to the involvement of DEGs in plant growth and development, while neglecting the regulation of DASGs. Through statistical and enrichment analysis, it was found that the regulation of post-transcriptional level was also closely related to plant growth and development, and is relatively independent from the transcriptional regulation.

Features of DAS events in transcript isoforms during the transition from fishhook (P2) to lobed (P7) leaves

To specifically investigate the characteristics of DAS events during the initial stages of leaf development, we further focused on the comparison between P2 and P7. Proportion spliced-in (PSI) is utilized to quantify the relative abundances of the splicing events or transcript isoforms, while the difference of these relative abundances across conditions are denoted by differential splicing (dPSI). During the transition from fishhook (P2) to lobed (P7) leaves, dPSI tends to be lower than zero, encompassing 443 splicing events (55%, 443/804). This pattern indicated that developmental progression has resulted in an increased proportion of exclusion isoforms from differentially splicing genes (Fig. 5a).

Fig. 5
figure 5

The features of DAS events between fishhook (P2) and lobed (P7) leaves. a Volcano map showing the differential PSI (dPSI) levels. b Statistical analysis of dPSI value of different AS type. c-f Distribution of the dPSI value of different AS type, including AS33 (c), A5SS (d), RI (e) and SE (f). (g to j) Distribution of the sizes of variable region of different AS type, including AS33 (g), A5SS (h), RI (i) and SE (j)

Subsequently, we generated frequency distributions for inclusion or exclusion lengths and dPSI values, aiming to characterize different splicing types. The result showed that RI events predominantly exhibited dPSI values higher than zero, indicating a propensity for intron retention during development (Fig. 5b, e). Conversely, other splicing events, including SE, A3SS, A5SS and MEX, contain more events with dPSI less than zero, indicating that exons and 3’ and 5’ terminal sequences tend to be skipped during development (Fig. 5b-f). Furthermore, the length distribution of the alternative region is different for various splicing events. Skipped exons exhibited the broadest range, spanning from three to 8,367 nt, with a median length of 101 nt (Fig. 5g). In contrast, retained introns displayed a length distribution ranging from 25 to 5,817 nt, with a median length of 289 nt (Fig. 5h). The alternative splicing length of A3SS and A5SS were generally short, with a median length of 18 nt for A3SS and 38 nt for A5SS (Fig. 5i, j).

Alternative splicing events analysis of homologous gene related to leaf morphogensis

Based on previous research on leaf development, especially in the model plant Arabidopsis thaliana (L.) Heynh., we outlined a genetic mechanism network of leaf morphogenesis (Fig. 6a, b). Then, homologous genes in L. Chinense were obtained through protein sequence alignment, and some of these genes have been preliminarily validated for their functions in leaf development [35, 36]. Further survey revealed that some genes were regulated by alternative splicing post-transcriptional during leaf morphogenesis of L. Chinense, especially genes related to the leaf primordia initiation and the adaxial-abaxial polarity establishment (Supplementary Table S6). Two conservative mechanisms centered around KNOTTED1-like homeobox proteins (KNOXI) and PIN1 PIN-FORMED1 (PIN1) have been reported to be involved in leaf primordia, where the LcKNAT6 as homologous gene of KNOXI contains RI and A5 events, while A3 and SE events were identified in LcPIN1 gene. In addition, genes related to auxin homeostasis and transport, as well as plant growth regulation, such as AINTEGUMENTA-LIKE/PLETHORA (AIL/PLT), AUXIN RESPONSE FACTORS (ARF), GROWTH-REGULATING FACTOR (GRF) and basic/HELIX–LOOP–HELIX TRANSCRIPTION FACTOR (bHLH), are all involved in AS events. Regarding this phenomenon, it is speculated that during the initial stage of leaf development, external growth regulation and internal hormone homeostasis stimulate the expression of a substantial number of genes and the occurrence of AS events.

Fig. 6
figure 6

The genetic control of leaf morphogenesis based on previous reports. a Shoot apical meristems (SAM) maintenance and leaf primordia adaxial–abaxial polarity establishment. b Leaf growth along the proximal-distal and medio-lateral axis and leaf margin development. Red stars indicate homologous genes have undergone alternative splicing in L. Chinense, and details refer to Supplementary Table S6. c The sashimi plot of differentially spliced events of Lchi03252, which has been annotated as AIL5 from APETALA2 transcription factor family. d The sashimi plot of differentially spliced events of Lchi05740, which has been annotated as BOP1 acting on the adaxial side of leaf primordia

Subsequently, differential AS analysis was performed on these genes to investigate AS variation over development. The results showed that there were DAS events between the P2 and P7 stages for both AIL (Lchi13837, Lchi03252) and ROTUNDIFOLIA (ROT, Lchi23114), involving A5 and SE types (Supplementary Table S6). During the transition from P7 stage to mature leaves, DAS events were detected in BLOCK OF CELL PROLIFERATION (BOP, Lchi05740), BEL1-LIKE HOMEODOMAIN (BELL, Lchi23624) and PIN1 (Lchi16330). In addition, the SE event of ARF (Lchi18918) was detected differently between shoot apices and mature leaves. Notably, we conducted a detailed unravelling analysis of two genes, AIL5 (Lchi03252) and BOP1 (Lchi05740). AIL5 was significantly clustered into morphogenetic related pathways (Fig. 4f, Supplementary Table S5). This gene underwent an SE splicing event and had a significantly different between P2 and P7 stages (dPSI = -0.54, p-value = 0.005) (Fig. 6c). BOP1, a key gene that regulates the leaf’s proximal-distal patterning, produced an RI splicing event and significantly different between P7 and LE (dPSI = 0.43, p-value = 0.0025) (Fig. 6d). This spatiotemporal difference indicates that different isoforms could function at different development stages, these gene and its isoforms may be involved in leaf morphogenesis in ways that have not yet been explored.

Conservation of alternative splicing in LcAIL5 orthologs across plants

To further clarify the role of AIL5 in plant growth and development through AS, we conducted an analysis of the AS patterns of AIL5 orthologs across the plant kingdom. A total of protein sequences of 58 AIL5 genes and gene structure information of 70 isoforms from 48 plant species were obtained by conducting BLASTP search, involving 40 angiosperms, four gymnosperms and four ferns. Subsequently, we constructed a phylogenetic tree to analyze their phylogenetic relationship using the protein sequences encoded by these isoforms (Fig. 7a). Based on current database information, we discovered nine AIL5 genes with alternative splicing isoforms from eight species only. Since many sequences are retrieved from transcriptome data rather than a complete reference genome, the transcript information is not comprehensive and the actual number of isoforms may be far more than that. Nevertheless, an interesting phenomenon was found that seven out of the nine AIL5 genes had experienced SE events. Except for XP_021655431.1 gene, others all experienced the same forth exon skipped as LcAIL5. It is worth noting that the AIL5 gene (XP_020523523.1) from Amborella trichopoda, a basal angiosperm closely related to L. chinense, also underwent the same alternative splicing event. This further indicates that the AS patterns of AIL5 might exhibit evolutionarily conservation.

Fig. 7
figure 7

Phylogenetic and conservation analysis of AIL5. a Phylogenetic tree of LcAIL5 homologous proteins. The unrooted phylogenetic tree was constructed using the maximum-likelihood (ML) method in IQ-TREE. Isoform structure was constructed using GSDS 2.0. Blue structure, UTR region; yellow structure, CDS; red box, variable exon. b Conserved protein domain alignment of isoforms with SE splicing event. c Base bias analysis of alternative exon and splice site of AIL5 genes. d Relative expression of LcAIL5_AS1 and LcAIL5_AS2 in different developmental stages. The error bar shows Std. deviation

Through sequence alignment of proteins encoded by these genes and isoforms, it was found that the special exon that exhibited skipping encoded merely three amino acids, namely ‘VYL’. Notably, this motif is located in the first AP2 domain of AIL5 protein (Fig. 7b). Then, we performed a conservative analysis of the alternative exon and its 10 bp base near the two ends of the splice site (Fig. 7c). The result showed that these exon sequences were highly conserved except for the second and eighth bases. The acceptor and donor sites at the ends of the exon were ‘AG’ and ‘GT’ respectively, which conformed to the splicing rule [14]. Collectively, the same AS pattern exists in LcAIL5 and its homologous genes, indicating a potentially evolutionarily conserved role in governing biological processes.

Subsequently, we designed primers for RT-qPCR to detect the expression patterns of LcAIL5_AS1 (TCONS_00079661) and LcAIL5_AS2 (TCONS_00079662) produced by SE (Fig. 7d). The result showed that the relative expression level of LcAIL5_AS1 was higher in the early stage of leaf development (P2), and then gradually decreased. On the contrary, the relative expression level of LcAIL5_AS2 was higher in late stages (P7). This is consistent with the transcriptome, suggesting that the two transcripts may be involved in leaf morphogenesis at different stages.

Discussion

As an important post-transcriptional regulatory mechanism, alternative splicing (AS) is widely involved in various biological processes in plants [62, 63]. In addition, AS greatly enhances the coding capacity of the genome and enriches the diversity of the transcriptome and proteome [3]. Nowadays, high-throughput transcriptome sequencing and abundant genomic resources provide convenience for computational analysis of alternative splicing. AS has been widely explored in many plant species, including Arabidopsis thaliana [2, 64], rice (Oryza sativa) [20], maize (Zea mays) [16], cassava (Manihot esculenta) [6] and poplar (Populus trichocarpa) [8] and so on. Based on genome-wide integration analysis of RNA-seq data from multiple tissues and developmental stages, we detected more AS events in L. Chinense than previously reported [37]. In this study, approximately 32.9% (10,685/32,445) of genes were identified to undergo alternative splicing, involving a total of 50,259 AS events (Fig. 2b). This means that it might lead to encode more protein products that perform different functions [14]. A statistical classification of genome-wide AS events revealed that RI was the most abundant splicing event, aligning with the general characteristics of previous plant studies [65].

The process of alternative splicing not only significantly enhances the genome coding capacity, but also serves as an efficient regulatory approach towards achieving developmental plasticity in plants [59]. Transcriptome is widely used to explore various biological mechanisms of plants, including leaf morphogenesis [34, 66]. In the present study, we conducted an in-depth analysis of transcriptome data at different developmental stages of leaves in L. Chinense. In order to explore the effects of transcriptional and post-transcriptional regulation on leaf morphogenesis, we compared differentially alternative splicing genes (DASGs) and differentially expressed genes (DEGs) among different leaf development stages. The results suggested that less overlap was observed between DASGs and DEGs (Fig. 4a-c). The enrichment analysis results also showed that DASG and DEG were rarely enriched in common pathways (Fig. 4f). According to previous studies, the post-transcriptional level and the transcriptional level may be independent regulatory mechanisms [6, 42, 67]. In fact, AS generates different mRNA splicing isoforms encoding protein products with different functions without affecting gene expression. Differentially expressed genes encode the same protein product but change its expression abundance. Therefore, both DASGs and DEGs may also play important roles in the establishment of leaf morphology. In addition, the protein products caused by AS may serve as regulatory signals to affect the expression patterns of other genes [68, 69]. Collectively, these findings indicate that the interaction between AS and transcriptional modulations is intricate and need further investigation.

Leaves originate from the shoot apical meristem (SAM) and undergo cells proliferation and expansion in three-dimensional polarity, ultimately forming its final size and shape [70]. Previous studies have shown that genes functioning at leaf initiation stage determine leaf morphology [71], which consistent with the results of GO enrichment analyses. DASGs and DEGs were significantly enriched in the GO terms associated with morphogenesis during the transition from fishhook (P2) to lobed (P7) leaves at leaf initiation stage, including ‘formation of plant organ boundary’, ‘formation of anatomical boundary’, ‘shoot system morphogenesis’ and more (Fig. 4f). Subsequently, in-depth analysis of the enriched genes revealed that the four DASGs related to leaf morphogenesis were derived from different transcription factor families, including BLH (BEL1-LIKE HOMEODOMAIN), ARF, (AUXIN RESPONSE FACTORS), and AIL/PLT (AINTEGUMENTA-LIKE/PLETHORA), which have been reported to be involved in leaf development or morphogenesis [27, 72, 73]. Besides, a comprehensive investigation of genes related to leaf morphogenesis revealed that some of them underwent AS events during leaf development involving various AS types, and eight of them identified as DASGs with difference between developmental stages (Fig. 6a, b, Supplementary Table S6). In summary, the findings indicated that key DASGs involved in leaf morphogenesis, weaving a complex regulatory network.

Our comprehensive analysis revealed that AINTEGUMENTA-LIKE (AIL) transcript factor, LcAIL5, undergoes AS regulation during leaf initiation, with a SE event of the fourth exon (Fig. 6c). AIL proteins belong to the APETALA 2/ETHYLENE RESPONSE FACTOR (AP2/ERF) family and contain two AP2 conserved domains. This variable exon, encoding the three amino acids residues ‘VYL’, is characteristic of microexons and is located in the first conserved AP2 domain of the AIL5 protein (Fig. 7b). Members of the AIL clade are commonly expressed in all dividing tissues of plants and play central roles in various developmental processes, including embryogenesis, stem cell niche specification, meristem maintenance, organ positioning and growth [74]. Prasad et al. investigated the function loss mutants of Arabidopsis PLT3 (AIL3), PLT5 (AIL5) and PLT7 (AIL7) and discovered that these AIL proteins are key regulators of PIN1 activity in controlling phyllotaxis [75]. Additionally, the AIL5/PLT5 gene is expressed in the shoot tissues, including in vegetative and inflorescence shoot apical meristems (SAMs) [76]. These evidences indicate the important role of AIL5 gene in early morphogenesis and meristem maintenance of young organs such as leaf and flowers.

Despite the confirmed function of AIL5 in leaf morphogenesis in the model plant A. thaliana, the alternative splicing of this gene has not been explored. This study constructed a landscape of alternative splicing events in L. chinense through integrated transcriptome analysis during leaf shape development and discovered AIL5 as a key candidate gene undergoing a SE event of 'VYL'. This is the first report of alternative splicing of AIL5. Additionally, analysis of AIL5 homologous proteins across the plant kingdom revealed similar alternative splicing events in other species such as A. trichopoda and poplar (Populus trichocarpa Torr. & Gray) (Fig. 7a). The ‘VYL’ motif is not only present in the AIL5 protein, but also in other AP2 proteins. Krizek et al. used a yeast system to test the ability of a randomly mutagenized population of Arabidopsis AINTEGUMENTA (ANT) to bind to ANT target sequences in yeast. Site-directed mutation of either “Y” or “L” residue leads to complete impairment of transcriptional activation, whereas mutation of the “V” residue results in reduced ANT activation [77]. This indicated that VYL plays an important role in binding to the target DNA sequence. Another AP2 family gene, WRI1, a key transcription factor in the regulation of plant oil synthesis, has the same exon skipping pattern in Arabidopsis, oil palm (Elaeis guineensis Jacq.), oat (Avena sativa L.) and castor (Ricinus communis L.). The study demonstrated that site-directed mutagenesis within the VYL motif in AtWRI1 failed to restore the oil content in the wri1 mutant, indicating that VYL is important for its function [78]. In this study, the VYL domain exhibited differential expression at different developmental stages. The ‘VYL’ inclusion isoform was upregulated during the P2 stage, while the ‘VYL’ exclusion isoform displayed highly expression at the P7 stage. This expression pattern suggests a potential phased and alternating involvement of these two transcripts in leaf morphogenesis. However, their specific functions require further experimental investigation.

Conclusion

In summary, this study revealed the genome-wide alternative splicing (AS) landscape during leaf development in L. Chinense based on RNA-seq data and comprehensively deciphered the AS regulation mode of key genes involved in leaf morphogenesis. Based on the identified transcripts that are regulated by AS during leaf development, further functional characterization stands to enrich our understanding of the modulatory role of AS in leaf morphogenesis. Our findings also provide new perspectives and strategies for further deciphering the mechanism of plant leaf morphogenesis.

Availability of data and materials

The Supplementary File S1 and File S2 are freely accessible through the Figshare database (https://doi.org/10.6084/m9.figshare.25047311). The PacBio long reads and Illumina short reads were obtained from NCBI database with accession numbers of PRJNA559687, SRR8101043, SRR8101042, SRR8101041 and SRR8101040. All data generated or analyzed in this study are provided in this manuscript or supplementary information files.

References

  1. Chen W, Moore MJ. Spliceosomes. Curr Biol. 2015;25(5):R181–3.

    Article  CAS  PubMed  Google Scholar 

  2. Zhu D, Mao F, Tian Y, Lin X, Gu L, Gu H, Qu L-J, Wu Y, Wu Z. The features and regulation of co-transcriptional splicing in Arabidopsis. Mol Plant. 2020;13(2):278–94.

    Article  CAS  PubMed  Google Scholar 

  3. Kelemen O, Convertini P, Zhang Z, Wen Y, Shen M, Falaleeva M, Stamm S. Function of alternative splicing. Gene. 2013;514(1):1–30.

    Article  CAS  PubMed  Google Scholar 

  4. Verta J-P, Jacobs A. The role of alternative splicing in adaptation and evolution. Trends Ecol Evol. 2022;37(4):299–308.

    Article  CAS  PubMed  Google Scholar 

  5. Szakonyi D, Duque P. Alternative splicing as a regulator of early plant development. Front Plant Sci. 2018;9:1174.

    Article  PubMed  PubMed Central  Google Scholar 

  6. Li S, Yu X, Cheng Z, Zeng C, Li W, Zhang L, Peng M. Large-scale analysis of the cassava transcriptome reveals the impact of cold stress on alternative splicing. J Exp Bot. 2020;71(1):422–34.

    CAS  PubMed  Google Scholar 

  7. Li Y, Dai C, Hu C, Liu Z, Kang C. Global identification of alternative splicing via comparative analysis of SMRT- and Illumina-based RNA-seq in strawberry. Plant J. 2017;90(1):164–76.

    Article  CAS  PubMed  Google Scholar 

  8. Zhu F, Chen X, Song Y, Lam LPY, Tobimatsu Y, Gao B, Chen M, Cao F. SWATH-MS-based proteogenomic analysis reveals the involvement of alternative splicing in poplar upon lead stress. Genome Res. 2023;33(3):371–85.

    Article  PubMed  PubMed Central  Google Scholar 

  9. Wang R, Yin P, Ruixia Y, Liu X, Luo L, Xu J. Genome-wide profiling of alternative splicing genes in hybrid poplar (P.albaxP.glandulosa cv.84K) leaves. PLoS One. 2020;15(11):e0241914.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  10. Yan X, Bai D, Song H, Lin K, Pang E. Alternative splicing during fruit development among fleshy fruits. BMC Genomics. 2021;22:762.

    Article  PubMed  PubMed Central  Google Scholar 

  11. Filichkin S, Priest HD, Megraw M, Mockler TC. Alternative splicing in plants: directing traffic at the crossroads of adaptation and environmental stress. Curr Opin Plant Biol. 2015;24:125–35.

    Article  CAS  PubMed  Google Scholar 

  12. Reddy ASN, Marquez Y, Kalyna M, Barta A. Complexity of the alternative splicing landscape in plants. Plant Cell. 2013;25(10):3657–83.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  13. Wahl MC, Will CL, Lührmann R. The spliceosome: design principles of a dynamic RNP machine. Cell. 2009;136(4):701–18.

    Article  CAS  PubMed  Google Scholar 

  14. Lee Y, Rio DC. Mechanisms and regulation of alternative pre-mRNA splicing. Annu Rev Biochem. 2015;84:291–323.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  15. Laloum T, Martín G, Duque P. Alternative splicing control of abiotic stress responses. Trends Plant Sci. 2018;23(2):140–50.

    Article  CAS  PubMed  Google Scholar 

  16. Chen Q, Han Y, Liu H, Wang X, Sun J, Zhao B, Li W, Tian J, Liang Y, Yan J, et al. Genome-wide association analyses reveal the importance of alternative splicing in diversifying gene function and regulating phenotypic variation in maize. Plant Cell. 2018;30(7):1404–23.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  17. Barbazuk WB, Fu Y, McGinnis KM. Genome-wide analyses of alternative splicing in plants: opportunities and challenges. Genome Res. 2008;18(9):1381–92.

    Article  CAS  PubMed  Google Scholar 

  18. Wang ET, Sandberg R, Luo S, Khrebtukova I, Zhang L, Mayr CFS, Kingsmore GPS, Burge CB. Alternative isoform regulation in human tissue transcriptomes. Nature. 2008;456:470–6.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  19. Zhu F, Chen M, Ye N, Shi L, Ma K, Yang J, Cao Y, Zhang Y, Yoshida T, Fernie AR, et al. Proteogenomic analysis reveals alternative splicing and translation as part of the abscisic acid response in Arabidopsis seedlings. Plant J. 2017;91(3):518–33.

    Article  CAS  PubMed  Google Scholar 

  20. Chen M, Zhu F, Gao B, Ma K, Zhang Y, Fernie RA, Chen X, Dai L, Ye N, et al. Full-length transcript-based proteogenomics of rice improves its genome and proteome annotation. Plant Physiol. 2020;182(3):1510–26.

    Article  CAS  PubMed  Google Scholar 

  21. Hu H, Yang W, Zheng Z, Niu Z, Yang Y, Wan D, Liu J, Ma T. Analysis of alternative splicing and alternative polyadenylation in Populus alba var. pyramidalis by single-molecular long-read sequencing. Front Genet. 2020;11:48.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  22. Shen Y, Zhou Z, Wang Z, Li W, Fang C, Wu M, Ma Y, Liu T, Kong L, Peng D, et al. Global dissection of alternative splicing in paleopolyploid soybean. Plant Cell. 2014;26(3):996–1008.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  23. Wang H, Kong F, Zhou C. From genes to networks: the genetic control of leaf development. J Integr Plant Biol. 2021;63(7):1181–96.

    Article  CAS  PubMed  Google Scholar 

  24. Richardson A. Plant development: coordinating across space and time. Curr Biol. 2020;30(24):R1492–4.

    Article  CAS  PubMed  Google Scholar 

  25. Lodha M, Marco CF, Timmermans MCP. The ASYMMETRIC LEAVES complex maintains repression of KNOX homeobox genes via direct recruitment of Polycomb-repressive complex2. Genes Dev. 2013;27(6):596–601.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  26. Bilsborough GD, Runions A, Barkoulas M, Jenkins HW, Hasson A, Galinha C, Laufs P, Hay A, Prusinkiewicz P, Tsiantis M. Model for the regulation of Arabidopsis thaliana leaf margin development. Proc Natl Acad Sci U S A. 2011;108(8):3424–9.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  27. Mizukami Y, Fischer RL. Plant organ size control: AINTEGUMENTA regulates growth and cell numbers during organogenesis. Proc Natl Acad Sci U S A. 2000;97(2):942–7.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  28. Mitchum MG, Yamaguchi S, Hanada A, Kuwahara A, Yoshioka Y, Kato T, Tabata S, Kamiya Y, Sun T. Distinct and overlapping roles of two gibberellin 3-oxidases in Arabidopsis development. Plant J. 2006;45(5):804–18.

    Article  CAS  PubMed  Google Scholar 

  29. Liu L, Zhou Y, Mao F, Gu Y, Tang Z, Xin Y, Liu F, Tang T, Gao H, Zhao X. Fine-tuning of the grain size by alternative splicing of GS3 in rice. Rice. 2022;15:4.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  30. Li M, Wang R, Song B, Cao B, Sun M, Liu Y, Chen X, Zhang S, Li J, Zhao K, et al. Alternative splicing analysis provides insight into the domestication and improvement of pear (Pyrus pyrifolia) fruit. Sci Hortic-amsterdam. 2022;297:110976.

    Article  CAS  Google Scholar 

  31. Klepikova AV, Kasianov AS, Gerasimov ES, Logacheva MD, Penin AA. A high resolution map of the Arabidopsis thaliana developmental transcriptome based on RNA-seq profiling. Plant J. 2016;88(6):1058–70.

    Article  CAS  PubMed  Google Scholar 

  32. Guo W, Calixto CPG, Brown JWS, Zhang R. TSIS: an R package to infer alternative splicing isoform switches for time-series data. Bioinformatics. 2017;33(20):3308–10.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  33. Vaneechoutte D, Estrada AR, Lin Y, Loraine AE, Vandepoele K. Genome-wide characterization of differential transcript usage in Arabidopsis thaliana. Plant J. 2017;92(6):1218–31.

    Article  CAS  PubMed  Google Scholar 

  34. Ma J, Wei L, Li J, Li H. The analysis of genes and phytohormone metabolic pathways associated with leaf shape development in Liriodendron chinense via de novo transcriptome sequencing. Genes. 2018;9(12):577.

    Article  PubMed  PubMed Central  Google Scholar 

  35. Ma J, Mei G, Liu H, Li H. Overexpression of a novel LcKNOX transcription factor from Liriodendron chinense induces lobed leaves in Arabidopsis thaliana. Forests. 2019;11(1):33.

    Article  CAS  Google Scholar 

  36. Wen S, Li J, Hao Z, Wei L, Ma J, Zong Y, Li H. Overexpression of the LcCUC2-like gene in Arabidopsis thaliana alters the cotyledon morphology and increases rosette leaf number. PeerJ. 2022;10:e12615.

    Article  PubMed  PubMed Central  Google Scholar 

  37. Tu Z, Shen Y, Wen S, Zong Y, Li H. Alternative splicing enhances the transcriptome complexity of Liriodendron chinense. Front Plant Sci. 2020;11:578100.

    Article  PubMed  PubMed Central  Google Scholar 

  38. Chen J, Hao Z, Guang X, Zhao C, Wang P, Xue L, Zhu Q, Yang L, Sheng Y, Zhou Y, et al. Liriodendron genome sheds light on angiosperm phylogeny and species-pair differentiation. Nat Plants. 2018;5(1):18–25.

    Article  PubMed  PubMed Central  Google Scholar 

  39. Kim D, Pertea G, Trapnell C, Pimentel H, Kelley R, Salzberg SL. TopHat2: accurate alignment of transcriptomes in the presence of insertions, deletions and gene fusions. Genome Biol. 2013;14:R36.

    Article  PubMed  PubMed Central  Google Scholar 

  40. Pertea M, Pertea GM, Antonescu CM, Chang T-C, Mendell JT, Salzberg SL. StringTie enables improved reconstruction of a transcriptome from RNA-seq reads. Nat Biotechnol. 2015;33(3):290–5.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  41. Li H, Birol I. Minimap2: pairwise alignment for nucleotide sequences. Bioinformatics. 2018;34(18):3094–100.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  42. Jin Z, Lv X, Sun Y, Fan Z, Xiang G, Yao Y. Comprehensive discovery of salt-responsive alternative splicing events based on Iso-Seq and RNA-seq in grapevine roots. Environ Exp Bot. 2021;192:104645.

    Article  CAS  Google Scholar 

  43. Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014;15:550.

    Article  PubMed  PubMed Central  Google Scholar 

  44. Trincado JL, Entizne JC, Hysenaj G, Singh B, Skalic M, Elliott DJ, Eyras E. SUPPA2: fast, accurate, and uncertainty-aware differential splicing analysis across multiple conditions. Genome Biol. 2018;19:40.

    Article  PubMed  PubMed Central  Google Scholar 

  45. Foissac S, Sammeth M. ASTALAVISTA: dynamic and flexible analysis of alternative splicing events in custom gene datasets. Nucleic Acids Res. 2007;35(2):W297–9.

    Article  PubMed  PubMed Central  Google Scholar 

  46. Hao Z, Lv D, Ge Y, Shi J, Weijers D, Yu G, Chen J. RIdeogram: drawing SVG graphics to visualize and map genome-wide data on the idiograms. PeerJ Comput Sci. 2020;6:e251.

    Article  PubMed  PubMed Central  Google Scholar 

  47. Wickham H. ggplot2: elegant graphics for data analysis. New York: Springer-Verlag; 2016.

    Book  Google Scholar 

  48. Wagih O. ggseqlogo: a versatile R package for drawing sequence logos. Bioinformatics. 2017;33(22):3645–7.

    Article  CAS  PubMed  Google Scholar 

  49. Huerta-Cepas J, Szklarczyk D, Heller D, Hernández-Plaza A, Forslund SK, Cook H, Mende DR, Letunic I, Rattei T, Jensen LJ, et al. eggNOG 5.0: a hierarchical, functionally and phylogenetically annotated orthology resource based on 5090 organisms and 2502 viruses. Nucleic Acids Res. 2019;47(D1):D309–14.

    Article  CAS  PubMed  Google Scholar 

  50. Wu T, Hu E, Xu S, Chen M, Guo P, Dai Z, Feng T, Zhou L, Tang W, Zhan L, et al. clusterProfiler 4.0: A universal enrichment tool for interpreting omics data. Innovation. 2021;2(3):100141.

    CAS  PubMed  PubMed Central  Google Scholar 

  51. Liu H, Wang X, Wang G, Cui P, Wu S, Ai C, Hu N, Li A, He B, Shao X, et al. The nearly complete genome of Ginkgo biloba illuminates gymnosperm evolution. Nat Plants. 2021;7(6):748–56.

    Article  CAS  PubMed  Google Scholar 

  52. Fu F, Song C, Wen C, Yang L, Guo Y, Yang X, Shu Z, Li X, Feng Y, Liu B, et al. The Metasequoia genome and evolutionary relationships among redwoods. Plant Commun. 2023;4(6):100643.

  53. Nystedt B, Street NR, Wetterbom A, Zuccolo A, Lin Y-C, Scofield DG, Vezzi F, Delhomme N, Giacomello S, Alexeyenko A, et al. The Norway spruce genome sequence and conifer genome evolution. Nature. 2013;497(7451):579–84.

    Article  CAS  PubMed  Google Scholar 

  54. Wegrzyn JL, Lee JM, Tearse BR, Neale DB. TreeGenes: a forest tree genome database. Int J Plant Genomics. 2008;2008:412875.

    Article  PubMed  PubMed Central  Google Scholar 

  55. Katoh K, Standley DM. MAFFT multiple sequence alignment software version 7: improvements in performance and usability. Mol Biol Evol. 2013;30(4):772–80.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  56. Nguyen L-T, Schmidt HA, von Haeseler A, Minh BQ. IQ-TREE: a fast and effective stochastic algorithm for estimating maximum-likelihood phylogenies. Mol Biol Evol. 2015;32(1):268–74.

    Article  CAS  PubMed  Google Scholar 

  57. Letunic I, Bork P. Interactive Tree Of Life (iTOL): an online tool for phylogenetic tree display and annotation. Bioinformatics. 2007;23(1):127–8.

    Article  CAS  PubMed  Google Scholar 

  58. Hu B, Jin J, Guo A-Y, Zhang H, Luo J, Gao G. GSDS 2.0: an upgraded gene feature visualization server. Bioinformatics. 2015;31(8):1296–7.

    Article  PubMed  Google Scholar 

  59. Muhammad S, Xu X, Zhou W, Wu L. Alternative splicing: An efficient regulatory approach towards plant developmental plasticity. Wiley Interdiscip Rev RNA. 2023;14(3):e1758.

    Article  CAS  PubMed  Google Scholar 

  60. Sheth N, Roca X, Hastings ML, Roeder T, Krainer AR, Sachidanandam R. Comprehensive splice-site analysis using comparative genomics. Nucleic Acids Res. 2006;34(14):3955–67.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  61. Moyer DC, Larue GE, Hershberger CE, Roy SW, Padgett RA. Comprehensive database and evolutionary dynamics of U12-type introns. Nucleic Acids Res. 2020;48(13):7066–78.

    CAS  PubMed  PubMed Central  Google Scholar 

  62. Chen M, Zhang K, Zhang M, Das D, Fang Y, Dai L, Zhang J, Zhu F. Alternative splicing and its regulatory role in woody plants. Tree Physiol. 2020;40(11):1475–86.

    Article  CAS  PubMed  Google Scholar 

  63. Wang J, Eulgem T. The Arabidopsis RRM domain proteins EDM3 and IBM2 coordinate the floral transition and basal immune responses. Plant J. 2023;116(1):128–43.

    Article  CAS  PubMed  Google Scholar 

  64. Martín G, Márquez Y, Mantica F, Duque P, Irimia M. Alternative splicing landscapes in Arabidopsis thaliana across tissues and stress conditions highlight major functional differences with animals. Genome Biol. 2021;22:35.

    Article  PubMed  PubMed Central  Google Scholar 

  65. Romanowski A, Schlaen RG, Perez-Santangelo S, Mancini E, Yanovsky MJ. Global transcriptome analysis reveals circadian control of splicing events in Arabidopsis thaliana. Plant J. 2020;103(2):889–902.

    Article  CAS  PubMed  Google Scholar 

  66. Yao X, Qi Y, Chen H, Zhang B, Chen Z, Lu L. Study of Camellia sinensis diploid and triploid leaf development mechanism based on transcriptome and leaf characteristics. PLoS One. 2023;18(2):e0275652.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  67. Liu Z, Qin J, Tian X, Xu S, Wang Y, Li H, Wang X, Peng H, Yao Y, Hu Z, et al. Global profiling of alternative splicing landscape responsive to drought, heat and their combination in wheat (Triticum aestivum L.). Plant Biotechnol J. 2018;16(3):714–26.

    Article  PubMed  Google Scholar 

  68. Ye L, Wu Y, Zhang J, Zhang J, Zhou H, Zeng R, Zheng W, Qiu M, Zhou J, Xie Z, et al. A bZIP transcription factor (CiFD) regulates drought- and low-temperature-induced flowering by alternative splicing in citrus. J Integr Plant Biol. 2023;65(3):674–91.

    Article  CAS  PubMed  Google Scholar 

  69. Wang Z, Hong Y, Yao J, Huang H, Qian B, Liu X, Chen Y, Pang J, Zhan X, Zhu J, et al. Modulation of plant development and chilling stress responses by alternative splicing events under control of the spliceosome protein SmEb in Arabidopsis. Plant Cell Environ. 2022;45(9):2762–79.

    Article  CAS  PubMed  Google Scholar 

  70. Bhatia N, Runions A, Tsiantis M. Leaf shape diversity: from genetic modules to computational models. Annu Rev Plant Biol. 2021;72:325–56.

    Article  CAS  PubMed  Google Scholar 

  71. Wilson-Sánchez D, Bhatia N, Runions A, Tsiantis M. From genes to shape in leaf development and evolution. Curr Biol. 2022;32(21):R1215–22.

    Article  PubMed  Google Scholar 

  72. Wen L, Kong Y, Wang H, Xu Y, Lu Z, Zhang J, Wang M, Wang X, Han L, Zhou C. Interaction between the MtDELLA-MtGAF1 complex and MtARF3 mediates transcriptional control of MtGA3ox1 to elaborate leaf margin formation in Medicago truncatula. Plant Cell Physiol. 2021;62(2):321–33.

    Article  CAS  PubMed  Google Scholar 

  73. Kumar R, Kushalappa K, Godt D, Pidkowich MS, Pastorelli S, Hepworth SR, Haughn GW. The Arabidopsis BEL1-LIKE HOMEODOMAIN proteins SAW1 and SAW2 act redundantly to regulate KNOX expression spatially in leaf margins. Plant Cell. 2007;19(9):2719–35.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  74. Horstman A, Willemsen V, Boutilier K, Heidstra R. AINTEGUMENTA-LIKE proteins: hubs in a plethora of networks. Trends Plant Sci. 2014;19(3):146–57.

    Article  CAS  PubMed  Google Scholar 

  75. Prasad K, Grigg SP, Barkoulas M, Yadav RK, Sanchez-Perez GF, Pinon V, Blilou I, Hofhuis H, Dhonukshe P, Galinha C, et al. Arabidopsis PLETHORA transcription factors control phyllotaxis. Curr Biol. 2011;21(13):1123–8.

    Article  CAS  PubMed  Google Scholar 

  76. Nole-Wilson S, Tranby TL, Krizek BA. AINTEGUMENTA-like (AIL) genes are expressed in young tissues and may specify meristematic or division-competent states. Plant Mol Biol. 2005;57(5):613–28.

    Article  CAS  PubMed  Google Scholar 

  77. Krizek BA. AINTEGUMENTA utilizes a mode of DNA recognition distinct from that used by proteins containing a single AP2 domain. Nucleic Acids Res. 2003;31(7):1859–68.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  78. Ma W, Kong Q, Arondel V, Kilaru A, Bates PD, Thrower NA, Benning C, Ohlrogge JB. Wrinkled1, a ubiquitous regulator in oil accumulating tissues from Arabidopsis embryos to oil palm mesocarp. PLoS One. 2013;8(7):e68887.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

Download references

Funding

This work is supported by the National Natural Science Foundation of China (31770718, 32371910) and by the Priority Academic Program Development of Jiangsu Higher Education Institutions (PAPD).

Author information

Authors and Affiliations

Authors

Contributions

Y. Z. performed the experiments, analyzed the data and wrote the manuscript. F. Z., H. W. and H. X. re-annotated the genes and analyzed the data. J. W., Z. T., and L. Y. performed the experiments and drew the figures. H. L. designed the project and revised the manuscript. All authors contributed to the article and approved the submitted version.

Corresponding author

Correspondence to Huogen Li.

Ethics declarations

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Additional file 1:

 Supplementary Table S1. Mapping rates of all RNA-seq data.

Additional file 2:

 Supplementary Table S2. Types of splice junctions.

Additional file 3:

 Supplementary Table S3. Differentially alternative splicing genes (DASGs) and events during leaf development.

Additional file 4:

 Supplementary Table S4. Differentially expressed genes (DEGs) during leaf development.

Additional file 5:

 Supplementary Table S5. DASGs and DEGs enriched into morphogenetic related pathways.

Additional file 6:

 Supplementary Table S6. List of genes involved in leaf morphogenesis.

Additional file 7: Supplementary

Table S7. The information aboutAIL5 homologous genes from 58 plant species.

Additional file 8: Supplementary File S1. Updated annotation file in GFF3 format.

Additional file 9: Supplementary File S2. Genome wide alternative splicing events in L. Chinense.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Zong, Y., Zhang, F., Wu, H. et al. Comprehensive deciphering the alternative splicing patterns involved in leaf morphogenesis of Liriodendron chinense. BMC Plant Biol 24, 250 (2024). https://doi.org/10.1186/s12870-024-04915-x

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s12870-024-04915-x

Keywords