Alternative Polyadenylation Results in mRNA Transcript Instability in Gestational Diabetes Mellitus

Introduction

Gestational diabetes mellitus (GDM) is defined as glucose intolerance of various degrees that is first detected during pregnancy.¹ GDM is a common chronic disease during pregnancy, damaging the health of millions of women around the world.^2,3 Nowadays, the incidence rate of gestational diabetes is increasing, and it has reached 9–25% worldwide due to different ethnic groups, regions and other elements.^4,5 The obvious feature of GDM is impaired glucose tolerance as a result of maternal pancreatic β-cell dysfunction, which leads to imbalance of glucose and lipid metabolism during pregnancy.⁶ GDM has short-term and long-term adverse effects on health of mother and fetus. However, until now, the pathophysiology of GDM is not fully understood.

In eukaryotes, it is estimated that 30%–70% of the genes will produce mRNA with multiple selectable 3’UTRs, which is called selective polyadenylation (APA). The APA modification of mRNA has become an important mechanism of post-transcriptional gene regulation in higher eukaryotes. The 3’UTR regions with variable length may lead to different levels of protein expression. In order to better explore the mechanism of APA’s role in GDM, we chose a method based on deep sequencing, called Poly(A) site sequencing (PAS-seq), to quantitatively analyze RNA polyadenylation at the transcriptome level.⁷ When the RNA chain extends to the transcription termination site (TTS), the process of transcription is terminated. The change of TTS site leads to the change of 3’UTR length, and further changes the interaction with microRNAs or other regulatory RNA binding proteins (RBPs) to affect the stability, location, transport and even translation of mRNA. The target site of miRNA exists in the 3’UTR region. The change of the length of the 3’UTR region during APA affects the function of miRNA. RBPs have similar effects. They bind to the 3’UTR region and exist as mRNA unstable elements.⁸

Nowadays, in the exploration of GDM, numerous sequencing research have been published, which provide new ideas and directions for the pathogenesis and clinical potential application of GDM. For the first time, we introduced the change of APA events in GDM, and pointed out some differential genes that may play a key role in it. In conclusion, we believe that exploring the mechanism of APA events participating in the occurrence and development of GDM is of positive significance for the study of GDM, including its risk factors, biomarkers and health management.

Materials and Methods

Sample Collection

GDM pregnant women and healthy pregnant women who were delivered by Caesarean section at Shengjing Hospital of China Medical University from January 2022 to June 2022 were selected as the study subjects. Each diagnosis of GDM was based on the Guidelines for the Prevention and Treatment of Type 2 diabetes in China (2020). A 75-g oral glucose tolerance test (OGTT) was performed as required between 24 and 28 weeks during pregnancy. A diagnosis of GDM was made if the blood glucose of the patient met the following criteria at any point: 5.1 mmol/L≤fasting blood glucose <7.0 mmol/L, OGTT 1h blood glucose ≥10.0 mmol/L, or ≤8.5 mmol/L OGTT 2h blood glucose <11.1 mmol/L. The control group comprised healthy pregnant women who were hospitalized for caesarean section during the same period as those in the GDM group. The pregnant women in two groups all had single pregnancies without pregnancy hypertension, heart disease, infectious disease, serious liver and kidney disease, or endocrine disease, or other pregnancy complications. Ultimately, three normal and three GDM women were selected as the study subjects. The general clinical data of the subjects is shown in Table 1.

Table 1 Clinical Information of the Subjects

Within 10 minutes of the placenta being delivered, tissue was collected (about 1cm³ from the maternal side of the placenta within 2cm of the diameter of the umbilical cord root) and was immediately placed into liquid nitrogen for preservation after washing with normal saline. The study was conducted in accordance with the Declaration of Helsinki (as revised in 2013). This research scheme was approved by the Ethics Committee of Shengjing Hospital of China Medical University [2022PS313K].

RNA Extraction and Sequencing

