Signatures of mRNA Alternative Polyadenylation in Arabidopsis Leaf Development

Alternative polyadenylation (APA) of pre-mRNA is an important co-transcriptional mechanism that modulates gene expression, leading to transcriptomic and functional diversities. The role of APA in Arabidopsis leaf development, however, remains elusive. We applied a poly(A)-tag sequencing (PAT-seq) technique to characterize APA-mediated regulation events in cotyledon and in five stages of true leaf development. Over 60% APA was identified in genes expressed in leaves, consistent with the results in previous publications. However, a reduced APA level was detected in younger leaves, reaching 44% in the 18th true leaf. Importantly, we also found that >70% of the poly(A) site usages were altered in the second true leaf relative to the cotyledon. Compared with the cotyledon, more genes in the second true leaf tended to use the distal site of 3′UTR, but this was not found in pairwise comparison among other true leaves. In addition, a significant APA gene was found to be decreased in a pairwise comparison among true leaves, including differentially expressed genes. The APA genes identified herein were associated with specific biological processes, including metabolic and cellular processes and response to stimuli and hormones. These results provide a new insight into the regulation of Arabidopsis leaf development through APA.

Introduction

Leaves are the main plant organ responsible for photosynthesis through light energy harvesting and driving circulation through transpiration. Leaf morphogenesis in eudicots is divided into four stages. First, the founder cells are recruited from the peripheral zone of the shoot apical meristem (SAM) at the site of the incipient leaf primordium to undergo extensive cell division. Second, the growth occurs along the adaxial–abaxial, proximal–distal, and medio–lateral axes. Third, the leaf blade and petiole regions are established. Finally, across the entire blade, cell proliferation and expansion occur, which leads to both distal and lateral leaf expansion (Du et al., 2018). In this dynamic and precise process, extensive coordination and regulation of gene activities are crucial, which requires the whole program of gene expression circuitry, including mRNA processing control, e.g., pre-mRNA splicing and polyadenylation.

Indeed, the leaf development is controlled by a complex mechanism (Bar and Ori, 2014). During leaf initiation, the expression of Shoot Meristemless (STM) is required for maintaining the fate of the SAM cells in Arabidopsis (Long et al., 1996). Once a leaf primordium has been initiated, a peg-like structure grows and proliferates from the flank of the SAM along its adaxial–abaxial, proximal–distal, and medio–lateral axes to establish leaf polarity. In Arabidopsis thaliana, members of the class III homeodomain-leucine zipper (HD-ZIP III) family, including revoluta (REV), phabulosa, and phavoluta, specify the cell identity on the adaxial surface of the leaf (Manuela and Xu, 2020). In fact, HD-ZIP III genes are regulated by miR165 and miR166 (Zhu et al., 2011). In addition, the auxin response transcription factors ARF3 and ARF4 play a key role in adaxial/abaxial polarity (Sarojam et al., 2010). Then, growth along the proximal–distal axis determines the length of the leaf. The genes Blade-On-Petiole1 (BOP1) and BOP2 belong to the BTB family, which regulates the proximal–distal patterning of the leaf by inducing AS2 and repressing KNOXI (Ha et al., 2003; Jun et al., 2010). Finally, the growth along the medio-lateral axis determines the width of the leaf. Pressed flower/wuschel-like homeobox3 (PRS/WOX3) and WOX1, the two most important regulators, promote medio–lateral growth (Nakata et al., 2012). Interestingly, WOX1 and PRS genes are activated by the HD-ZIP III genes. The expression of the HD-ZIP III gene REV overlaps with the expression of PRS (Guan et al., 2019). Therefore, the steps in the establishment of leaf polarity may overlap with each other.

Following initiation and polarity establishment, the leaf primordium begins to expand until it has acquired its final shape and size. The molecular machinery underlying this progression balances cell proliferation with cell expansion. Auxin regulates both the AP2-domain transcription factor aintegumenta and CYCD3; 1 to influence both cell proliferation and cell expansion (Mizukami and Fischer, 2000). The deficiency of ARF2 promotes cell proliferation and leads to an increase in the leaf size (Okushima et al., 2005). In addition, gibberellin and brassinosteroids (BR) also promote the growth of the leaf blade because the deficient mutants of these genes show small leaves, whereas the overexpression of these genes produces larger than normal leaves (Mitchum et al., 2006; Oh et al., 2011; Zhiponova et al., 2013). As abovementioned, although it is well known that leaf development is modulated by transcription factors, hormones, miRNAs, and so on, the mechanism of leaf ontogeny remains to be elucidated. The overall process of leaf ontogeny proceeds along a fairly long timescale, but the leaf morphology is predictable. The continuous process of typical dicot leaf ontogeny and morphogenesis has been monitored using anatomical markers (Poethig, 1997). Furthermore, analyses of transcriptome dynamics during leaf development have identified multiple expression patterns of gradual increase, decrease, or a single expression peak with age. From the analyses of leaf development transcriptome, a quantitative definition of the dynamic leaf was established, termed as a digital differentiation index, to predict the leaf age (Efroni et al., 2008). However, the role of co-transcriptional regulation in leaf development has not been reported.

