The small interfering RNA (siRNA) duplexes for LGALS1 transcripts were purchased from Gemma (Suzhou, China). Non-targeting control siRNA (si-Negative) sequence was: 5’-UUCUCCGAACGUGUCACGUTT-3’ (sense), while the siRNA targeting LGALS1 (siLGALS1-3) sequence was: 5’-GCAAAGACAGCAACAACCUTT-3’ (sense). Three biological replicates were prepared for the subsequent experiments in this study.

HUVECs were originally obtained from ScienCell (8000, ScienCell, Carlsbad, CA, USA), and provided by Shanghai Zhong Qiao Xin Zhou Biotechnology Co.,Ltd. (Shanghai, China). The cells were authenticated by immunofluorescence analysis of the positive marker gene of vWF, and tested for the free of mycoplasma contamination through PCR.

HUVECs were cultured at 37 °C with 5% CO${}_2$ in endothelial cell medium (ECM, ScienCell) containing 10% fetal bovine serum (FBS), 100 $\mu$g/mL streptomycin, and 100 U/mL penicillin. The siRNA transfection targeting LGALS1 was performed using Lipofectamine 2000 (Invitrogen, Carlsbad, CA, USA) according to the manufacturer’s protocol. After 48 h, the RT-qPCR and western blot (WB) experiments were conducted to validate the repressed expression of LGALS1. Transfected HUVECs were harvested using a previously published procedure and the 2${}^{-\mathrm{\Delta }\mathrm{\Delta }\text{CT}}$ method was employed [19]. The primers for RT-qPCR are shown in Supplementary Table 1. For WB experiments, samples were boiled for 10 min in 1$\times$ sodium dodecyl sulfate (SDS) sample buffer and separated by 10% SDS-PAGE. They were then incubated with LGALS1 primary antibody (1:1000 dilution, A18040, ABclonal, Wuhan, China) or GAPDH primary antibody (1:2000, 60004-1-IG, Proteintech, Wuhan, China), followed by HRP-conjugated secondary antibody. Bound secondary antibody (anti-mouse or anti-rabbit 1:10,000) (Abcam, London, England) was then detected using enhanced chemiluminescence (ECL) reagent (Bio-Rad, 170506, Hercules, CA, USA).

After harvesting LGALS1-silenced (siLGALS1) and control cells, genomic DNA was removed with RQ1 DNase (Promega, Madison, WI, USA) to obtain total RNA. Qualified RNA, as measured by absorbance at 260 nm/280 nm (A260/A280) using SmartSpec Plus (BioRad, Hercules, USA), was used in the following experiments. One $\mu$g RNA was captured from the siLGALS1 and control samples (three replicates of each) using VAHTS mRNA capture Beads (Vazyme, N401, Nanjing, China). RNA-seq libraries were prepared on the Illumina® Platform using polyadenylated RNAs and oligo (dT) with the KAPA Stranded Kit. After size selection (300–500 bp), the cDNA products were processed, stored at –80 °C, then used for high-throughput sequencing (RNA-seq) with the Illumina NovaSeq 6000 system (Illumina, San Diego, CA, USA). Following that, 150 nt paired-end sequencing reads were obtained.

Raw RNA-seq reads were filtered by removing N-containing reads, adaptors and low-quality bases, and discarding short reads (16 nt) using the FASTX-Toolkit (Version 0.0.13) (New York, USA), resulting in high-quality data. High-quality reads were then aligned onto the human GRCh38 genome by HISAT2 (Baltimore, MD, USA) [20], with no more than 4 mismatches. Aligned reads with unique genomic location were used to calculate the gene expression level as FPKM (fragments per kilobase of transcript per million fragments mapped) [21]. The DESeq2 package [22] from R Bioconductor was then used to predict DEGs, with criteria set as a corrected p-value (FDR) 0.05 and fold-change $>$2. Selected DEGs with altered expression levels were validated by RT-qPCR experiments (primers shown in Supplementary Table 1). Glyceraldehyde-3-phosphate dehydrogenase (GAPDH) was used as the internal control.

The ABLas pipeline was used to predict and quantify alternative splicing events (ASEs) and galectin-1-regulated ASEs in siLGALS1 and Ctrl groups as described previously [23, 24]. Based on the splicing reads, 9 types of ASEs were analyzed: alternative 3’ splice site (A3SS), alternative 5’ splice site (A5SS), cassette exon, exon skipping (ES), A3SS&ES, A5SS&ES, mutually exclusive exons (MXE), mutually exclusive 5’UTRs (5pMXE), and mutually exclusive 3’UTRs (3pMXE). Student’s t-test was used to calculate the statistical significance of the ASEs’ ratio differences in order to identify regulated ASEs (RASEs). ASEs with a p-value 0.05 were considered RASEs. RT-qPCR experiments were conducted to validate the altered expression levels for selected RASEs. Primers for RT-qPCR are shown in Supplementary Table 1.