Total RNA of tissue was extracted by the TRIZOL (Ambion). The RNA was further purified with two phenol-chloroform treatments and then treated with RQ1 DNase (Promega, Madison, WI, USA) to remove DNA. The quality and quantity of the purified RNA were redetermined by measuring the absorbance at 260 nm/280 nm (A260/A280) using Nanodrop one (Thermo). RNA integrity was further verified by 1.5% agarose gel electrophoresis. For each sample, 5 μg of total RNA was treated with RQ1 DNase (Promega) to remove DNA, fragmented RNAs were used for RNA-seq library preparation. mRNAs were captured with mRNA Capture Beads kit (Vazyme, N401-01). Fragmented mRNAs were used for directional PAS-seq library preparation by KAPA Stranded mRNA-Seq Kit for Illumina^® Platforms (KK8544). Fragmented mRNAs were transcribed into cDNA with pas-RT primer; following end repair and A tailing, the DNAs were ligated to Diluted Roche Adaptor (KK8726). After purification of ligation product and size fractioning to 300–500 bp, the ligated products were amplified and purified, quantified and stored at -80°C before sequencing. The strand marked with dUTP (the second cDNA strand) is not amplified, allowing strand-specific sequencing. For high-throughput sequencing, the libraries were applied to Illumina Novaseq 6000 system for 150 nt paired-end sequencing.

RNA-Seq Raw Data Clean and Alignment

Raw reads containing more than 2-N bases were first discarded. Then, adaptors and low-quality bases were trimmed from raw sequencing reads using FASTX-Toolkit (Version 0.0.13). The short reads less than 16nt were also dropped. After that, clean reads were aligned to the GRCh38 genome by HISAT2 allowing four mismatches. Uniquely mapped reads were used for gene reads number counting and FPKM calculation (fragments per kilobase of transcript per million fragments mapped).

Peak Calling

After mapping, pA site information was extracted. In pas-seq data analysis, the localization of the 3’ end of effective poly(A) reads that can be compared to the reference genome was defined as poly(A) site. CAGEr was used for screening and statistical analysis of these poly(A) sites. Values of shifting score are in range between -Inf and 1. Positive values can be interpreted as the proportion of transcription initiation in the sample with lower expression that is happening “outside” (either upstream or downstream) of the region used for pA sites in the other sample. In contrast, negative values indicate no physical separation. The shift distance is at least greater than 10nt.

Key Genes Analysis

The R Bioconductor package DESeq2 was utilized to screen out the differentially expressed genes (DEGs). The P value <0.05 and fold change >1.5 or <0.67 were set as the cut-off criteria for identifying DEGs. Target scan (8.0), a miRNA prediction database, was used to predict miRNAs associated with key genes. RBPs of key gene interactions were predicted by starbase database, and DEGs in RNA-seq were analyzed for overlap.

Results

Preliminary Analysis of PAS-Seq Sequencing Data

We conducted high-throughput sequencing and database building for the samples. The database was built using illumina Nextseq 500 sequencing platform for double-ended sequencing, and high-quality transcriptome data were obtained. We take the peak of poly(A) site as the center, extend 40bp upstream, extract the sequence within this range, and then conduct sequence feature enrichment analysis. It is found that the most enriched motif near poly(A) site in the sample is AAUAAA, which is consistent with the reported classical poly(A) signal sequence, as shown in Figure 1.

Figure 1 Sequence analysis of poly(A) loci.

The clinical information of the sample is shown in Table 1. Relevant results of RNA-seq are shown in Figure S1.

PAS-Seq Reveals Alternative Polyadenylation Events in GDM

Through the analysis of PAS-seq database, 202 TTS loci were found to have different changes, including 112 genes with TTS locus in advance and 109 genes with TTS locus in delay. There were 57 genes with significant difference in TTS (P<0.05). The genes in advance of TTS locus were enriched in biological processes such as cell development, protein transport and phosphorylation, signal transduction, etc. The delayed genes at TTS site are enriched in such biological processes as transcriptional regulation, cell migration and cycle, DNA repair and damage. These results are shown in Figure 2. By searching the documents related to GDM and TTS locus genes with significant differential changes, 13 authentic and credible genes, including ADHFE1, CD44, CHURC1, COLEC11, KDELR2, LGR4, LIPT2, VN1R1, ZNF701, DDX58, VCPIP1, LYPD3 and GIMAP5, were found worthy of attention.

Figure 2 Poly(A) site sequencing revealed alternative polyadenylation events in gestational diabetes mellitus. (A) PCA base on FPKM value of all detected genes in PAS-seq. The ellipse for each group is the confidence ellipse; (B) scatter plot showing the number of genes changing TTS. Gray dots refer to genes with no statistical significance; (C) enrichment of around transcriptional start site (TSS) and terminal site (TTS). Normalized reads count (RPM) in a ± 1000 bp window around TSS and TTS for all transcripts. RPM, reads per million; (D) the Scatter plot exhibiting the most enriched GO biological process results of the difference with TTS.