Polyadenylation is an important co-transcriptional processing step, which is essential for the maturation of eukaryotic precursor mRNAs (pre-mRNAs) before they are exported from the nucleus to the cytoplasm (Elkon et al., 2013). This process is a two-step reaction, cleavage of the nascent mRNAs, and addition of a poly(A) tail, which requires the recognition of the poly(A) signal located in the pre-mRNA sequence with a multi-protein complex, and regulates mRNA stability, nuclear export, and translation (Thomas et al., 2012). Expressed sequence tags and high-throughput studies have found that over 70% of genes in mammals, yeast, and plants harbor multiple cleavages and polyadenylation sites (poly(A) sites), and this phenomenon is termed as “alternative polyadenylation” (APA) (Shen et al., 2008; Wu et al., 2011; Derti et al., 2012). Poly(A) sites can be located at different regions of a gene, which may impact the length of the encoded protein or alter the presence of the regulatory elements in the 3′UTR to modulate mRNA stability. Thus, APA enriches the complexity and diversity of transcriptome and proteome and plays a key role in the regulation of gene expression (Xing and Li, 2011). Accumulating evidence has demonstrated that APA is involved in numerous biological processes and it is also subject to many regulations, including the epigenetic process (Lin et al., 2020). In humans, APA is related to disease occurrence, immunity regulation, cancer formation, and cell reprogramming (Mayr and Bartel, 2009; Carpenter et al., 2014). In plants, APA functions in flowering time regulation, seed dormancy, root development, and stress response (Zhang et al., 2008; Bruggeman et al., 2014; Zhang et al., 2015; Cyrek et al., 2016; Lin et al., 2017). The role of APA in leaf development has not yet been reported.

In this study, we applied a poly(A) tag sequencing (PAT-seq) approach to study dynamic APA during Arabidopsis leaf development. We found that the ratio of APA genes shows a downward trend in the progression of leaf development, suggesting that APA might be involved in this process. In addition, the number of significant APA genes and DE genes was found to decrease in the true leaf comparison. GO terms of significant APA genes were related to metabolic processes, cellular processes, and response to a stimulus. In addition, during true leaf development, there is a tendency that more genes use proximal sites in 3′UTR. Our results show that genome-wide APA mediates Arabidopsis leaf development and contributes to the molecular mechanism of the regulatory process of gene expression.

Materials and Methods

Plant Growth Conditions

Arabidopsis thaliana ecotype Columbia (Col-0, ABRC stock # CS60000) plants were grown on soil in a growth chamber condition under a light/dark cycle of 12/12 h at 22°C. Then cotyledon and different true leaves, including the second true leaf (2 TL), 6, 10, 14, and 18TL, were collected containing three biological replicates (each repeat contained ten leaves, each leaf from an individual plant) in 5-week-old plants, and the leaves were placed in liquid nitrogen and stored at −80°C.

RNA Extraction and PAT-Seq Library Construction

Total RNA extraction was performed using the TRIzol reagent (Invitrogen). Following isolation, the DNA in RNA samples was removed using DNase I (Takara) and was purified in an RNA clean column, and then 2 µg of the total RNA was used to construct PAT-seq libraries as described previously (Yu et al., 2019). Briefly, the total RNA was fragmented at 94°C for 4 min, followed by poly(A) enrichment with oligo (dT)₂₅ magnetic beads (New England Biolabs). Next, the fragments of poly(A) enriched RNA was reverse transcribed using barcoded oligo (dT) primers with the SMARTScribe enzyme (Clontech) for 2 h, and then 5’ adapter for template switching was added after 2 h. The cDNAs were purified using AMPURE XP beads (Beckman), and then were amplified with PCR for 18 cycles with Phire II (Thermo Fisher Scientific). The PCR product was run on a 2% agarose gel, and 300- to 500 bp fragments were purified to produce PAT-seq libraries. Qualified and quantified libraries were checked by Agilent Bioanalyzer 2100, Qubit 2.0, and qPCR.

High-Throughput Sequencing and Data Processing

The libraries were sequenced on the Illumina HiSeq 2500 platform at the facility located in the College of the Environment and Ecology, Xiamen University. The sequencing data were processed as described previously (Wu et al., 2011; Yu et al., 2019). Briefly, the low-quality reads, barcodes, and poly(T) stretch of raw reads were filtered out and trimmed. Then valid poly(A) reads were mapped to the Arabidopsis reference genome using Bowtie2 software. Next, a list of poly(A) clusters (PACs) was generated, where poly(A) sites (PAs) within 24 nucleotides were grouped into one poly(A) cluster (PAC). For the subsequent analyses, total reads of PACs among all samples with <10 were discarded to avoid uncertainty from low read counts.

Poly(A) Site Usage Analysis

Poly(A) site usage (PAU) represents the ratio of reads in one PAC relative to a gene total reads (Ha et al., 2018). The average of PAU among three biological replicates were used for calculating cumulative distribution function (CDF) and plotted using the mountain plot package in R. The details of the methods were described previously (Yu et al., 2019). Then genes passing through the following filtering criteria were considered as APA switching events: 1) total read count >10 for a gene; 2) absolute value (PAU of one age leaf–PAU of another age leaf)≥0.1; 3) at least one PAC belongs to significant differences PAC.

The 3′UTRs of TAIR9 genome annotation were extended by 120 nt. For those protein-coding genes without an annotated 3′UTR, their 3′UTRs lengths were defined as 338 nt, which is the average 3′UTR length 218 nt in Arabidopsis plus 120 nt extended lengths after the annotated stop codons, as defined previously (Wu et al., 2011). For 3′UTR APA analyses, the average weighted length of each 3′UTR of a gene was calculated as described previously (Fu et al., 2016). The 3′UTR length of the proximal poly(A) site was defined as a1 nt (the distance from the poly(A) site to the start of the 3′UTR), PAT number is b1; the 3′UTR length of distal poly(A) site was defined as a2 nt, PAT number is b2. Then, the average weighted length of 3′UTR of each gene was (a1*b1+a2*b2)/(b1+b2). A cut-off p value of 0.01 was used to justify significantly longer or shorter 3′UTR.

Genome-Wide DE-PACs, DE Genes, and GO Analysis

Differential expression PACs (DE-PAC) were identified using the DESeq2 package. The criterion of a statistically significant difference was an adjusted p value < 0.05. All PATs of a gene were summed for representing gene expression levels. Similarly, DE genes were analyzed using the DESeq2 package and an adjusted p value < 0.05 threshold was considered to mark significant differences. GO analyses of specific DE-PAC APA genes among different comparison groups were carried out using agriGO (v.2.0), singular enrichment analyses (SEA) was chosen and the Arabidopsis Information Resource 10th annotation acted as the background. False discovery rate (FDR) corrected p value < 0.05 were selected as statistically significant.