Functional cluster analyses were conducted using the Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) databases to analyze DEGs and RASGs for enriched functional pathways. Briefly, DEGs and RASGs were mapped to all terms in the GO database to count gene numbers, and the hypergeometric distribution test was used to obtain the significance of each GO term for DEGs and RASGs’ against a background of GO annotation for the whole genome. The analysis method was the same as for the KEGG analysis mentioned above.

The iRIP-seq data for galectin-1 was downloaded from the NCBI GEO database with accession number GSE147565. This was used to investigate how galectin-1 regulates DEGs and RASEs. Galectin-1-bound peaks and genes were re-analyzed according to a previously published method [16]. Analyses of identified DEGs and RASGs were performed after the overlapping peaks and genes were obtained from two iRIP-seq replicates.

Data were analyzed with GraPhPad prism application software (San Diego, USA) and were presented as the mean $\pm$ standard deviation. Statistical differences between groups were evaluated by Student’s t-test, or multiple comparison single factor difference analysis for two independent samples. * p-value 0.05; ** p-value 0.01; *** p-value 0.001; **** p-value 0.0001. Hypergeometric test was performed to analyze the overlapping genes.

Our previous study demonstrated that silencing of galectin-1 could inhibit capillary tube formation in HUVECs [16]. However, it is not known how galectin-1 regulates angiogenesis, and therefore the LGALS1 knockdown experiment was repeated using siLGALS1. Successful knockdown of LGALS1 was confirmed by RT-qPCR and WB experiments (Fig. 1A,B). RNA-seq was then performed to identify LGALS1-regulated genes at the transcriptional and post-transcriptional levels. Three biological replicates with high-quality sequencing data were obtained for siLGALS1 and negative control (Ctrl) samples (Supplementary Table 2). After aligning quality-filtered reads to the human genome and calculating FPKM for the gene expression level, LGALS1 was found to be highly expressed, with an average FPKM value of 1731. This was only slightly lower than expression of the house-keeping genes ACTB and GAPDH (Fig. 1C). SiLGALS1 samples had a lower FPKM value compared to Ctrl samples (Fig. 1C). Principal component analysis (PCA) for all expressed genes revealed a clear separation between siLGALS1 and Ctrl samples, with a short inner distance for the three biological replicates (Fig. 1D). The sharp contrast indicated the global changes in transcriptome profiles following siLGALS1. These results indicate that galectin-1 is highly expressed in HUVECs and may play a critical role in regulating the transcriptome in these cells.

Fig. 1.

SiLGALS1 has a global influence on the transcriptome profile of HUVECs. (A) Bar plot showing the relative expression levels of LGALS1 in siLGALS1 and Ctrl samples, as determined by the RT-qPCR method. N = 3. (B) Western blot result showing a decreased level of LGALS1 protein in siLGALS1 samples compared with Ctrl samples. N = 3. (C) Scatter plot showing the correlation between siLGALS1 and Ctrl samples using the FPKM values for all expressed genes. The FPKM values for LGALS1 are shown in the bar plot. (D) PCA result showing that siLGALS1 samples are clearly separated from Ctrl samples at the first principal component. *** p-value 0.001.

We used siLGALS1 to identify dysregulated genes that could participate in angiogenesis. With the criteria set to a 2-fold change and 0.05 false discovery rate (FDR), 1451 DEGs were identified using the edgeR package [25], comprising 604 up-regulated and 847 down-regulated genes (Fig. 2A, Supplementary Table 3). A hierarchical cluster heatmap of all DEGs showed consistent expression changes for the three replicates (Fig. 2B). Interestingly, more down-DEGs were detected than up-DEGs, suggesting that galectin-1 is likely to promote gene expression. We then investigated the functions of up- and down-DEGs using GO and KEGG databases. The top 5 enriched GO biological process (BP) terms included muscle filament sliding, the adenylate cyclase-modulating G-protein coupled receptor signaling pathway, sodium ion transmembrane transport, the G-protein coupled receptor signaling pathway, and cell-cell signaling (Fig. 2C). We next focused on the galectin-1-promoting genes (down-DEGs). Functional enrichment analysis revealed that the top 5 terms were signal transduction, angiogenesis, ECM organization, inflammatory response, and positive regulation of angiogenesis (Fig. 2D). Down-DEGs were significantly enriched in angiogenesis and in positive regulation of angiogenesis pathways. These two terms included 18 and 12 genes, respectively. A heatmap of the two terms showed that all of the genes were down-regulated in siLGALS1 samples (Fig. 2E). To further validate these genes, 10 were selected from the two terms and RT-qPCR was performed. The results confirmed that all of the genes were down-regulated in siLGALS1 samples (Fig. 2F). These DEG results demonstrate that siLGALS1 decreases the expression of angiogenesis-related genes, which could be associated with the promotion of angiogenesis.

