Pre-eclampsia intronic polyadenylation enriched in VEGFA-VEGFR2 signaling pathway.

Junhua Zhang: Writing – review & editing, Writing – original draft, Conceptualization. Yingying Lu: Writing – review & editing, Writing – original draft, Resources, Conceptualization. Lei Li: Writing – review & editing, Writing – original draft, Data curation. Xiongying Li: Writing – review & editing, Writing – original draft, Formal analysis. Jingxia Ying: Writing – review & editing, Writing – original draft, Formal analysis. Sicong Li: Writing – review & editing, Writing – original draft, Formal analysis. Lingling Wu: Writing – review & editing, Writing – original draft, Formal analysis. Lijing Li: Writing – review & editing, Writing – original draft, Supervision.

This work was supported by a grant from by Jinhua Key program (No. 2022-3-051 ).

We collected transcriptome samples from women with pre-eclampsia to investigate intronic polyadenylation (IPA) selection using IPAFinder. Additionally, we obtained transcriptome data from normal pregnancy samples and documented the sample information in Table S1 . All data used in this study were mRNA samples that were enriched for polyadenylated RNA. We analyzed the correlation between IPA events and gene regulation of mRNA and miRNA in pre-eclampsia. Moreover, we identified miRNAs that were specific for IPA events present in the transcriptome of women with pre-eclampsia. We predicted IPA event regulation and conducted validation in the circulating RNA dataset of pre-eclampsia patients ( Fig. 1 ). Fig. 1 Flow-chart. IPA, intronic Polyadenylation. CfRNA, Cell free RNA. Fig. 1 Flow-chart. IPA, intronic Polyadenylation. CfRNA, Cell free RNA. In PE-related IPA events, we observed only a small number of up-regulation events (genes RPS29 etc.), while more IPA down-regulation events were observed ( Fig. 2 A). Additionally, cluster analysis revealed that IPA events, except for two samples, were clustered into distinct sub-specific IPA down-regulated classes. The IPA analysis showed that YPEL5 was the most significantly down-regulated gene (p-value = 1.76e-117), and the gene with the largest fold change was RIPOR2 (FC = 0.75) ( Fig. 2 B). Furthermore, our findings demonstrate that the VEGFA-VEGFR2 signaling pathway (WP3888) was functionally enriched with PE-IPA ( Fig. 2 C). Other significantly enriched pathways include the trbp containing complex, actin filament-based process (GO: 0030029). Gene network module analysis revealed that the VEGFA-VEGFR2 signaling pathway ( Table S2 ). And the regulation of mRNA metabolic process (GO: 1903312) were significantly clustered networks ( Fig. 2 D). Compared to only four genes in the regulation of mRNA metabolic process (GO: 1903312), the VEGFA-VEGFR2 signaling pathway had ACTG1 as the central node, and TPM3 , TLN1 , TPM4 , MYL6 , ARPC5 , WIPF1 as the network core genes. Collectively, our results establish that IPA changes in the VEGFA-VEGFR2 signaling pathway greatly influences transcriptome changes in pre-eclampsia. Fig. 2 Identification of pre-eclampsia related IPA event. (A) Heatmap of IPA score in pre-eclampsia transcriptome compared to normal samples. (B) Volcano plot of log2-fold change and p-value from IPA-finder. (C) Functional enrichment of IPA related pathways. (D) Network analysis of IPA genes. Fig. 2 Identification of pre-eclampsia related IPA event. (A) Heatmap of IPA score in pre-eclampsia transcriptome compared to normal samples. (B) Volcano plot of log2-fold change and p-value from IPA-finder. (C) Functional enrichment of IPA related pathways. (D) Network analysis of IPA genes. To determine the effect of IPA events on transcriptional regulation, we analyzed the differential gene expression of pre-eclampsia and normal samples. Our analysis identified 2540 downregulated genes and 3690 upregulated genes (| log (FC) | > 1, p-adjusted <0.05). Functional annotation of these differential genes showed that the down-regulated DEG genes in pre-eclampsia samples were primarily enriched in female pregnancy pathways (GO: 0007565), as well as development-related pathways ( Fig. 3 A), such as GO: 0060485 (mesenchyme development). Conversely, the up-regulated DEG genes were mainly enriched in humoral immune response and cell membrane ion potential pathways ( Fig. 3 B). The pathway enrichment results of the VEGFA-VEGFR2 signaling pathway did not show significant results (p-adjusted = 0.052). Therefore, we checked the correlation between gene expression in the VEGFA-VEGFR2 signaling pathway and IPA events to determine whether this result was affected by too many differential genes and caused overly harsh fold change cutoff. By examining the VEGFA-VEGFR2 signaling pathway genes listed in the DEG table with adjusted p-value less than 0.05, We found VEGFA-VEGFR2 signaling pathway genes expression are up-regulation in the PE samples (log(DE_FC) > 0) ( Fig. 3 C).A significant negative correlation between the degree of IPA (log(IPA_FC)) and the degree of differential expression (log(DE_FC)) of this pathway (Spearman r = 0.27, p-value = 0.002). Considering distribution of log (IPA_FC) in VEGFA-VEGFR2 signaling pathway DEGs are closely around default IPAfinder cutoff ( Fig. 3 C), We relaxed the restriction of the cutoff to |log (IPA_FC) | > 0.01 of IPA FC. We observed that the pattern of IPA down-regulation remained maintained in the VEGFA-VEGFR2 signaling pathway DEGs ( Fig. 3 D). The core gene ACTG1 had the most significant p-value (p-value < 1e-50) in this pathway, indicating that the IPA of this gene may have a significant change effect. Fig. 3 Gene expression change and GO enrichment of down-regulated (A) and up-regulated (B) genes. (C) Correlation of IPA fold change to gene expression change levels in VEGFA-VEGFR2 pathway. (D) Volcano plot of IPA events in VEGFA-VEGFR2 signaling pathway. Fig. 3 Gene expression change and GO enrichment of down-regulated (A) and up-regulated (B) genes. (C) Correlation of IPA fold change to gene expression change levels in VEGFA-VEGFR2 pathway. (D) Volcano plot of IPA events in VEGFA-VEGFR2 signaling pathway. We continued to confirm the presence of IPA events in circulating RNA by examining their occurrences in blood samples. We observed only a slight signal indicating that IPA events of ACTG1 were expressed in the circulating blood samples ( Fig. S1 ). However, due to the lower sequencing quality and shorter sequencing length of circulating RNA, the signal of IPA events in blood may be underestimated. To investigate the potential functional consequences of IPA events in pre-eclampsia, we analyzed genome sequences in the 50bp upstream and downstream of IPA changes sites and identified RNA-binding proteins (RBPs) that were likely to bind ( Table 1 ). Our results confirmed with the IPA result and revealed an enrichment of poly-A sequences near the IPA events ( Fig. 4 ). The top RBPs that were predicted to bind to these sequences were LIN28B (Lin-28 Homolog B) and AGO2 (Argonaute RISC Catalytic Component 2), with AGO2 having the highest coverage (coverage = 0.75). The most significant motifs of AGO2 were the PAZ and Piwi motifs (p-value = 7.00e-14). Considering the involvement of AGO2 in the IPA-enriched pathway trbp containing complex in miRNA regulation, our findings suggest a potential link between miRNA binding to IPA related 3′ untranslated region and the pathogenesis of pre-eclampsia. Table 1 Enriched motifs of IPA sites. Table 1 Protein Domains p-value Coverage LIN28B CSDx1; CSDx1; Znf_CCHCx2; Znf_CCHCx2; 5.10e-15 0.52 AGO2 PAZ; Piwi; 7.00e-14 0.75 ALKBH5 Alpha-ketoglutarate binding; 3.50e-12 0.56 HuRMNase RRMx1; RRMx2; RRMx3; 2.90e-11 0.59 RBM10 G-patch; RRMx1; RRMx2; 2.40e-10 0.54 EWSR1 RRMx1; Znf_RanBP2; 4.80e-10 0.53 Fig. 4 Enriched motif of IPA-related sequence of LIN28B (left) and AGO2 (right). Fig. 4 Enriched motifs of IPA sites. Enriched motif of IPA-related sequence of LIN28B (left) and AGO2 (right). The analysis of sequences near IPA events and miRNA binding sites identified four miRNAs with distinct network cores in IPA-related miRNA targeting networks. Of these, miR-193b and miR-365a had the most potential binding sites ( Fig. 5 ). Interestingly, miRNA targeting ACTG1 (Actin Gamma 1) was not found among the top four core miRNAs. In the VEGFA-VEGFR2 signaling pathway, SH3BGRL3 (SH3 Domain Binding Glutamate Rich Protein Like 3) was regulated by hsa-miRNA-33b-3p, and the IPA event of SH3BGRL3 was down-regulated with a significant increase in expression, possibly due to the decrease in hsa-miRNA-33b-3p bound to it. Other examples of miRNA targeting in VEGFA-VEGFR2 signaling pathway included SIAH2 (Siah E3 Ubiquitin Protein Ligase 2, regulated by miR-365a-5p and miR-193b-5p), MYH9 (Myosin Heavy Chain 9, regulated by miR-193b), and RPLP2 (Ribosomal Protein Lateral Stalk Subunit P2, regulated by miR-33b-3p), which suggest that IPA events are accompanied by abnormalities in miRNA targeting regulation. Fig. 5 Network of miRNA target in IPA related UTR regions. Fig. 5 Network of miRNA target in IPA related UTR regions.