Results

Dynamic APA in Arabidopsis Leaf Development

We set out to understand the involvement of genome-wide APA during Arabidopsis leaf development. To this end, we employed PAT-seq to profile dynamic polyadenylation across cotyledon (abbreviated as Coty) and different true leaf (or TL) according to leaf development sequence, including the second true leaf (2 TL), the sixth true leaf (6 TL), Supplementary Figure S1 gives an overall view of the status of each leaf in a pool of wild-type plants at 5-week post-germination. The principal component analysis indicated that three independent repetitions showed good similarity (Supplementary Figure S2), indicating that the sequence data were reliable. Then, a summary of the raw reads, mapped poly(A) tag reads (PATs), and poly(A) site clusters (PACs) for each library are presented in Table 1, suggesting that sufficient data were obtained for further analyses. The distribution of PACs and PATs at the different genic regions was analyzed (Supplementary Figure S3). Altogether, 60% of PACs and over 80% of PATs fell within the annotated 3′UTRs. However, no significant difference was seen in terms of PAT and PAC distribution between different leaves. In total, 56,620 PACs were obtained, as shown in Supplementary Data 1. These PACs mapped to 19,612 genes (69% of total genes in Arabidopsis). Of these, the percentage of more than one PAC (defined as APA genes) is 63% (12,357/19,612) (Figure 1A), which is similar to previous research (Wu et al., 2011). Interestingly, we found that the ratio of APA genes shows a downward trend during different leaf development, in particular, it drops to 44% in the 18 TL (Figure 1B), suggesting that APA is related to different leaf development.

TABLE 1

TABLE 1. Summary of PAT mapping results.

FIGURE 1

FIGURE 1. Analysis of APA gene numbers and ratios. (A) Numbers of genes with different counts of poly(A) sites. PA site, Poly(A) site. (B) Ratio of APA and non-APA genes in each leaf stage. Cotyledon, coty; true leaf, TL; number before TL indicates the order of the leaves, for example, 2TL, the second true leaf.

To further illustrate the role of APA in leaf development, the relative abundance of each isoform within one gene was calculated and referred to as poly(A) site usage (PAU). PAU values were plotted using a cumulative distribution function (CDF) at the genic region level (Figure 2). Compared to cotyledon, 70% PAU of 2 TL was altered, because CDF curves merge 0.7 on the right y-axis (Figure 2A). In addition, the start site of CDF and the median (the point where the curves are folded) differed among 2 TL and cotyledon, indicating unique poly(A) site was used between these two samples. During true leaf development, the maximum distance of adjacent leaf comparison of CDF had a gradually increased tendency (Figure 2B) indicating that APA might play an important role in mediating leaf development. However, the CDF curves almost coincided except for slight differences at the start and end sites. These could be due to having a similar leaf shape (14 TL and 18 TL), or no overall difference in the cumulative fashion (2 TL and 10 TL). Nevertheless, there seem to have three different patterns among these comparisons (Figure 2B).

FIGURE 2

FIGURE 2. Curve of cumulative distribution function (CDF) based on poly(A) site usage (PAU). PAU values were calculated as the ratio of its expression in all isoforms of a gene. The x-axis shows the log PAU based on the average of three biological replicates. The y-axis shows folded empirical CDF. The curve of CDF was based on a mountain plot to examine PAU distribution. The Kolmogorov–Smirnov (K–S) test was used to detect the differences in both location and shape of the empirical CDF between different samples. The numbers in the figure insert represent the maximum distance in paired comparison between two CDF curves. **p value < 0.01. (A) Comparison of cotyledon and the second true leaf. (B) Comparison between true leaves. Leaf annotations are the same as those in Supplementary Figure S1.

Identification of PACs and Genes Displaying Differential Usages

To further delve into the role of APA in leaf development, differentially utilized PAC APA genes were characterized using the DESeq2 package. DE-PACs tables were as shown in Supplementary Data 2–6. A large number of differentially expressed (DE) PAC APA genes were identified. Compared to cotyledon, over 1,000 DE-PAC APA genes were differentially utilized in 2 TL (Figure 3A). In the true leaf comparisons, the largest number of DE-PAC APA genes was found between 2 TL and 18TL, possibly because these true leaves were the farthest apart. Interestingly, the number of DE-PAC APA genes showed a reduction trend with adjacent leaf comparison (Figure 3B). DE genes appeared to show a similar tendency to DE-PAC APA genes (Figure 3C,D). These results indicated that APA might be involved in Arabidopsis leaf development.

FIGURE 3

FIGURE 3. Analysis of differentially expressed poly(A) cluster (DE-PAC) APA genes and DE genes. Differential expression of PACs and genes were identified using the DESeq2 package. PAT read differences with padj < 0.05 were selected as DE genes and DE-PAC. (A) DE-PAC APA genes in pairwise comparisons. (B) The number of DE-PAC APA genes in adjacent true leaf comparison. (C) DE genes in pairwise comparisons. (D) The number of DE genes in adjacent true leaf comparison.

To assess the possible impact of APA on the expression of their associated genes, overlap analyses of DE-PAC APA genes and DE genes were explored (Figure 4). A large proportion of overlap (mostly >60%) of DE genes belong to DE-PAC APA genes in the adjacent leaf comparison and between 2 TL and 18 TL (Figure 4A). These results suggest that APA might play a role in modulating the gene expression, particularly in early- and late-developing leaves.

FIGURE 4

FIGURE 4. Venn diagram of DE-PAC APA genes and DE genes (A) and DE genes versus switch genes (B). The blue circle represents DE genes; the orange circle represents DE-PAC APA genes; and the pink one represents switch genes. The number in the circles shows gene counts. The percentages represent the ratio of DE genes belonging to DE-PAC APA genes or the ratio of switch genes belonging to DE genes.

