Rat primary aortic endothelial cells (RAECs) were isolated from male rats and those with the purity greater than 95% (determined by both positive staining for CD31 and morphologically based on their classical “hill and valley” appearance) were kept for the subsequent study, as performed according to the methods described previously (Yang et al., 2018). Cells that experienced three to five passages were used for Ang II stimulation experiments. Quantified RAECs were treated with Ang II (2 μM) or vehicle (sodium chloride) for 48 hours (h) as described before (Yang et al., 2018).

PolyA⁺ RNA was enriched by oligo (dT)25 Dynabeads (Invitrogen) from total RNA. The dUTP-based strand-specific RNA-seq libraries for control and Ang II-induced RAECs were constructed according to the method previously described (Yang et al., 2020). RNA-seq libraries were sequenced on Illumina HiSeq platform, and paired-end reads of 150 nucleotide (nt) lengths were obtained.

PolyA⁺ RNA after double-stranded cDNA preparation, PCR amplification, product purification with AMPurePB magnetic beads, verification of DNA size, followed by qualitative and quantitative analysis using a Bioanalyzer. The resulting double-stranded DNA was used to generate the SMRTbell™ libraries using PacBio Template Prep Kit. SMRTbell templates are then sequenced in the PacBio Sequel System.

We performed initial data processing using IsoSeq v3 pipeline (IsoSeq/isoseq-clustering.md at master·PacificBiosciences/IsoSeq·GitHub) for raw subreads. The circular consensus sequencing (CCS) reads were obtained through the consensus generation step. Then, we determine whether they are classified as full-length reads based on the presence or absence of the 5′and 3′primers, and use lima v2.0.0 (https://lima.how/) to remove the primers. The obtained full-length reads, were refined by trimming poly (A) tails and removing concatemer, then full-length non-concatemer (FLNC) reads are finally obtained. To find transcript clusters, isoform-level clustering algorithm ICE (Iterative Clustering for Error Correction) (Gordon et al., 2015) was applied to all full-length (FL) transcripts to obtain the consensus sequence for each cluster. Quiver (Chin et al., 2013) was used for error correction to obtain high-quality (HQ) (accuracy, ≥99%) isoforms. To get non-redundant isoforms, HQ isoforms were mapped to rat genome (mRatBN7.2) using minimap2 (Li, 2018) with parameters -ax splice -uf--secondary = no -C5 -O6,24 -B4, and collapsed by a third-party module Cupcakes Python script (collapse_isoforms_by_sam.py) with the following parameters “-c 0.99 -i 0.85”. The 5′degraded isoforms were filtered away because the Clontech SMARTer cDNA kit used to create the full-length cDNA does not do cap trap, these may be 5′RNA degradation products and not biologically meaningful. In that case, we filter them away according to the method described in Cupcake (https://github.com/Magdoll/cDNA_Cupcake/wiki/Cupcake:-supporting-scripts-for-Iso-Seq-after-clustering-step#filtersubset). Then the resulted isoforms were finally annotated using SQANTI3 (Tardaguila et al., 2018). SQANTI3 marked the internal priming resulted PAS by identifying downstream six consecutive “A” or the percentage of “A” is greater than 60% in the 20 nucleotides downstream of the potential PAS, the transcript is considered to be an IntraPriming transcript and is filtered out.

The coding potential of the resulted isoforms was also predicted with SQANTI3 which uses GeneMarkS-T (GMST) algorithm (Tang et al., 2015). GeneMarkS-T utilizes iterative self-training and a hidden semi-Markov model to predict coding regions in eukaryotic transcripts.

All FLNC reads were merged and the FMLRC (Wang J. R. et al., 2018) was then used for error correction with RNA-seq to get high-quality FLNC reads. The corrected FLNC reads were then mapped to mRatBN7.2 genome using minimap2 (Li, 2018) with parameters -ax splice -uf--secondary = no -C5 -O6,24 -B4, then using the sam_to_gff3.py from cDNA_Cupcake to generate an FLNC gff3 file. This file was then compared with non-redundant HQ isoforms using inhouse scripts to calculate the FLNC reads supporting HQ isoforms in Control and Ang II samples independently.

First, RNA-seq raw reads were cleaned by trimming using Cutadapt (Martin, 2011) with the default parameter. The clean reads were mapped to the mRatBN7.2 reference genome using STAR (version 2.7.7a) (Dobin et al., 2013) with parameters--twopassMode Basic--outFilterMultimapNmax 1 --outSAMstrandField intronMotif. Differential expression analysis between two conditions was performed using the DESeq2 R package (version 1.30.1) (Love et al., 2014). The resulting p values were adjusted using Benjamin and Hochberg’s approach for controlling for the false discovery rate (FDR). Genes with an adjusted p-value <0.05, and |FoldChange (FC)| > 1.5 identified by DESeq2 (Love et al., 2014) were assigned as differentially expressed. Transcripts per million (TPM) of reference transcriptome defined using Iso-Seq was used in calculating the expression level of transcripts by Salmon (v1.8.0) (Patro et al., 2017).

For differentially expressed genes (DEGs) and differentially spliced genes, GO and KEGG enrichment analysis were performed with clusterProfiler (Yu et al., 2012), which supports statistical analysis and visualization of functional profiles for genes and gene clusters.

Active promoter sites were identified by proActiv (v0.1.0) (Demircioğlu et al., 2019) based on transcriptome generated by Iso-Seq. Short-reads of RNA-seq were aligned to the reference genome using STAR (Dobin et al., 2013). The aligned RNA-seq bam files and GTF files generated from Iso-Seq were used as input by proActive to estimate promoter activities in each sample. We divided the promoter set into three categories depending on their relative promoter activity calculated by proActive, major, minor and inactive promoters. Promoters with the highest relative activity for each gene across the sample cohort were defined as major promoters. Besides, promoters with relative promoter activity less than 0.1 constitute inactive promoters whereas the rest promoters of a gene constitute minor promoters. Known promoter was defined if known genome annotation also uses this promoter. Novel promoter was defined only if Iso-Seq isoforms use this promoter. Differentially regulated promoters (DRPs) were identified using DESeq2. The resultant p values were adjusted using BH approach for controlling the false discovery rate (FDR). The promoters with an activity change level of |FoldChange| > 1.5 and adjusted p-value <0.05 were considered significant DRPs. We identified alternative promoters (APs) by considering both promoter expression and promoter activity. The criteria were as follows: (1) mean relative promoter activity >0.1 and were significantly changed between two conditions (FDR <0.05, BH adjusted); (2) |Fold Change| > 1.5 of promoter expression. (3) Gene expression FPKM >1.

Through the use of SUPPA2 (Trincado et al., 2018), a tool designed to find splicing events and determine whether there are differences between samples, alternative splicing patterns were found with the PacBio full-length transcripts. By scanning the ‘exon’ lines and categorizing each local event as a particular splicing type, the SUPPA2 “generateEvents” command was used to create local alternative splicing events from our high-quality PacBio transcript annotation. Then, using the Salmon (Patro et al., 2017) quantification results, the SUPPA2 “psiPerEvent” command was used to quantify the percentage of spliced-in (PSI) values for each local alternative splicing event. Differential expression of AS events was screened out with the adjusted p-value <0.05.

PAS cluster was using with modified TAPIS 1.2.1 (Abdel-Ghany et al., 2016) (https://bitbucket.org/comp_bio/tapis/src/master/) scripts that PASs within 50 base pair (bp) were clustered as one PAS. APA analysis between two conditions was performed by using the algorithm of DaPars (Xia et al., 2014). As DaPars predicts the proximal PAS using RNA-seq and Iso-Seq could provide validated PAS, we therefore supplied Iso-Seq supported PAS for APA analysis. These PAS were used as prior knowledge using a modified DaPars algorithm to calculate the distal PAS usage index (PDUI) based on RNA-seq coverage, and then Fisher exact test was used to determine differential APA events between the two groups, and the criteria for defining differential APA events between two groups were consistent with DaPars.

3′READS and Poly (A)-seq reads were mapped to mRatBN7.2 genome using Bowtie 2 with the default parameters (Langmead and Salzberg, 2012). Reads start pos were considered as potential PAS. PAS with at least five reads support is considered a potential validated PAS.

CAGE reads were mapped to mRatBN7.2 genome using Bowtie 2 (Langmead and Salzberg, 2012) with default parameters. Then peak calling was performed using MACS2 (Zhang et al., 2008). The resultant files with suffix of.treat_pileup.bdg and.summits.bed were used for further analysis.

RNA-seq is widely used for the quantification of gene expression and AS events. However, it is still challenging to accurately identify full-length splicing isoforms using RNA-seq alone (Abdel-Ghany et al., 2016). To reveal the complexity of the transcriptome of rat primary aortic endothelial cells (RAECs) and discover dynamic changes of alternative promoter (AP), alternative splicing (AS) and alternative polyadenylation (APA) in Ang II-induced endothelial cell senescence, we sequenced the transcriptome of RAECs before (Control) and after Ang II treatment using both PacBio Iso-Seq and Illumina RNA-seq techniques.

We used PacBio’s official pipeline, IsoSeq3, to filter and process the raw sequencing files (details see: https://github.com/PacificBiosciences/IsoSeq/blob/master/isoseq-clustering.md). As a result, a total of 27,668,459 raw subreads with an average length of 1,838 bp were obtained in Control sample and 33,151,284 raw subreads with an average length of 1,700 bp in Ang II treatment samples (Table 1). The PacBio official pipeline we adopted for the analysis of raw data includes circular consensus sequencing (CCS) reads calling, primer removal and demultiplexing, and the refining steps (Supplementary Figure S1) (Details in Methods section). After processing with these steps, we got 1,469,285 FLNC reads with an average length of 2,361 bp in the Control sample and 1,676,143 FLNC with an average length of 2,486 bp in the Ang II sample, respectively (Table 1). Then, polish step was performed on all the full-length transcripts to generate consensus isoforms. It created a consensus isoform for each cluster after classifying them into groups based on their sequence similarity. Subsequently, the full insert sequence that had not been tested was aligned to the adjusted to produce high-quality (HQ) isoforms (125,900) with a mean length of 2,603 from the consensus isoform and low quality (LQ) isoforms (5,323) with a mean length of 2,807 (Table 1). To get non-redundant transcripts, HQ isoforms were mapped to rat genome (mRatBN7.2) using minimap2 (Li, 2018) and collapsed with a third-party module Cupcakes Python script (collapse_isoforms_by_sam.py). Then we got 36,278 transcripts derived from 10,145 gene loci. These transcripts were annotated by SQANTI3 (Tardaguila et al., 2018) for further analysis.

To validate the reliability of these isoforms identified by Iso-Seq, we used a series of public available datasets generated by independent methods including cap analysis gene expression (CAGE) data from FANTOM (Lizio et al., 2015), RNA-seq short junction reads, rat 3′READ raw data used in polyA_DB3 (Wang R. et al., 2018), rat PolyA-Seq data (Derti et al., 2012), and ncbiRefSeq known rat genome annotation (REF) to construct reliable 5′-ends, splicing junction, and 3′-ends pool, respectively (Figure 1Ai). If the isoforms have the same splicing pattern and PAS (within 50 nt) but different potential TSS, we only kept theose whose TSS are closest to the validated 5′-ends pool (Figure 1Aii, iii). Since Iso-Seq protocol uses the similar strategy as the Clontech SMARTer cDNA kit to generate full-length cDNA, which does not perform cap trap, and thus the obtained 5′ends can be possibly derived from degradation products or reverse transcriptase drop offs (Schaarschmidt et al., 2020). To filter out possible artificial isoforms, Iso-Seq defined transcription start site (TSS) and polyadenylation site (PAS) within 50 nt of validated 5′-ends and 3′-ends pool respectively, were called supported. As for splicing, if an Iso-Seq defined splicing event could be supported by at least two short junction reads in the pool, it was also called supported. We found considerable high support rates for TSS (96.5%), splicing junction (86.4%) and PAS (93.0%) (Figure 1B), suggesting the reliability of Iso-Seq detected isoforms.

After confirming the reliability of the detected isoforms, we next asked whether full-length Iso-Seq could improve the gene annotation of rat genome. We used UCSC ncbiRefSeq as a reference for comparison, because this annotation is the most complete for rat genome (Ji et al., 2020). We adopted the transcript structure classification strategy of SQANTI3 (Tardaguila et al., 2018), which divides transcripts into eight categories compared to reference, to classify the Iso-Seq identified transcripts (Supplementary Figure S2). Among all these transcripts annotated to known genes (n = 34,642), 61.59% (n = 21,337) were categorized as full-splice match (FSM, corresponds to a transcript with all the same splice junctions as a reference transcript) in catalog, 15.12% (n = 5,238) were categorized as novel not in catalog (NNC; transcripts use novel donors and/or acceptors), 13.54% (n = 4,692) were categorized as novel in catalog (NIC, transcripts contain new combinations of annotated splice junctions or novel splice junctions formed from annotated donors and acceptors), 8.42% (n = 2,916) were categorized as incomplete splice match (ISM, transcripts matching consecutive, but not all, splicing junctions of the REF), 0.59% (n = 204) were categorized as genic, 0.32% (n = 112) were categorized as intergenic transcripts, 0.23% (n = 80) were categorized as fusion, and 0.18% (n = 63) were categorized as antisense transcripts [structure categories defined by (Tardaguila et al., 2018)] (Figure 1C).

We proceeded to characterize transcripts in these eight structure categories above. First, we cross-validated the reliability of transcripts in each category. PacBio transcripts with 5′-ends within 50 nt to CAGE or REF TSSs peaks were considered as validated TSSs, and PacBio transcripts 3′-end within 50 nt to PASs (generated from PolyA-Seq, 3′READ or REF) were considered as cross-validated PASs. Junctions of PacBio transcripts supported with RNA-seq junction reads were used for percent calculation. In line with our speculation, FSM category has high validation rate regarding TSS (100%), splicing junction (99.87%) and PAS (92.98%) (Supplementary Figure S3). Noteworthy, two categories that contain novel sites (NNC) and (NIC) also showed consistently high validation rates (>95% for TSS, junction and PAS) (Supplementary Figure S3), suggesting the novel isoforms detected by Iso-Seq are reliable. For categories that have lower validation rates (such as junction in genic and antisense categories), experimental or technical errors could possibly account for them, as discussed previously (Tardaguila et al., 2018).

After demonstrating the authenticity of these novel isoforms identified by Iso-Seq, we further explored the potential underlying contributing factors. We found AP, AS and APA could all contribute to the generation of novel isoforms, and APA contributed most, followed by AS and AP (Figure 1D). Further, 92.83% of the APA-derived novel PASs are within 24 nt distance to validated PAS clusters (Supplementary Figure S4), further supporting the biggest contribution of APA to novel isoforms. There were several good examples showing novel isoforms annotated by our data. For instance, the gene Arhgap21 has two novel TSSs not annotated in ncbiRefSeq but supported by CAGE data (Figure 1E). The gene Khdc4 having intron retention (a type of AS), is also a contributor to novel isoforms (Figure 1F). The gene Arhgef2 in the ncbiRefSeq annotation has only one PAS, however, we identified several APA isoforms of this gene based on our PacBio data, with both the shortest (novel proximal PAS) and the longest (annotated distal PAS) isoforms supported by strong PolyA-Seq/3′READs signals (Figure 1G). Together, these above results support the notion that a combination of Iso-Seq and RNA-seq could largely improve the rat genome annotation.

Previous analyses showed that Iso-Seq could identify full-length isoforms qualitatively but less quantitatively. In contrast, RNA-seq is more accurate in quantifying gene expression rather than distinguishing the different isoforms from the same gene (Huang et al., 2021). Having confirmed the improvement of rat genome annotation by the combination of Iso-Seq and RNA-seq, we further applied this strategy to explore the expression changes at both gene and isoform levels in Ang II-induced senescent RAECs. First, we used RNA-seq to perform differential gene expression analysis and identified a total of 2,671 differentially expressed genes (DEGs) between Ang II-induced RAECs and the Control ones, including 1,249 upregulated and 1,422 downregulated DEGs (Supplementary Figure S5A). The downregulated DEGs were most significantly enriched in Cell Cycle and DNA replication (Supplementary Figure S5B), consistent with the cell cycle arrest phenotypes in Ang II-induced senescent RAECs (Yang et al., 2020). The upregulated genes were highly enriched in ECM-receptor interaction, PI3K-Akt signaling pathway and TNF signaling pathway (Supplementary Figure S5C), which were known to promote senescence (Saito and Berk, 2002; Yin et al., 2003; Chen et al., 2012). Notably, a substantial number of DEGs were involved in Cell-Cycle pathway, including Cyclin- Dependent Kinase (CDK) gene family such as Cdk2, encoding a well-known cell cycle regulator, were decreased upon Ang II stimulation (Supplementary Figure S5D). DEGs of PI3K-Akt signaling pathway, such as those encoding the inflammatory cytokine IL-6 (IL-6), and vascular endothelial growth factors A and D (Vegfa, Vegfd) were increased upon Ang II stimulation (Supplementary Figure S5E). Overall, the DEG analysis with RNA-seq data indicates that the Ang II-induced cell senescence model is reliable and could be used for subsequent analysis.

Next, we used FLNC reads generated from Iso-Seq to perform isoform level analysis. FLNC reads were chosen because they can be considered as full-length isoforms (Wang et al., 2019). Although not as quantitative as RNA-seq, FLNC reads in Iso-Seq were still found to be semi-quantitative, especially for those full-length reads that derived from relatively highly expressed genes (Uapinyoying et al., 2020; Wyman et al., 2020). Therefore, we defined the isoform abundance as a percentage of isoform-supportive FLNC reads to the sum of all FLNC reads of a gene. The isoform with the largest percentage is called the major isoform, the isoform with abundance <0.1 is called inactive isoform, and the rest are called minor isoforms (Figure 2A). The abundance distribution of these three isoform types is in line with their definition (Figure 2B). The same trend was observed when quantification was performed with RNA-seq short reads (Figure 2C). We next compared the isoform abundance between Ang II and Control conditions. The isoform abundance between these two conditions was overall correlated (Figure 2D, Pearson’s correlation coefficient = 0.91; p < 2.2e-16). The abundance of most isoforms (78.62%, |Δ abundance| < 0.1) was not changed (Figure 2D). However, if we do not consider isoforms that are inactive before and after Ang II stimulation (n = 14,329) (Table 2), 36.2% of the isoforms undergo change (|Δ abundance| ≥ 0.1) in the proportion of usage before and after stimulation. Specifically, 3,548 isoforms showed upregulated expression and 3,470 isoforms showed downregulated expression upon Ang II treatment (Figure 2D). KEGG pathway enrichment analysis showed that genes with changed isoform usage were more likely to be enriched in senescence or aging-related pathways such as p53 signaling pathway, mTOR signaling pathway and MAPK signaling pathway when compared to those genes with unchanged isoform usage (Supplementary Figures S6A, B). We also used RNA-seq to examine the difference in isoform abundance and obtained similar results (Figure 2E; Table 3; Supplementary Figures S6C, D).

We selected Mxi1 and Cdh2 as examples to show the isoform usage difference upon Ang II treatment, both genes were associated with aging/senescence diseases (Schreiber-Agus et al., 1998; Mayosi et al., 2017; Kruse et al., 2019; Zhuo et al., 2019). Iso-Seq detected three isoforms of Mxi1, and two of these isoforms showed differential expression between Ang II and Control, despite no differential gene expression (Log2FC = 0.16), highlighting the importance of examining isoform level changes (Supplementary Figure S7). Cdh2 gene is widely studied in CVDs and other diseases (Mayosi et al., 2017; Kruse et al., 2019; Zhuo et al., 2019). Cdh2 showed PAS preference after Ang II stimulation, the major isoforms of Cdh2 in control tend to use the distal PAS while major isoforms in Ang II tend to use the proximal one (Supplementary Figure S8).

As alternative promoter (Demircioğlu et al., 2019; Huang et al., 2021), alternative splicing (Roy et al., 2013) and alternative polyadenylation (Chen et al., 2017) are three major contributors to isoform diversity in mammalian cells, we therefore respectively explored their contributions in Ang II-induced RAECs senescence. To explore the alternative promoter (AP) usage, we applied proActive (Demircioğlu et al., 2019) to carry out the promoter analysis, wherein the promoter usage was defined based on the first junction from transcriptome GTF file (here we used the GTF file generated from our Iso-Seq data), and the activity of the promoter was evaluated based on the number of junction reads supported by RNA-seq (Figure 3A). The promoter usage predicted by proActive is consistent with CAGE and H3K4me3 histone ChIP-seq data (Demircioğlu et al., 2019). We identified 1,109 genes with multiple promoters and 9,013 genes with single promoter (Figure 3B) using Iso-Seq data of RAECs. proActive classifies promoters into Major, Minor and Inactive according to their relative promoter activity. Promoters were also classified into known and novel promoters (Figure 3C, see Methods for details), and the known promoter accounted for the largest proportion, while the novel ones accounted for a small proportion (Figure 3D). We then applied DESeq2 to analyze differentially regulated promoters (DRPs) (Huang et al., 2021) between Control and Ang II-treated cells, and found 188 and 197 known promoters were upregulated and downregulated respectively in Ang II-treated RAECs. Notably, 18 and 10 novel promoters were upregulated and downregulated respectively upon Ang II stimulation (Figure 3E). Then we compared the genes with DRPs (DRPGs) with differentially expressed genes (DEGs) and found that 79.8% downregulated promoters overlapped with downregulated genes (Supplementary Figure S9A). These genes with downregulated promoters were enriched in functional categories such as cell division, nuclear division, negative regulation of cell cycle (Supplementary Figure S9B). 74.6% upregulated promoters overlapped with upregulated genes (Supplementary Figure S9C), and genes with upregulated promoters were enriched in Ang II stimulation-related terms such as cellular response to external stimulus, wound healing, regulation of vascular development, regulation of angiogenesis, and reactive oxygen species metabolic process (Supplementary Figure S9D). These results above indicated that the changed promoter expression may probably function in AngII-induced senescence in RAECs.

Next, we asked how many promoters are experiencing dynamic changes between Control and Ang II based on the promoter usage and relative promoter activity (Figure 3F). Relative promoter activity difference >0.1 and |FoldChange| > 1.5 were considered as alternative promoter (AP) between two conditions (Details in Methods section). A total of 239 APs derived from 121 genes were identified, including 34 genes with both upregulated and downregulated APs, which means switched promoter usage upon Ang II stimulation (Figure 3G). Three genes were selected as representatives for the promoter usage switch (Figure 3H). Map1b showed to have two promoters (Prmtr.29507 and Prmtr.29508) in Control, and after Ang II stimulation, Prmtr.29507 was downregulated (log2FC = −2.56) and Prmtr.29508 was upregulated (log2FC = 1.21) (Figure 3H left panel). In Cdk16, the inactive promoter (Prmtr.4399) in Control turned to be active after Ang II stimulation (log2FC = 1.95), though the major promoter (Prmtr.4398) kept not changed (log2FC = −0.21) (Figure 3H middle panel). As to Hpcal1, the activated promoter (Prmtr.11774) originally in Control turned to be inactive after Ang II treatment (log2FC = −1.84), but the major promoter (Prmtr.11773) did not show obvious change (log2FC = 0.35) (Figure 3H right panel). These observations suggest that APs serve as one contributor to the transcriptome diversity in Ang II-induced senescence in RAECs.

Alternative splicing (AS) is a well-known contributor to transcriptome diversity, and participates in diverse biological processes (Wright et al., 2022). Therefore, we examined whether alternative splicing also showed dynamic changes upon Ang II stimulation in RAECs using SUPPA2 (Trincado et al., 2018). AS patterns were generated based on the GTF file generated from Iso-Seq and quantified by Salmon quantification results (See Methods for details). All the AS events were classified into seven types: alternative 3′-acceptor (A3), alternative 5′-donor (A5), alternative first exon (AF), alternative last exon (AL), mutually exclusive exon (MX), intron retention (IR) and exon skipping (SE) (Figure 4A). Among the expressed 10,122 genes based on PacBio Iso-Seq, 7,321 genes expressed more than one isoform. The most frequent AS type was SE (6,968), followed by AF (4,618), A3 (3,701), A5 (3,026), AL (593), IR (455) and MX (356) (Figure 4B). Interestingly, the most frequent novel AS type found was AF (3,877), and IR (2,047) was the second (Figure 4B), suggesting these two types of AS events are especially less annotated in the ncbiRefSeq. Next, we compared the percent spliced-in (PSI) values between known and novel in the seven AS types. We found that the expression of known AS events are generally higher than the novel ones. In novel AS events, the expression of AF, AL, MX was relatively higher than the expression of other events (Figure 4C). SUPPA2 was next used for identification of differential AS events between Ang II and Control, which was evaluated with ΔPSI, the difference of the mean PSI between Ang II and Control. There were 666 significant different AS events between Ang II and Control, among of which, AF counts the largest proportion (n = 223, 33.5%) (Figure 4D). Of note, AF was the most frequent novel differentially spliced AS type and SE was the most frequent known differentially spliced AS type (Figure 4E).

We functionally categorized the differential AS genes above and found that they were enriched in Gene Ontology (GO) terms such as wound healing, regulation of intrinsic apoptotic signaling pathway, and regulation of cell growth (Figure 4F). The importance of differential AS in Ang II-induced RAEC senescence is indicated by the following two example genes. Castor1, a regulator of mTORC1 (Li et al., 2021), which plays a key role in cell proliferation in response to nutrition and growth stimuli, showed decreased intron retention level (a novel AS event for Castor1) upon Ang II stimulation (Figure 4G). Atp5mc1, a gene associated with CVDs (Hartmann et al., 2022), also displayed a usage shift of alternative first exon upon Ang II treatment (Figure 4H). These results indicate that AS can contribute to Ang II-induced RAECs senescence through isoform changes in responsible genes.

Alternative polyadenylation (APA) fine-tunes gene expression and largely enhances the complexity of mammalian transcriptome. To explore whether and to what extent APA contributes to the isoform changes in response to Ang II stimulation in RAECs, we performed APA analysis by integrating both Iso-Seq and RNA-seq data. We used Iso-Seq to identify Polyadenylation sites (PASs) and RNA-seq to quantify APA changes. PASs within 50 nt were clustered into one PAS, and PASs with distance >50 nt were considered as different PASs, as defined previously (Tang et al., 2022). If a PAS is within 50 nt of a known ncbiRefSeq PAS, it is considered as known, otherwise defined as novel (Figure 5A). We found that Iso-Seq detected a large number of PASs, including thousands of novel PASs (Figure 5B). To evaluate whether these novel PASs are authentic, we compared them with the experimentally verified 3′-ends pool generated by 3′READ, PolyA-Seq and polyA_DB3 of rat samples. We found that most PASs were within 50 nt to those in the validated 3′-ends pool for both known and novel PASs (Figure 5C). We analyzed the nucleotide composition of the 50 nt around the identified PASs and found that the 25 nt region upstream PAS is predominantly A-enriched, while the downstream is predominantly U-enriched (Supplementary Figure S10A), consistent with known literature reports (Chen et al., 1995). Moreover, we used MEME (Bailey et al., 2015) to discover motif and found the canonical polyadenylation signal (AAUAAA) also over-represented in the PASs (Supplementary Figure S10B) (Elkon et al., 2013). These above lines of evidence suggested that these PASs identified by Iso-Seq are of satisfied quality and reliable for further analysis. A total of 5,755 (57.32%) genes had at least two PASs based on Iso-Seq of RAECs (Figure 5D). We next analyzed the differential usage of PAS between Ang II and control RAECs using the modified algorithm of DaPars (Xia et al., 2014) for the differential PAS usage between Ang II and control RAECs. DaPars uses mean squared error (MSE) to predict the proximal PAS with mean RNA-seq coverage. Here Iso-Seq can provide reliable but not predicted PASs, we therefore used these PASs as prior knowledge and calculated their differential PAS usage between Ang II and Control cells by Fisher’s exact test (Figure 5E, see details in Method section). We found that 278 and 18 genes showed a preference usage of proximal and distal PASs, respectively, showing a global 3′UTR shortening upon Ang II stimulation (Figure 5F). Most of the length-changed 3′UTRs are derived from novel PAS (Figure 5F), suggesting Iso-Seq provides previously-unknown information for studying the APA regulation in Ang II-induced RAEC senescence. Further analysis showed that genes tending to use proximal PASs were enriched in pathways related to senescence such as p53 signaling pathway (Ciardiello and Tortora, 2008; Huang and Fu, 2015) (Figure 5G), implying that APA may fine-tunes certain key pathways to participate in regulation of Ang II-induced senescence.

As mentioned above, Iso-Seq and RNA-seq have their advantages and disadvantages. PacBio Iso-Seq data provides more reliable PASs even at relatively low coverage, and DaPars was originally designed to use Illunima’s RNA-seq data to predict proximal PASs usage and analyze global 3′UTR length changes between two conditions. Here we ran both the original DaPars program solely on RNA-seq data (Hereinafter referred to as RNA-seq) and the modified DaPars on both Iso-Seq and RNA-seq data (Hereinafter referred to as Iso-Seq), the results run with original DaPars used as comparison references. The DaPars result based on RNA-seq data also showed a global shortening trend of 3′UTR length after Ang II stimulation (Supplementary Figure S11A). DaPars using RNA-seq predicted a total of 260 differentially expressed proximal PASs, of which 119 were known ones (located within 50 nt to the reliable 3′-ends pool) (Supplementary Figure S11B). We found that most of the differential APA genes predicted by RNA-seq were also supported by Iso-Seq (Supplementary Figure S11C). In contrast, most of the differentially used APA events identified by Iso-Seq were not covered by RNA-seq, suggesting the necessity of integrating Iso-Seq and RNA-seq to achieve a more comprehensive APA results. Pten and Pnkp are two genes that were detected only in Iso-Seq but not in RNA-seq, their proximal PASs were supported by PolyA-Seq and 3′READ (Supplementary Figure S11D). Pten is a well-known tumor suppressor gene that functions in preventing cell survival and proliferation (Song et al., 2012). Our previous study discovered that SRSF3-mediated 3′UTR shortening of Pten could lead to senescence-associated phenotypes in both HUVEC and 293T human cells (Shen et al., 2019). Interestingly, this phenomenon is also evident in our Iso-Seq data in Ang II treated RAECs comparing to control ones (Supplementary Figure S12). Pnkp is a gene involved in DNA repair and in response to oxidative damage, it tends to use proximal PAS and therefore shortened 3′UTR after Ang II stimulation in RAECs (Supplementary Figure S11D). Mpc1 was only detected to display 3′UTR shortening only in RNA-seq but not in Iso-Seq. Both the proximal and distal PASs of Mpc1 were detected by PolyA-Seq and 3′READ, but Iso-Seq could only detect the distal one, possibly due to its limited sequencing coverage (Supplementary Figure S11D). These results demonstrate the power of integrating Iso-Seq and RNA-seq in comprehensive understanding of APA regulation in Ang II-induced RAEC senescence.