VCPIP1 and LGR4 Change TTS Site and mRNA Expression in GDM

We found that a total of eight genes changed simultaneously in PAS-seq and RNA-seq. Among them, the mRNA of VCPIP1 and LGR4 had significant differential expression. In the sequencing, the mRNA expression of LGR4 and VCPIP1 was decreased in the GDM group. The TTS locus of LGR4 was shortened in advance, that is, 3’UTR, while the TTS locus of VCPIP1 was delayed, that is, 3’UTR was prolonged, as the results can be seen in Figure 3.

Figure 3 Alternative polyadenylation in the mRNAs results in transcript instability. (A) The TTS distribution of LGR4, Left panel: IGV-sashimi plot showing the distribution of TTS locus reads, the green panels represent the position of changing. Reads distribution is plotted in the up panel and the transcripts of each gene are shown below. Right panel: The Bar plot showing the expression pattern and statistical difference of DEGs. Error bars represent mean±SEM. (B) The TTS distribution of VCPIP1, Left panel: IGV-sashimi plot showing the distribution of TTS locus reads, the green panels represent the position of changing. Reads distribution is plotted in the up panel and the transcripts of each gene are shown below. Right panel: The Bar plot showing the expression pattern and statistical difference of DEGs. Error bars represent mean±SEM.

Prediction miRNAs and RBPs Combined to Key Genes

The miRNA and RBP of VCPIP1 and LGR4 were predicted. It was found that VCPIP1 bound has-miR-499a-5p, while LGR4 predicted two miRNAs, hsa-miR-520f-3p and hsa-miR-302c-3p. Through literature search, it is found that has-miR-499a-5p has a high correlation with GDM. Then, look for other targets of miR-499a-5p, and conduct overlapping analysis in DEGs. It is found that 24 target mRNA have significant differential expression in RNA-seq of this project, of which 2 genes are up-regulated and 22 genes are down-regulated (Figure 4).

Figure 4 miR-499a-5p binds to the VCPIP1 3’UTR. VCPIP1, valosin-containing protein interacting protein 1. (A) The network plot showing the miR-499a-5p and its target mRNA; (B) Venn diagram showing the overlapped genes between DEGs and target mRNA; (C) Hierarchical clustering heat map showing expression correlations of overlap genes with VCPIP1.

The starbase database was used to predict the RBPs interacting with LGR4 and VCPIP1, and the overlap analysis with DEGs showed that FMR1 and DICER1 were the most reliable, and they were differentially expressed in RNA-seq that we can see in Figure 5.

Figure 5 RBP binds to VCPIP1, and LGR4 influences mRNA stability. RBP, RNA binding protein; VCPIP1, valosin-containing protein interacting protein 1; LGR4, leucine-rich repeat-containing G-protein coupled receptor 4. (A) Venn diagram showing the overlapped genes between DEGs and target RBP of LGR4 and target RBP of VCPIP1; (B) Bar plot showing the expression pattern and statistical difference of DEGs for overlap genes. Error bars represent mean±SEM.

Discussion

GDM refers to diabetes with normal glucose metabolism before pregnancy or potential impaired glucose tolerance, which occurs or is diagnosed during pregnancy. More than 80% of abnormal glucose metabolism in pregnant women is GDM. The glucose metabolism of GDM patients can return to normal after delivery, but the risk of T2DM will increase in the future.^6,9 In the second and third trimesters of pregnancy, the increase of anti-insulin-like substances in pregnant women, such as placental lactogen, estrogen, progesterone, cortisol and placental insulinase, makes the sensitivity of pregnant women to insulin decrease with the increase of gestational weeks. In order to maintain normal glucose metabolism, insulin demand must be increased accordingly.¹⁰

GDM has short-term and long-term adverse effects on health. For mothers, the short-term will increase the probability of adverse pregnancy events, and the long-term will increase women’s risk of T2DM and cardiovascular disease (CVD).^11,12 For the fetus, the probability of macrosomia is significantly increased, and infant hypoglycemia occurs.^13–16 If not handled properly, it will also cause brain damage.¹⁷ In the long term, infants born with GDM pregnancy are at increased risk of obesity, T2DM, CVD and related metabolic diseases.⁶