To find representative PE transcriptome data, we searched Gene Expression Omnibus (GEO) database with keyword “Pre-eclampsia”. Dataset with sample less than six sample or missing disease annotation were filtered. In total messenger RNA (mRNA) sequencing data of 90 samples from dataset GSE148241 [ 17 ] and GSE192902 [ 18 ] were extracted from the GEO ( https://www.ncbi.nlm.nih.gov/geo/ ) using SRAtoolkit v2.9.2. The patient information was retrieved from the Sequence Read Archive (SRA) metadata. We collected metadata table of sample in “severe pre-eclampsia” samples in “disease” column were choose as PE samples in disease column (N = 57), “Normal” samples noted in “subject_status” column were choose as normal samples (N = 33). Sample with empty disease note were excluded from following analysis. All FASTQ files extracted from SRA were pre-processing using fastp v0.20.1 [ 19 ] and aligned to the human genome (hg38) with STAR v2.7.7a [ 20 ]. IPA events were analyzed using IPAfinder software (version20230731) from alignment files using IPAFinder_DetectIPA.py with default parameter [ 13 ]. Differential usage of IPA terminal exon between PE and normal groups was inferred with Infer_DUIPA module with a cutoff of abs(IPUI_diff) > 0.1. Alignment BAM files of samples listed in the IPUI_diff results were loaded into Integrative Genomics Viewer (IGV). Genome positions with the highest IPA scores were visually inspected for read coverage supporting IPA events. Gene expression count tables were estimated utilizing featureCount [ 21 ] with hg38.ensGene GTF downloaded from UCSC goldenPath website: https://genome.ucsc.edu . Differentially expressed genes (DEG) were evaluated using edgeR [ 22 ] with cutoff of |log2 (fold change) | > 1 and adjusted p-value <0.05. Function enrichment was performed through Clusterprofiler [ 23 ] with cutoff of false discovery rate (FDR) < 0.05. Miranda software (version 3.3a) was used to perform miRNA target gene prediction [ 24 ]. The regulatory networks were visualized using Cytoscape software with Kappa similarities above 0.3 [ 25 ]. Additionally, the BRIO web server MOTIF Prediction tool [ 26 ] was used for RNA sequence and structure motif scan to identify RNA-binding proteins enrichment around IPA event in our analysis. All statistical analyses were performed using R version 4.0.3. The spearman correlation between IPA fold change (FC) and DEG fold change (FC) was estimated using the lm model, and significance was determined using the function cor.test. Statistical significance was set at p-value <0.05. Function enrichment p-values were adjusted using false discovery rate (FDR).

Pre-eclampsia (PE) is a multifaceted pregnancy disorder and research on Pre-eclampsia biomarkers remains a topic of great interest. Recent studies have suggested that transcriptome-level abnormalities may play a crucial role in the pathogenesis of PE [ 6 ]. In this study, we analyzed Intronic polyadenylation (IPA) events in transcriptome data from PE samples and found significant differences in IPA events compared to normal samples. In addition, we observed a negative correlation between gene expression levels and IPA events in VEGFA-VEGFR2 pathway. Our analysis of miRNA binding sites associated with these IPA events revealed the potential regulatory roles of miR-193b and miR-365a in IPA genes. Previous studies have reported alterations in transcriptome polyadenylation sites in pre-eclampsia, with changes in the FLT1 gene being particularly noteworthy [ 27 ]. In our study, we found a negative correlation between IPA events and gene expression for the VEGFA-VEGFR2 signaling pathway. Notably, IPA of ACTG1 showed the most significant p-value, despite only rare variations being detected. On the other hand, SH3BGRL3 showed maximum fold change in IPA throughout the study and the most significant differential gene expression. Consistent with our findings, SH3BGRL3 was reported significant serum protein changes in PE [ 28 ]. Our IPA study has potential implications for further research in this area. Our motif enrichment analysis suggested that IPA events may be associated with miRNA binding. Specifically, our miRNA binding network analysis revealed that miR-193b and miR-365a play central roles in the regulation of IPA events. Previous research has validated the function of miR-193b in PE [ 29 ]. Within the VEGFA-VEGFR2 signaling pathway, the SIAH2 gene IPA is significantly regulated by miR-193b and miR-365a. Moreover, miR-365a-5p has been found to be significantly enriched in exosomes from tumors and osteonecrosis diseases [ 30 ]. SH3BGRL3 is also regulated by miR-33b-3p. Up-regulation of miR-33b has been shown to promote endometriosis through Zinc finger E-box binding homeobox 1 ( ZEB1 ) expression via the Wnt/β-catenin signaling pathway [ 31 ]. These findings suggest potential clinically relevant connections between IPA-associated miRNA and gene expression change in pre-eclampsia. Our future work will focus on developing an extracellular vesicle-based therapy [ 32 ] utilizing the IPA-related miRNAs identified in this study. It is important to note that our study had a small sample size, and the current IPA events reported in the literatures were primarily analyzed in cell lines due to the difficulty of obtaining clinical samples [ 11 , 33 ]. We are conducting an ongoing collection of clinical blood sample RNA-sequencing in our hospital [Zhang, manuscript under preparation]. As we continue our research, we aim to expand the available datasets to enhance the reliability and broad applicability of our findings.

Pre-eclampsia (PE) is a significant complication of pregnancy that leads to maternal and perinatal morbidity and mortality [ 1 ]. While the exact cause of pre-eclampsia remains unclear, it has been established that this condition results from the elevated presence of antiangiogenic proteins in the mother's blood, particularly the soluble fms-like tyrosine kinase 1 (sFlt1s) protein [ 2 ]. The excessive amount of sFlt1 affects the mother's ability to respond to vascular endothelial growth factor ( VEGF ), which is crucial for maintaining healthy blood vessels during pregnancy [ 3 ]. One possible treatment involves using short interfering RNA (siRNA) to target specific sFlt1 mRNA sequences to effectively silence sFlt1 expression [ 4 ]. However, possible therapies warrant further examination of whether other VEGF-related proteins might be structurally truncated or regulated by other small RNA [ 5 , 6 ]. Polyadenylation is a post-transcriptional modification that involves the addition of poly (A) tails to the 3′ end of mRNA. This process is essential for maintaining mRNA stability and enabling efficient translation regulated by microRNA [ 7 , 8 ]. The mechanisms responsible for defining poly (A) sites include poly (A) signals that are typically found in pre-mRNAs and consist of AAUAAA hexamers.Recent alternative polyadenylation landscape has revealed that nearly half of all human genes feature multiple 3′ UTRs that allow for tissue-specific post-transcriptional regulation in universally expressed genes [ 9 ]. However, both fms-like tyrosine kinase 1 ( FLT1 ) analysis and transcript-wide polyadenylation site analysis in pre-eclampsia failed to find any significant changes in alternative polyadenylation associated with PE in past decades [ 10 ]. Recent studies indicate polyadenylation events can occur in introns instead of 3′ UTRs, resulting in transcripts with truncated coding regions and C-terminal domains loss [ 11 ]. Genes responsible for producing intronic polyadenylation (IPA) isoforms are found to have competitive splicing and polyadenylation signals with diverse functional impact [ 12 ]. For instance, recently matured bioinformatics tools have emerged that can identify the widespread expression of IPA isoforms in messenger RNA sequencing data [ 13 ]. Similar tools have been used in examination of transcriptome IPA in tumor cells, immune cells and other contexts [ 14 , 15 ]. This process allows us to thoroughly examine IPA-related changes in genes of the VEGF receptor pathway. Understanding these changes is crucial for comprehending the basics of the PE transcriptome and exploring potential therapeutic developments surrounding this molecular mechanism [ 16 ]. This study aims to 1) utilize the latest transcriptome data from public databases to investigate IPA events in patients with pre-eclampsia. 2) search for potential regulation of miRNA on IPA events.

Argonaute RISC Catalytic Component 2 Differentially expressed genes Fold change Intronic polyadenylation Lin-28 Homolog B micro RNA Pre-eclampsia untranslated regions vascular endothelial growth factor

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Data will be made available on request.