In order to exclude the effect of transcription on transcript level, poly(A) site usage (PAU) was applied to calculate switch genes. We found that over 200 genes switched APA site usage during leaf development (Figure 4B), except the pair of 18 TL and 14 TL comparisons that has only 72 switch genes which could be because they are both the newest leaves. In addition, over 40% switch genes altered their expression. More interestingly, over 50% of them regulated the transcript shifts and did not change the level of gene expression. This result further indicated that APA might mediate the gene expression involved in leaf development, regardless of the steady-state transcript level.

Specific APA Genes Associated With Metabolic Processes and the Response of Hormones

The abovementioned results showed that the expression of over 1,000 APA genes was altered. Overlap analyses were performed to reveal the specific DE-PAC APA genes among adjacent leaf comparisons. A total of 34 DE-PAC APA genes were found in all comparison groups, termed as common genes (Supplementary Figure S4). Interestingly, we also found hundreds of specific APA genes (Supplementary Figure S4). Gene ontology (GO) analyses were performed to explore the functions of specific APA genes. The top 10 GO terms for each pairwise comparison are represented in Figure 5. These include metabolic processes, cellular processes, and response to the stimulus. We found an ATP-binding microtubule motor family protein, which is related to the leaf size (Fleury et al., 2007). Distal transcript (PA2) abundance has an increased tendency with a new leaf growing up, but no expression in cotyledon (Figure 6A). Conversely, another tetratricopeptide repeat (TPR)-like superfamily protein, involving pathogen response and cell cycle (Ascencio-Ibanez et al., 2008), showed a proximal transcript (PA1) expression level that decreased with leaf development (Figure 6B). However, this gene was not expressed in the 18th leaf since it was not found in the PAT-seq data of this leaf stage.

FIGURE 5

FIGURE 5. GO enrichment analyses of specific DE-PAC APA genes. GO analyses were carried out using agriGO (v.2.0). The top 10 GO terms of the biological process are shown. Solid line: FDR = 0.01. The colors represent different adjacent leaf comparisons, and the sizes of the circles represent gene numbers.

FIGURE 6

FIGURE 6. Sequencing coverage of genes visualized using integrative genomics viewer software. For the gene models on the bottom of each panel, the boxes denote exons, thin boxes denote untranslated regions, lines denote introns, and arrows denote transcription directions. PA indicates alternative polyadenylation site. (A) An ATP-binding microtubule motor family protein (AT4G38950). (B) tetratricopeptide repeat (TPR)-like superfamily protein (AT3G26580). (C) THESEUS1 (THE1) (AT5G54380), a receptor kinase. (D) SAUR16 (AT4G38860). (E) APUM23 (AT1G72320). (F) ARF2 (AT5G62000).

Apart from the common GO terms, some specific GO terms were also detected. The numbers of specific GO terms were 4, 4, 3, 8, and 10 among 2 TL vs. cotyledon, 6 TL vs. 2TL, 10 TL vs. 6TL, 14 TL vs. 10TL, and 18 TL vs. 14 TL pairwise comparisons, respectively, suggesting that DE-PAC APA genes possess diversified functions at different leaf development stages. Notably, one of the specific GO terms among 2 TL vs. cotyledon was related to the response of hormones (Figure 5). Under this term, the expression of THESEUS1 (THE1) encoding a receptor kinase, is regulated by BR and is required for cell elongation (Fujikura et al., 2014). Relative to the cotyledon, the distal transcript (PA2) abundance was greater in 2TL, while the proximal transcripts (PA1) were lower (Figure 6C). However, the expression of PA1 and PA2 was not obviously different among the true leaves. Another gene, small auxin upregulated RNA 16 (SAUR16), a SAUR-like auxin-responsive protein family, contributes to light-induced cotyledon opening (Dong et al., 2019). Notably, the PA1 and PA2 levels, particularly PA2, increased dramatically in the 2 TL relative to that of the cotyledon. In addition, the PA1 and PA2 abundance showed a downward trend during true leaf development (Figure 6D). While the expression level of SAUR16 was in a dynamic mode, the differential APA could aid this change. Importantly, we found that the PA2 expression of Arabidopsis Pumilio (APUM23), functioning in leaf development (Huang et al., 2014), had a downward trend (Figure 6E). In addition, the PA1 abundance of auxin response factor 2 (ARF2) involving in leaf shape (Okushima et al., 2005), also showed a downward trend during leaf development (Figure 6F). These results demonstrated that APA assists to regulate the expression of key genes associated with leaf development.

3′UTR Tends to Shorten During True Leaf Development

As shown in Supplementary Figure S3, majority of the APA are located in 3′UTR, reflecting that 3′UTR APA may be more significant in regulating the gene expression. Thus, we focused on 3′UTR length analyses in different leaf comparisons. Interestingly, more genes used longer 3′UTR in 2 TL than cotyledon (Figure 7A), indicating a preference to use distal poly(A) sites in 3′UTR in 2 TL. However, compared to the 2TL, more genes used shorter 3′UTR in the 18TL, and the ratio of genes using proximal sites to those using distal ones is over five times (Figure 7A). Interestingly, during true leaf development, there is a tendency for more genes to use proximal sites in 3′UTR. Strangely, the 6 TL and 10 TL comparison groups showed an opposite trend (Figure 7B).

FIGURE 7

FIGURE 7. Analysis of 3′UTR length distributions. The p value < 0.01 was defined as the threshold for significance. (A) 3′UTR length changes between 2 TL and cotyledon, 2 TL and 18 TL. (B) 3′UTR length changes during true leaf development. (C) Median weighted length of 3′UTR in significantly lengthened genes in different comparison groups. (D) Median weighted length of 3′UTR in significantly shortened genes in different comparison groups.

The median 3′UTR length of significantly lengthened and shortened genes were analyzed. In significantly lengthened genes of 3′UTR, the average 3′UTR length was 26 nucleotides longer in 2 TL than that in the cotyledon (Figure 7C). The most significant change in average 3′UTR length occurred between 2 TL and 18 TL, with an extension of 65 nucleotides. Similarly, in significantly shortened genes of 3′UTR, 32 and 26 nucleotides were shorter in 2 TL vs. cotyledon and 2 TL vs. 18 TL, respectively (Figure 7D). These results suggested that 3′UTR APA is associated with leaf development.