When the RNA chain extends to the TTS site, the RNA polymerase will no longer form a new phosphodiester bond, the RNA DNA hybrid will be separated, the transcription bubble will collapse, and the DNA will return to the double-stranded state, while the RNA polymerase and RNA chain will be released from the template, which is the termination of transcription. The change of TTS site makes the transcript 3’UTR lengthen or shorten. The median length of 3’UTR varies from 140nt to 2300nt. RBPs bound by 3’UTR interact with cell movement or other proteins, which may guide variable 3’UTR to be located in different regions to perform different functions. 3’UTR binding to AU-rich RBPs or miRNAs will lead to instability of mRNA and shorten the half-life of mRNA. 3’UTR binding RBPs interact with translation initiation factors to regulate mRNA translation. In addition, it has also been reported that the function of 3’UTR not only affects the relevant attributes of mRNA. It can also affect newborn proteins, which recruit proteins to translation sites by acting as scaffolds, so that proteins can form complexes with newborn peptide chains, and regulate the formation and positioning of membrane proteins. It has also been reported that RBP can regulate the variable 3’UTR events of pre-mRNA.¹⁸

Through our analysis of sequencing data, we found that there are two differentially expressed genes, VCPIP1 and LGR4, which are most likely to cause changes in mRNA due to changes in TTS sites.

VCPIP1, also known as VCIP135, is a kind of ubiquitination enzyme, which is related to the reorganization of Golgi apparatus and endoplasmic reticulum after mitosis, and is also involved in DNA protein cross-linking (DPC) repair. It is reported that VCPIP1 can also stably catalyze light-chain metalloproteinases. Other studies have shown that VCPIP1 −/− mice are born normally, and they are prone to genomic instability and premature aging. VCPIP1 knockout also causes hypersensitivity to DPC inducers, Top1cc accumulation, chromosomal abnormalities and micronucleus formation, indicating that VCPIP1 is crucial in DPC repair and maintaining genomic stability.¹⁹ It has been found that in GDM, the expression of Has-miRNA-499a-5p is increased, and it may be related to hypertension in pregnancy, preeclampsia, etc.²⁰ These studies suggest that VCPIP1 in GDM may lead to the prolongation of 3’UTR region through the change of TTS site, leading to the increase of binding miRNA-499a-5p, thus affecting the stability of VCPIP1 mRNA, further affecting its protein expression, and reducing its expression.

LGR4 is a G protein-coupled receptor-4 rich in leucine sequence, belonging to the B subfamily of G protein-coupled receptor superfamily, which controls multiple pathways and biological processes. At the molecular level, LGR4 is involved in different intracellular signal transduction pathways, including Wnt, NFκB, TGFβ and cAMP, which indicates that LGR4 is involved in excessive mechanisms that control cell events in a tissue-specific manner. LGR4 also plays an important role in metabolic regulation, coordinating various metabolic signals.²¹ Recently, Dong et al^22,23 described RSPO2/LGR4 as a key regulator for the differentiation of precursor cells into adipocytes. In addition, Li et al observed that the level of LGR4 in diabetes patients was significantly lower than that in non-diabetes patients (P<0.05), and the level of LGR4 in hypertensive patients was lower than that in non-hypertensive patients with higher blood pressure.

To sum up, this study screened five signal pathways and eight key genes through high-throughput sequencing of GDM placental tissue and bioinformatics analysis of data. Through literature search, we further found two genes that may play a key role in the occurrence and development of GDM, VCPIP1 and LGR4, which are expected to become the susceptible genes predicted by GDM, new molecular markers for diagnosis, and new targets for the development of targeted drugs for treatment.

Conclusion

The abnormality of APA process may be involved in the occurrence and development of GDM. The genes with significantly different changes in TTS locus may become biomarkers or predictors for GDM to assess the incidence, disease progression and disease severity, and may also become potential targets for GDM treatment.

We speculate that LGR4 and VCPIP1 may be related to TTS site changes in the placenta of pregnant women with diabetes, and miRNAs or RBPs bind to 3’UTR of LGR4 and VCPIP1 to inhibit their expression that to affect the stability of its mRNA, and ultimately down regulate its expression.