Fig. 2.

SiLGALS1 reduces the expression of angiogenesis-related genes. (A) Volcano plot showing the expression pattern of DEGs. Red and blue points indicate up- and down-regulated DEGs in siLGALS1 samples, respectively. (B) Hierarchical clustering heatmap showing the expression pattern of all DEGs. (C) Bubble plot showing the top enriched GO BP pathways of the up-regulated DEGs. (D) Bubble plot showing the top enriched GO BP pathways of the down-regulated DEGs. (E) Hierarchical clustering heatmap showing the down-regulated genes selected from the two angiogenesis pathways. (F) Bar plot showing RT-qPCR results for the 10 genes selected from the two angiogenesis terms. N = 3. **** p-value 0.0001.

Our previous study found extensive interactions between galectin-1 and thousands of RNA molecules [16], hinting at potential functions for galectin-1 in RNA regulation. AS is an important post-transcriptional regulation mechanism in all eukaryotes. Galectin-1-regulated AS events in HUVECs were systematically analyzed. The ABLas program [24] was applied to identify ASEs and calculate statistical differences between siLGALS1 and Ctrl samples. In total, we identified 2327 RASEs that were significantly different between the two groups. After classifying these RASEs into 10 types, ES and IR events were found to be the most common RASE types (Fig. 3A). More ES (ES, A3SS&ES, and A5SS&ES) and less cassette exon (opposite to ES) events were identified in siLGALS1 samples (Fig. 3A). Therefore, siLGALS1 showed a trend to induce ES events in HUVECs. Moreover, the number of included RASEs (up, 1523) was much higher than the number of excluded RASEs (down, 804), suggesting that galectin-1 has a tendency for AS regulation. A hierarchical clustering heatmap of all RASE ratios demonstrated a clear difference between the siLGALS1 and Ctrl groups (Fig. 3B), suggesting the AS events for these genes were consistently regulated by siLGALS1. We next assessed the enriched functions of these RASGs. The top 5 enriched GO BP terms were found to be autophagy, cellular protein metabolic process, viral reproduction, positive regulation of protein targeting to membrane, and DNA repair (Fig. 3C). KEGG pathway analysis also revealed that RASGs were highly enriched in focal adhesion and VEGF signaling pathways (Fig. 3D). Previous studies have reported that RASGs play important roles in angiogenesis [26, 27]. Of the genes involved in focal adhesion and VEGF signaling pathways, 8 RASEs were selected to confirm the altered AS ratios using RT-qPCR experiments (ADRM1, EGFL7, FGFR1, FLNC, SPHK1, MAPK3, FLT1, and MAPK14). Specific primers for these RASEs (shown in Supplementary Table 1) were designed according to previously published methods [28]. The RNA-seq data revealed that 7 of the 8 RASEs showed consistent ratio changes, with MAPK3, FGFR1, FLNC, and SPHK1 presented as examples in Fig. 3E,F and Supplementary Fig. 1A,B. Together, these results indicate that galectin-1 is a strong regulator of the AS patterns for angiogenesis-associated genes in HUVECs.

Fig. 3.

Galectin-1 is a strong regulator of the AS profiles of VEGF signaling pathway genes in HUVECs. (A) Bar plot showing the number of each RASE type. Up and down indicates the inclusion and exclusion pattern, respectively, of each RASE type in siLGALS1 samples. (B) Hierarchical clustering heatmap showing the ratio patterns of all RASEs. (C) Bubble plot showing the top enriched GO BP pathways of all RASGs. (D) Bubble plot showing the top enriched KEGG pathways of all RASGs. (E) Reads plot and RT-qPCR results of ES event for MAPK3 in the VEGF signaling pathway. N = 3. (F) Reads plot and RT-qPCR results for A3SS event of FGFR1 in the VEGF signaling pathway. N = 3. * p-value 0.05; ** p-value 0.01; **** p-value 0.0001.