Previous studies demonstrated that shortening of 3′UTR in mammalian cells results in an increase in protein production, with the exception of miRNA-targeting genes (Sandberg et al., 2008). However, we found that there is no obvious correlation between the 3′UTR length switching and gene expression during leaf development, which is shown in Supplementary Figure S5 where the trend lines are all flat, indicative of no correlation.

Discussion

The leaf is a basic organ that captures light energy and synthesizes organic compounds through photosynthesis. Leaves are the ultimate source of most food. Thus, an understanding of the underlying mechanism of leaf development will provide an insight into the basic processes in biology, and it may have a profound significance regarding strategies to improve the crop performance. Recent data indicate that some genes show a specific trend, either gradually increasing or decreasing during leaf development, indicating that leaf ontogeny is continuously changing, involving sequential and transient activation of regulators (Schmid et al., 2005; Efroni et al., 2008). However, the underlying molecular mechanism of leaf development remains to be elucidated.

Polyadenylation is an important co-transcription and mediates a set of biological processes (Liu et al., 2014). Previous research showed that APA of Ethylene Response Factor4 (ERF4) is involved in the regulation of Arabidopsis leaf senescence (Riester et al., 2019). In our data, we found that the proximal transcript (PA1) of ERF4 had a downward trend with a new leaf growing up (Supplementary Figure S6A). However, the distal transcript (PA2) had an anomaly change. In addition, Hydroxypyruvate Reductase (HPR), encoding a leaf-peroxisomal protein, has two transcripts by alternative splicing (Mano et al., 2000). Interestingly, we found that HPR contained different poly(A) sites, and the expression level of primary transcript (PA1) showed increased tendency during leaf development, suggesting that APA and alternative splicing may show a crosstalk relationship during leaf development.

In this study, we detected a dramatic change of genome-wide APA during Arabidopsis leaf development and found that the APA ratio, DE-PAC APA genes, and DE genes show a tendency to be reduced during leaf development, suggesting that APA might be involved in leaf development. In addition, APA is tightly associated with the expression of DE genes during leaf aging. Significant changes in 3′UTR APA occur during leaf development. Compared to cotyledon, the second true leaf prefers to use the distal poly(A) site in 3′UTR. However, an opposing trend is seen in terms of 3′UTR APA during true leaf development. Strangely, compared to 6TL, more genes in 10 TL showed a preference to use distal poly(A) site in 3′UTR, indicating a complex mechanism in addition to APA that is involved in leaf development.

Previous research showed that the deficiency of CFIm 68 and CFIm25 leads to a shift toward proximal poly(A) site usage in human HEK293 cells (Gruber et al., 2012; Martin et al., 2012). Interestingly, we also found that the expression level of CFIS1 and CFIS2, orthologous to CFIm 68 and CFIm25, respectively, had a downward trend during leaf development (Supplementary Figure S7A,B), indicating that the expression of CFIS1 and CFIS2 may contribute to the usage of poly(A) site in 3′UTR during leaf development. On the contrary, the loss of Pcf11 enhances distal poly(A) site usages in 3′UTR (Li et al., 2015). We also found that the expression of PCFS1 and PCFS5, the Arabidopsis orthologous of hPcf11 complex, showed an increased tendency during leaf development, except in 14 TL and 18 TL (Supplementary Figure S7A,B).

Previous studies have found that the length of 3′UTR modulates the expression level of genes in mammals (Sandberg et al., 2008). While there are about 40% of switch genes associated with the change of steady-state levels of their corresponding total transcripts, many other APA switch genes (>50%) did not show co-relationship with the steady-state level of transcripts. For the latter case, at the total transcript level, there is no obvious or little correlation between the length switching of 3′UTR (switch genes) and the steady-state gene expression levels (transcript copy numbers of genes) during leaf development (Supplementary Figure S5). This result points to a scenario where, in >50% of genes, different lengths of the APA switched transcripts may impact the translation efficiency, not the stability of these genes. Of course, this is a general trend line, not representing every gene, as many examples can be found in Supplementary Figure S5 where the data points are above or under the trend line. In fact, there are few miRNA targets that can be found between two APA sites in Arabidopsis (QQLi, unpublished observation) to justify the action of miRNA in mRNA degradations. On the other hand, what proportion of these hundreds of APA genes in fact participates in leaf development is an open question. It requires further studies to sort this out. In addition, our previous work also did not find an obvious negative correlation between 3′UTR length variation and gene expressions in different rice (Oryza sativa) tissues and different poly(A) factor mutants (Shen et al., 2008; Yu et al., 2019). Our results confirm that there may exist different mechanisms between what happened in plants and in mammals with respect to the correlation of 3′UTR length variation and the level of gene expression.

Leaf morphogenesis is an active process, controlled by plant hormones, transcriptional regulators, miRNAs, and mechanical forces, among factors (Wang et al., 2021). Significantly, we found that some key APA genes are related to leaf development. Interestingly, the distal transcript (PA2) abundance of an ATP-binding microtubule motor family protein, which is related to leaf size (Fleury et al., 2007), has an increased tendency toward growth with the appearance of a new leaf. The size and shape of a leaf are related to cell proliferation and differentiation. We also found that the proximal transcript (PA1) expression level of tetratricopeptide repeat (TPR)-like superfamily protein, involved in the cell cycle (Ascencio-Ibanez et al., 2008), decreased along with leaf development. In addition, the accumulating evidence demonstrated that plant hormones, including auxin, cytokinin, and BR have important roles during leaf development. GO enrichment analyses of the specific DE-PAC APA gene were carried out in this study. One of the terms, response to hormone, was found between the cotyledon and the 2 TL comparisons. SAUR16, a SAUR-like auxin-responsive protein family, contributes to light-Induced cotyledon opening (Dong et al., 2019). Notably, the PA1 and PA2 levels, particularly PA2, increased dramatically in the 2 TL relative to the levels in the cotyledon. In addition, PA1 and PA2 abundance show a downward trend during true leaf development. These results indicate that APA regulated the expression of key genes associated with leaf development. Although some genes have been identified as participating in leaf development in recent years, the molecular mechanism is not fully clear. With the development of genomic approaches, it will help to bring insights into mechanisms of the dramatic leaf development process.

Data Availability Statement

The datasets presented in this study can be found in online repositories. Raw high-throughput sequence data for this study were deposited in the National Center for Biotechnology Information SRA database under the code PRJNA778569.

Author Contributions

QQL and LH designed the research. ZY performed experiments. QQL and ZY contributed to the data analyses and wrote and revised the manuscript. All authors approved the manuscript.

Funding

This work was supported in part by a grant from the National Key R&D Project of China (2016YFE0108800).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors, and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Acknowledgments

The authors thank Haidong Qu, Xiuxiu Wang, and Xiaoxuan Zhou for their technical assistance and Charlesworth Author Services for language editing.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.2022.863253/full#supplementary-material

References

Ascencio-Ibáñez, J. T., Sozzani, R., Lee, T.-J., Chu, T.-M., Wolfinger, R. D., Cella, R., et al. (2008). Global Analysis of Arabidopsis Gene Expression Uncovers a Complex Array of Changes Impacting Pathogen Response and Cell Cycle during Geminivirus Infection. Plant Physiol. 148, 436–454. doi:10.1104/pp.108.121038

PubMed Abstract | CrossRef Full Text | Google Scholar

Bar, M., and Ori, N. (2014). Leaf Development and Morphogenesis. Development 141, 4219–4230. doi:10.1242/dev.106195

PubMed Abstract | CrossRef Full Text | Google Scholar

Bruggeman, Q., Garmier, M., De Bont, L., Soubigou-Taconnat, L., Mazubert, C., Benhamed, M., et al. (2014). The Polyadenylation Factor Subunit Cleavage and Polyadenylation Specificity Factor30: A Key Factor of Programmed Cell Death and a Regulator of Immunity in Arabidopsis. Plant Physiol. 165, 732–746. doi:10.1104/pp.114.236083

PubMed Abstract | CrossRef Full Text | Google Scholar

Carpenter, S., Ricci, E. P., Mercier, B. C., Moore, M. J., and Fitzgerald, K. A. (2014). Post-transcriptional Regulation of Gene Expression in Innate Immunity. Nat. Rev. Immunol. 14, 361–376. doi:10.1038/nri3682

PubMed Abstract | CrossRef Full Text | Google Scholar

Cyrek, M., Fedak, H., Ciesielski, A., Guo, Y., Sliwa, A., Brzezniak, L., et al. (2016). Seed Dormancy in Arabidopsis Is Controlled by Alternative Polyadenylation of DOG1. Plant Physiol. 170, 947–955. doi:10.1104/pp.15.01483

PubMed Abstract | CrossRef Full Text | Google Scholar

Derti, A., Garrett-Engele, P., Macisaac, K. D., Stevens, R. C., Sriram, S., Chen, R., et al. (2012). A Quantitative Atlas of Polyadenylation in Five Mammals. Genome Res. 22, 1173–1183. doi:10.1101/gr.132563.111

PubMed Abstract | CrossRef Full Text | Google Scholar

Dong, J., Sun, N., Yang, J., Deng, Z., Lan, J., Qin, G., et al. (2019). The Transcription Factors TCP4 and PIF3 Antagonistically Regulate Organ-specific Light Induction of SAUR Genes to Modulate Cotyledon Opening during De-etiolation in Arabidopsis. Plant Cell 31, 1155–1170. doi:10.1105/tpc.18.00803

PubMed Abstract | CrossRef Full Text | Google Scholar

Du, F., Guan, C., and Jiao, Y. (2018). Molecular Mechanisms of Leaf Morphogenesis. Mol. Plant 11, 1117–1134. doi:10.1016/j.molp.2018.06.006

PubMed Abstract | CrossRef Full Text | Google Scholar

Efroni, I., Blum, E., Goldshmidt, A., and Eshed, Y. (2008). A Protracted and Dynamic Maturation Schedule UnderliesArabidopsisLeaf Development. Plant Cell 20, 2293–2306. doi:10.1105/tpc.107.057521

PubMed Abstract | CrossRef Full Text | Google Scholar

Elkon, R., Ugalde, A. P., and Agami, R. (2013). Alternative Cleavage and Polyadenylation: Extent, Regulation and Function. Nat. Rev. Genet. 14, 496–506. doi:10.1038/nrg3482

PubMed Abstract | CrossRef Full Text | Google Scholar

Fleury, D., Himanen, K., Cnops, G., Nelissen, H., Boccardi, T. M., Maere, S., et al. (2007). TheArabidopsis thalianaHomolog of YeastBRE1Has a Function in Cell Cycle Regulation during Early Leaf and Root Growth. Plant Cell 19, 417–432. doi:10.1105/tpc.106.041319

PubMed Abstract | CrossRef Full Text | Google Scholar

Fu, H., Yang, D., Su, W., Ma, L., Shen, Y., Ji, G., et al. (2016). Genome-wide Dynamics of Alternative Polyadenylation in rice. Genome Res. 26, 1753–1760. doi:10.1101/gr.210757.116

PubMed Abstract | CrossRef Full Text | Google Scholar

Fujikura, U., Elsaesser, L., Breuninger, H., Sánchez-Rodríguez, C., Ivakov, A., Laux, T., et al. (2014). Atkinesin-13A Modulates Cell-wall Synthesis and Cell Expansion in Arabidopsis thaliana via the THESEUS1 Pathway. Plos Genet. 10, e1004627. doi:10.1371/journal.pgen.1004627