Based on the above results, we deduced that galectin-1 has significant influence on both the expression and AS of genes in HUVECs. We next investigated whether any genes were co-regulated by siLGALS1 at both the expression and alternative splicing levels. A total of 54 genes were identified both in the 1451 DEGs and in the 1577 RASGs (Fig. 4A). Functional analyses of these 54 overlapping genes revealed they were enriched in the sphingolipid metabolism pathway (Fig. 4B, p-value 0.01), and included GBA2 and ARSA. To further explore the RNA level changes of RASGs, we investigated their expression patterns using a volcano plot. This showed that 756 genes (47.94% of all RASGs) were significantly dysregulated (FDR 0.05), including 356 up- and 400 down-regulated genes (Fig. 4C), without regard for the fold-change threshold. Interestingly, RASGs involved in the VEGF signaling pathway were also found to be dysregulated, including SPHK1, PLA2G4A, PLA2G4C, and HRAS (Fig. 4C). The expression patterns for ADRM1 and EGFL7 are presented in Fig. 4D,E to show their dysregulated RASEs and expression levels. Another important gene is FLT1, the receptor for VEGF and target for anti-angiogenic therapy [29]. FLT1 was found to be significantly dysregulated by siLGALS1 at the AS and expression levels (Supplementary Fig. 2). Based on previous findings that AS regulates gene expression [30], we further hypothesized that galectin-1 could regulate the expression and AS levels of angiogenesis-associated genes in HUVECs, but that the expression level changes were small because they did not reach the 2-fold threshold.

Fig. 4.

Analysis of the Overlap between DEGs and RASGs. (A) Venn diagram showing the number of overlapping genes between DEGs and RASGs. (B) Bubble plot showing the enriched GO BP pathways of the overlapping genes from A. (C) Volcano plot showing the expression patterns of RASGs. Red points indicate RASGs with significant FDR (0.05) for their expression levels. Several genes from the VEGF signaling pathway are labeled. (D) Reads plot, RT-qPCR result, and expression level for the ES event of EGFL7. The AS ratio, FPKM value, and RT-qPCR result of RASE are shown in the right panel. N = 3. (E) Reads plot, RT-qPCR result, and expression levels for the ES event of ADRM1. N = 3. *** p-value 0.001; **** p-value 0.0001.

We previously investigated the RNA binding ability of galectin-1 and its regulatory roles in angiogenesis using the iRIP-seq method [16]. Here, we focused on interpreting the relationship between galectin-1-bound RNAs and the galectin-1-regulated transcriptome. Using the overlapping peaks and genes bound by galectin-1 from two replicates of iRIP-seq, we identified 36 DEGs that were bound by galectin-1 (Fig. 5A). Although the overlap between DEGs and galectin-1-bound genes was not significant, 36 genes were found to be enriched in transcription regulation-associated pathways (Fig. 5B), including HR, ECM1, AHR, FOSL1, ZNF860, CREBRF, ATF6B, and SNAI1. These transcription-associated genes play critical roles in gene expression, indicating that galectin-1 may regulate gene expression profiles by affecting the expression level of these genes. Furthermore, galectin-1 may directly influence galectin-1-bound genes by modulating their AS patterns. The overlap analysis revealed that 352 RASGs were bound by galectin-1 (Fig. 5C, p-value = 5.646611 $\times$ 10${}^{-64}$, hypergeometric test). This result suggests that galectin-1 has a tendency to regulate AS, and this may be achieved by binding to pre-mRNAs. Functional enrichment analyses of the 352 overlapping genes revealed they were highly enriched in translation and alternative splicing-associated pathways (Fig. 5D). Two angiogenic pathways were also enriched (Fig. 5D). Proteins involved in translation and RNA splicing pathways are mainly RBPs that play crucial roles in transcriptional and post-transcriptional processes [31]. The fact that galectin-1 regulates RBPs at the AS level provides further evidence that galectin-1 is critical to transcriptome profiles. Finally, RTN4 and SRSF10 were selected from the angiogenesis and RNA splicing pathways, respectively, to illustrate the galectin-1-bound patterns and AS changes (Fig. 5E,F, and Supplementary Fig. 3A,B).

Fig. 5.

Results of Integration Analysis Between the LGALS1-regulated Transcriptome and RNA Interactome. (A) Venn diagram showing the overlapping genes between DEGs and LGALS1-bound RNAs. (B) Bubble plot showing the functionally enriched BP pathways for the overlapping genes shown in (A). (C) Venn diagram showing the overlapping genes between RASGs and LGALS1-bound RNAs. (D) Bubble plot showing the functionally enriched BP pathways for the overlapping genes in (C). (E) Reads plot showing the AS pattern of SRSF10. The green background indicates the region of the AS event. (F) Reads plot showing the LGALS1-bound profile of SRSF10. The green background indicates the region of LGALS1-bound peaks.