PubMed Abstract | CrossRef Full Text | Google Scholar

Gruber, A. R., Martin, G., Keller, W., and Zavolan, M. (2012). Cleavage Factor Imis a Key Regulator of 3′ UTR Length. RNA Biol. 9, 1405–1412. doi:10.4161/rna.22570

PubMed Abstract | CrossRef Full Text | Google Scholar

Guan, C., Du, F., Xiong, Y., and Jiao, Y. (2019). The 35S Promoter‐driven mDII Auxin Control Sensor Is Uniformly Distributed in Leaf Primordia. J. Integr. Plant Biol. 61, 1114–1120. doi:10.1111/jipb.12853

PubMed Abstract | CrossRef Full Text | Google Scholar

Ha, C. M., Kim, G.-T., Kim, B. C., Jun, J. H., Soh, M. S., Ueno, Y., et al. (2003). The BLADE-ON-PETIOLE 1 Gene Controls Leaf Pattern Formation through the Modulation of Meristematic Activity in Arabidopsis. Development 130, 161–172. doi:10.1242/dev.00196

PubMed Abstract | CrossRef Full Text | Google Scholar

Ha, K. C. H., Blencowe, B. J., and Morris, Q. (2018). QAPA: a New Method for the Systematic Analysis of Alternative Polyadenylation from RNA-Seq Data. Genome Biol. 19, 45. doi:10.1186/s13059-018-1414-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Huang, T., Kerstetter, R. A., and Irish, V. F. (2014). APUM23, a PUF Family Protein, Functions in Leaf Development and Organ Polarity in Arabidopsis. J. Exp. Bot. 65, 1181–1191. doi:10.1093/jxb/ert478

PubMed Abstract | CrossRef Full Text | Google Scholar

Jun, J. H., Ha, C. M., and Fletcher, J. C. (2010). BLADE-ON-PETIOLE1 Coordinates Organ Determinacy and Axial Polarity inArabidopsisby Directly Activating Asymmetric Leaves2. Plant Cell 22, 62–76. doi:10.1105/tpc.109.070763

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, W., You, B., Hoque, M., Zheng, D., Luo, W., Ji, Z., et al. (2015). Systematic Profiling of Poly(A)+ Transcripts Modulated by Core 3' End Processing and Splicing Factors Reveals Regulatory Rules of Alternative Cleavage and Polyadenylation. Plos Genet. 11, e1005166. doi:10.1371/journal.pgen.1005166

PubMed Abstract | CrossRef Full Text | Google Scholar

Lin, J., Hung, F.-Y., Ye, C., Hong, L., Shih, Y.-H., Wu, K., et al. (2020). HDA6-dependent Histone Deacetylation Regulates mRNA Polyadenylation in Arabidopsis. Genome Res. 30, 1407–1417. doi:10.1101/gr.255232.119

PubMed Abstract | CrossRef Full Text | Google Scholar

Lin, J., Xu, R., Wu, X., Shen, Y., and Li, Q. Q. (2017). Role of Cleavage and Polyadenylation Specificity Factor 100: Anchoring Poly(A) Sites and Modulating Transcription Termination. Plant J. 91, 829–839. doi:10.1111/tpj.13611

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, M., Xu, R., Merrill, C., Hong, L., Von Lanken, C., Hunt, A. G., et al. (2014). Integration of Developmental and Environmental Signals via a Polyadenylation Factor in Arabidopsis. PLoS One 9, e115779. doi:10.1371/journal.pone.0115779

PubMed Abstract | CrossRef Full Text | Google Scholar

Long, J. A., Moan, E. I., Medford, J. I., and Barton, M. K. (1996). A Member of the KNOTTED Class of Homeodomain Proteins Encoded by the STM Gene of Arabidopsis. Nature 379, 66–69. doi:10.1038/379066a0

PubMed Abstract | CrossRef Full Text | Google Scholar

Mano, S., Hayashi, M., and Nishimura, M. (2000). A Leaf-Peroxisomal Protein, Hydroxypyruvate Reductase, Is Produced by Light-Regulated Alternative Splicing. Cbb 32, 147–154. doi:10.1385/cbb:32:1-3:147

CrossRef Full Text | Google Scholar

Manuela, D., and Xu, M. (2020). Patterning a Leaf by Establishing Polarities. Front. Plant Sci. 11, 568730. doi:10.3389/fpls.2020.568730

PubMed Abstract | CrossRef Full Text | Google Scholar

Martin, G., Gruber, A. R., Keller, W., and Zavolan, M. (2012). Genome-wide Analysis of Pre-mRNA 3′ End Processing Reveals a Decisive Role of Human Cleavage Factor I in the Regulation of 3′ UTR Length. Cel Rep. 1, 753–763. doi:10.1016/j.celrep.2012.05.003

CrossRef Full Text | Google Scholar

Mayr, C., and Bartel, D. P. (2009). Widespread Shortening of 3′UTRs by Alternative Cleavage and Polyadenylation Activates Oncogenes in Cancer Cells. Cell 138, 673–684. doi:10.1016/j.cell.2009.06.016

PubMed Abstract | CrossRef Full Text | Google Scholar

Mitchum, M. G., Yamaguchi, S., Hanada, A., Kuwahara, A., Yoshioka, Y., Kato, T., et al. (2006). Distinct and Overlapping Roles of Two Gibberellin 3-oxidases in Arabidopsis Development. Plant J. 45, 804–818. doi:10.1111/j.1365-313X.2005.02642.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Mizukami, Y., and Fischer, R. L. (2000). Plant Organ Size Control: AINTEGUMENTA Regulates Growth and Cell Numbers during Organogenesis. Proc. Natl. Acad. Sci. U.S.A. 97, 942–947. doi:10.1073/pnas.97.2.942

PubMed Abstract | CrossRef Full Text | Google Scholar

Nakata, M., Matsumoto, N., Tsugeki, R., Rikirsch, E., Laux, T., and Okada, K. (2012). Roles of the Middle Domain-specific Wuschel-Related Homeobox Genes in Early Development of Leaves in Arabidopsis. Plant Cell 24, 519–535. doi:10.1105/tpc.111.092858

PubMed Abstract | CrossRef Full Text | Google Scholar

Oh, M.-H., Sun, J., Oh, D. H., Zielinski, R. E., Clouse, S. D., and Huber, S. C. (2011). Enhancing Arabidopsis Leaf Growth by Engineering the Brassinosteroid Insensitive1 Receptor Kinase. Plant Physiol. 157, 120–131. doi:10.1104/pp.111.182741

PubMed Abstract | CrossRef Full Text | Google Scholar

Okushima, Y., Mitina, I., Quach, H. L., and Theologis, A. (2005). Auxin Response Factor 2 (ARF2): a Pleiotropic Developmental Regulator. Plant J. 43, 29–46. doi:10.1111/j.1365-313X.2005.02426.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Poethig, R. S. (1997). Leaf Morphogenesis in Flowering Plants. Plant Cell 9, 1077–1087. doi:10.1105/tpc.9.7.1077

PubMed Abstract | CrossRef Full Text | Google Scholar

Riester, L., Köster-Hofmann, S., Doll, J., Berendzen, K., and Zentgraf, U. (2019). Impact of Alternatively Polyadenylated Isoforms of ETHYLENE RESPONSE FACTOR4 with Activator and Repressor Function on Senescence in Arabidopsis thaliana L. Genes 10, 91. doi:10.3390/genes10020091

PubMed Abstract | CrossRef Full Text | Google Scholar

Sandberg, R., Neilson, J. R., Sarma, A., Sharp, P. A., and Burge, C. B. (2008). Proliferating Cells Express mRNAs with Shortened 3' Untranslated Regions and Fewer microRNA Target Sites. Science 320, 1643–1647. doi:10.1126/science.1155390

PubMed Abstract | CrossRef Full Text | Google Scholar

Sarojam, R., Sappl, P. G., Goldshmidt, A., Efroni, I., Floyd, S. K., Eshed, Y., et al. (2010). Differentiating Arabidopsis Shoots from Leaves by Combined YABBY Activities. Plant Cell 22, 2113–2130. doi:10.1105/tpc.110.075853

PubMed Abstract | CrossRef Full Text | Google Scholar

Schmid, M., Davison, T. S., Henz, S. R., Pape, U. J., Demar, M., Vingron, M., et al. (2005). A Gene Expression Map of Arabidopsis thaliana Development. Nat. Genet. 37, 501–506. doi:10.1038/ng1543

PubMed Abstract | CrossRef Full Text | Google Scholar

Shen, Y., Ji, G., Haas, B. J., Wu, X., Zheng, J., Reese, G. J., et al. (2008). Genome Level Analysis of rice mRNA 3′-end Processing Signals and Alternative Polyadenylation. Nucleic Acids Res. 36, 3150–3161. doi:10.1093/nar/gkn158

PubMed Abstract | CrossRef Full Text | Google Scholar

Thomas, P. E., Wu, X., Liu, M., Gaffney, B., Ji, G., Li, Q. Q., et al. (2012). Genome-wide Control of Polyadenylation Site Choice by CPSF30 in Arabidopsis. Plant Cell 24, 4376–4388. doi:10.1105/tpc.112.096107

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, H., Kong, F., and Zhou, C. (2021). From Genes to Networks: The Genetic Control of Leaf Development. J. Integr. Plant Biol. 63, 1181–1196. doi:10.1111/jipb.13084

PubMed Abstract | CrossRef Full Text | Google Scholar

Wu, X., Liu, M., Downie, B., Liang, C., Ji, G., Li, Q. Q., et al. (2011). Genome-wide Landscape of Polyadenylation in Arabidopsis Provides Evidence for Extensive Alternative Polyadenylation. Proc. Natl. Acad. Sci. U.S.A. 108, 12533–12538. doi:10.1073/pnas.1019732108

PubMed Abstract | CrossRef Full Text | Google Scholar

Xing, D., and Li, Q. Q. (2011). Alternative Polyadenylation and Gene Expression Regulation in Plants. Wiley Interdiscip Rev RNA 2, 445–458. doi:10.1002/wrna.59

PubMed Abstract | CrossRef Full Text | Google Scholar

Yu, Z., Lin, J., and Li, Q. Q. (2019). Transcriptome Analyses of FY Mutants Reveal its Role in mRNA Alternative Polyadenylation. Plant Cell 31, 2332–2352. doi:10.1105/tpc.18.00545

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, J., Addepalli, B., Yun, K.-Y., Hunt, A. G., Xu, R., Rao, S., et al. (2008). A Polyadenylation Factor Subunit Implicated in Regulating Oxidative Signaling in Arabidopsis thaliana. PLoS One 3, e2410. doi:10.1371/journal.pone.0002410

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, Y., Gu, L., Hou, Y., Wang, L., Deng, X., Hang, R., et al. (2015). Integrative Genome-wide Analysis Reveals HLP1, a Novel RNA-Binding Protein, Regulates Plant Flowering by Targeting Alternative Polyadenylation. Cell Res 25, 864–876. doi:10.1038/cr.2015.77

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhiponova, M. K., Vanhoutte, I., Boudolf, V., Betti, C., Dhondt, S., Coppens, F., et al. (2013). Brassinosteroid Production and Signaling Differentially Control Cell Division and Expansion in the Leaf. New Phytol. 197, 490–502. doi:10.1111/nph.12036

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhu, H., Hu, F., Wang, R., Zhou, X., Sze, S.-H., Liou, L. W., et al. (2011). Arabidopsis Argonaute10 Specifically Sequesters miR166/165 to Regulate Shoot Apical Meristem Development. Cell 145, 242–256. doi:10.1016/j.cell.2011.03.024

PubMed Abstract | CrossRef Full Text | Google Scholar