RNA sequencing reveals the transcriptomic landscape and alternative splicing events induced by LGALS1 silencing in non-small cell lung cancer

Background. Non-small cell lung cancer (NSCLC) is a common clinical cancer with high mortality. The lectin galactoside-binding soluble 1 (LGALS1) is an RNA-binding protein (RBP) involved in NSCLC progression. Alternative splicing (AS) is a vital function of RBPs that contributes to tumor progression. It is unknown whether LGALS1 regulates NSCLC progression through AS events. Objectives. To profile the transcriptomic landscape and LGALS1-regulated AS events in NSCLC. Materials and methods. The A549 cells either with silenced LGALS1 (siLGALS1 group) or without them (siCtrl group) were subjected to RNA sequencing; differentially expressed genes (DEGs) and AS events were discovered and then the AS ratio was validated using reverse transcription-quantitative polymerase chain reaction (RT-qPCR). Results. High LGALS1 expression indicates poor overall survival (OS), first progression (FP) and post-progression survival (PPS). A total of 225 DEGs were identified, including 81 downregulated and 144 up-regulated in the siLGALS1 group compared to the siCtrl group. Differentially expressed genes were mainly enriched in interaction-related Gene Ontology (GO) terms and involved in cGMP-protein kinase G (PKG) and calcium signaling pathways. The RT-qPCR validation showed that the expressions of ELMO1 and KCNJ2 were upregulated, while HSPA6 was downregulated after LGALS1 silencing. The expressions of KCNJ2 and ELMO1 were upregulated to a peak at 48 h after LGALS1 knockdown, while HSPA6 expression decreased, after which their expressions returned to baseline. The overexpression of LGALS1 rescued the elevation in KCNJ2 and ELMO1 expression, and decrease in HSPA6 expression induced by siLGALS1. A total of 69,385 LGALS1-related AS events were detected, which produced 433 upregulated and 481 downregulated AS events after LGALS1 silencing. The LGALS1-related AS genes were mainly enriched in the apoptosis and ErbB signaling pathways. The LGALS1 silencing led to a decrease in the AS ratio of BCAP29 and an increase in CSNKIE and MDFIC. Conclusions. We characterized the transcriptomic landscape and profiled AS events in A549 cells following LGALS1 silencing. Our study provides abundant candidate markers and new insights into NSCLC.


Background
Lung cancer is a frequent malignant disease with high mortality, 1 and non-small cell lung cancer (NSCLC) accounts for the majority of lung cancer cases.This disease can be classified as lung adenocarcinoma (LUAD) or lung squamous cell carcinoma (LUSC). 2 Despite considerable progress in early diagnosis, there is still a high recurrence rate after treatment, 3 and the pathogenesis of NSCLC remains largely unknown.Thus, it is urgent to identify new driver genes and elucidate their mechanism to provide clinical advancements.
RNA-binding proteins (RBPs) play an important role in translation and gene regulation, 4,5 driving the development of many diseases and cancers. 6The RBP lectin galactoside-binding soluble 1 (LGALS1) is a β-galactosidebinding protein that has a carbohydrate recognition domain, 7 and an extensive body of research has revealed that LGALS1 participates in the progression of multiple tumors. 8,9For example, the upregulated expression of LGALS1 promotes the NSCLC progression by interacting with non-SMC condensin I complex subunit G (NCAPG). 10 Nevertheless, the mechanism and manner by which LGALS1 participates in the NSCLC progression are largely unknown.
Alternative splicing (AS) is a vital function of RBPs.It is reported that more than 90% of genes undergo AS to regulate gene expression, which ultimately results in proteome diversity, 11 but can also increase cancer risk. 12he LGALS1 participates in pre-mRNA splicing, and it has been reported to regulate gene expression. 13The LGALS1 was shown to reduce CaV1.2 calcium channel currents resulting in the regulation of vascular constriction via AS. 14owever, whether LGALS1 AS may regulate the NSCLC pathological processes remains unclear.

Objectives
We aimed to evaluate whether LGALS1 exerts its regulatory effects through the regulation of gene expression and AS to discover new prognostic and therapeutic targets for NSCLC.For this purpose, we knocked down LGALS1 in A549 cells and analyzed the differentially expressed genes (DEGs) and LGALS1-regulated AS events associated with NSCLC.

Public databases
The Kaplan-Meier plotter (http://kmplot.com/analysis/index.php?p=service&cancer=lung) was used to analyze the relationship between LGALS1 expression and survival, including overall survival (OS), first progression (FP) and post-progression survival (PPS), based on the set of NSCLC samples.Data were processed based on the study by Győrffy et al. 15 The Cancer Genome Atlas (TCGA, http://cancergenome.nih.gov) was used to analyze the expression pattern of LGALS1 in NSCLC.

Cell culture and transient transfection
Human A549 (#CL-0016) and H1650 (#CL-0166) cell lines were purchased from China Procell Life Science & Technology (Wuhan, China).Normal human lung epithelial BEAS-2B cells and human H1299 cells were provided by our laboratory.Cells were cultured in a minimum essential medium (OPTI-MEM) supplemented with 10% fetal bovine serum (FBS) and 1% penicillin/streptomycin, and maintained at 37°C with 5% CO 2 .
Transient transfection of siRNA sequences and shRNA vectors was performed using Lipofectamine™ 2000 (Invitrogen, Waltham, USA), according to the manufacturer's instructions.Briefly, 1 day before transfection, A549 cells were seeded in 6-well plates (Corning, Tewksbury, USA) at a density of 5×10 5 cells per well.The siRNA sequences were mixed with OPTI-MEM at a ratio of 5 μL:45 μL.The transfection reagent Lipofectamine™ 2000 (5 μL) was also diluted with 45 μL of OPTI-MEM.Next, 50 μL of diluted siRNA and 50 μL of diluted transfection reagent were mixed and incubated for 15 min at room temperature.After that, the entire mixture was added into the A549 cell sample to incubate for another 24 h until subsequent assays were performed.
After imaging, the expression of LGALS1 relative to GAPDH was calculated using ImageJ software (National Institutes of Health, Bethesda, USA).

RNA sequencing and data quality control
The genomic material was removed from A549 cells, and total RNA was separated using the Trizol reagent (Invitrogen) in both the siLGALS1 group (n = 3) and siCtrl group (n = 3).The RNA was purified using an RNA purification kit (Tiangen, Beijing, China) and the concentration and quality of each RNA sample was determined using SmartSpec Plus (BioRad, Hercules, USA) by calculating the absorbance at 260 nm/280 nm (A260/A280).The RNA sequencing (RNA-seq) libraries were prepared with the KAPA Stranded mRNA-seq Kit for Illumina® Platforms (#KK8401; KAPA Biosystems, Boston, USA) according to the manufacturer's instructions.Briefly, 1 μg of intact total RNA was used for library construction.Poly(A) RNA was captured with magnetic oligo-dT beads, which were then resuspended in 1X Fragment, Prime and Elute Buffer to elute the captured RNA.Next, the RNA was fragmented to the desired size by heating in the presence of Mg 2+ .Then, first-strand cDNA was synthesized with random oligo-dT primers followed by second-strand cDNA synthesis, and double-stranded cDNA was manufactured with RNA while marking the second-strand with dUTP.Next, dAMP was added to the 3'-end of double-strand cDNA fragments to obtain the 3'-dAMP library fragments, followed by 3'-dTMP adapter ligation, and the adapter-ligated library DNA was amplified using polymerase chain reaction (PCR).Finally, library fragment size distribution was confirmed with electrophoresis, and library concentration was determined with quantitative PCR (qPCR).Then, RNA-Seq libraries were sequenced on the Illumina NovaSeq 6000 platform using paired-end 150 nt reads.Raw reads were filtered using the FASTX-Toolkit (v.0.0.13)(http://hannonlab.cshl.edu/fastx_toolkit/index.html) to discard low-quality reads.The alignment of clean reads was then performed by tophat2 16 mapped to the GRCh38 human reference genome.Uniquely mapped results were used to obtain the read count of a gene as expression and then normalized using the FPKM algorithm. 17

DEGs analysis
The edgeR package 18 was used to distinguish DEGs between the siLGALS1 group and the siCtrl group employing criteria of the false discovery rate (FDR) less than 0.05 and fold change ≥2 or ≤0.5.To predict functional categories of siLGALS1-related genes, Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis were performed using the KOBAS 2.0 server. 19

RBP-related AS analysis
The ABLas pipeline was applied for identifying and quantifying the AS events and regulated AS events (RASEs). 20n short, 10 types of normative AS events were detected in 2 groups based on the splice junction reads, including cassette exons (CE), exon skipping (ES), mutual exclusive exon skipping (MXE), alternative 3' splice site (A3SS), A3SS and ES (A3SS & ES), MXE combined with an alternative polyadenylation site (3pMXE) or an alternative 5' promoter (5pMXE), alternative 5' splice sites (A5SS), A5SS and ES (A5SS & ES), and intron Retained (introR).Next, to determine RBP-related RASEs, the difference in AS events between the siLGALS1 group and the siCtrl group was evaluated using Student's t-test.The difference in AS events with a p-value cutoff corresponding to an FDR cutoff of 5% was considered RBP-related RASE.

RTqPCR verification of DEGs and AS
Primers of DEGs were designed using Primer Premier v. 6.0 (http://www.premierbiosoft.com/primerdesign/)validated with qPCR.For reverse transcription (RT)-qPCR validation (each gene was repeated 9 times), M-MLV Reverse Transcriptase (Vazyme Biotech, Nanjing, China) was used for RNA reverse transcription into cDNA.The SYBR Green PCR Reagents Kit (Yeasen Biotechnology, Shanghai, China) was used for the PCR reaction with 3 repetitions on the StepOne Realtime PCR System.The gene expression level was normalized to GAPDH using the 2 −ΔΔCT method. 21o verify AS events (each gene was repeated 9 times), a boundary-spanning primer was used to detect alternative isoforms for the sequence containing the junction of the alternative exon (according to "model exon" to detect model splicing or "altered exon" to detect altered splicing) and constitutive exon.The sequence containing the constitutive exon was detected with an opposing primer.The primer sequences are listed in Supplementary Table 1.

Statistical analyses
GraphPad Prism v. 9.0 (GraphPad Software, San Diego, USA) was used for statistical analysis.Data were obtained from at least 3 independent samples and indicated as means ± standard deviation (M ±SD).Considering that the sample size is less than 10, the data are all analyzed using a nonparametric test.The statistical comparison between the siLGALS1 group and the siCtrl group was performed using Mann-Whitney (M-W) test when the total sample >7; otherwise, the t-test was performed (data conformed to the normal distribution and variance homogeneity).The comparison between 3 or 4 groups was determined with Kruskal-Wallis (K-W) test followed by Dunn's multiple comparison test.A p-value <0.05 was considered statistically significant.

High expression of the LGALS1 gene corresponds to a poor prognosis of NSCLC
To determine whether LGALS1 expression was associated with NSCLC progression, the Kaplan-Meier plotter was used for drawing OS, FP and PPS curves based on the LGALS1 mRNA expression, which was stratified into low-and high-expression groups.The results showed that patients in the LGALS1 high-expression group exhibited a significant decrease in the OS (log-rank p = 6.6e-13),FP (log-rank p = 0.044) and PPS (log-rank p = 0.00026) compared with the LGALS1 low-expression group (Fig. 1), suggesting that high expression of the LGALS1 gene correlated with poor prognosis of NSCLC.However, unexpectedly, according to TCGA, LGALS1 expression was not statistically different in tumor and control tissues (Supplementary Fig. 1).

LGALS1 silencing alters the transcriptomic landscape of A549 cells
First, we performed RT-qPCR to determine LGALS1 expression in 4 NSCLC cell lines and normal cells.The results showed that compared with the BEAS-2B cell, LGALS1 was significantly overexpressed in 3 NSCLC cell lines, with its highest expression observed in A549 cells (K-W test, H = 9.462, p = 0.0067) (Fig. 2A).Therefore, to assess the molecular mechanisms of LGALS1 in NSCLC, we silenced LGALS1 expression in A549 cells using siRNA.The mRNA (M-W test, U = 0, p < 0.0001) and protein (t-test, t = 4.824, degrees of freedom (df) = 4, p = 0.0485) expression of LGALS1 in the siLGALS1 group were significantly lower than in the siCtrl group, suggesting that the knockdown was successful (Fig. 2B,C).Then, A549 cells containing silenced LGALS1 were subjected to transcriptome sequencing.The LGALS1 expression was also significantly reduced in the RNA-seq analysis, supporting the successful LGALS1 knockdown (t-test, p = 7.52e-14) (Fig. 2D).The Q30 and GC content of all 6 samples were larger than 94.9% and 50%, respectively (Supplementary Table 2), indicating that the quality control of transcriptome sequencing was good.Principal component analysis (PCA) showed that samples from the siLGALS1 group and the siCtrl group could be distinguished from each other (Fig. 2E), suggesting that the interference in LGALS1 expression leads to an alteration in transcriptomic profile in A549 cells.
To further investigate the genes regulated by LGALS1, DEGs (|fold change| >2 and FDR < 0.05) induced by LGALS1 silencing were characterized.A total of 225 DEGs were identified, including 144 upregulated and 81 downregulated DEGs in the siLGALS1 group compared to the siCtrl group (Fig. 2F).After hierarchical clustering of all the DEGs, a robust distinction was observed between the 2 groups, and a high coincidence among the changes was also observed in the 3 replicates (Fig. 2G).These results suggest that LGALS1 silencing alters the transcriptomic landscape of A549 cells.

Functional clustering analysis of DEGs
Next, we used the KOBAS 2.0 server to predict the molecular functions and metabolic pathways of up-and downregulated DEGs, respectively.The results showed that upregulated DEGs only affected small molecule metabolic processes and DNA transcription, and were mainly involved in cell adhesion molecules and chemical carcinogenesis pathways (Fig. 3A,B).In addition, downregulated DEGs were enriched to interaction-related GO terms such as extracellular matrix, steroid binding, laminin binding, and ankyrin binding (Fig. 3C).The KEGG analysis found that downregulated DEGs were mainly associated with protein digestion and absorption, cGMP-protein kinase G (PKG) signaling pathway, and calcium signaling pathways (Fig. 3D).These results imply that LGALS1 silencing alters numerous physiological processes in A549 cells via DEGs.
Given that gene regulation is time-dependent, we examined the expression of DEGs (KCNJ2, HSPA6 and ELMO1) at different timepoints after LGALS1 knockdown to further validate the regulation of DEGs by LGALS1.Gene expression in the shCtrl group did not change with time, while LGALS1 mRNA expression gradually decreased in the shLGALS1 group compared with the shCtrl group, reaching a minimum at 48 h, after which the transient knockdown effect gradually disappeared (Fig. 5A).As expected, shLGALS1 regulation of DEG expression showed temporal effectiveness for 48 h.The KCNJ2 and ELMO1    expressions were upregulated to a peak at 48 h after LGALS1 knockdown, while HSPA6 expression decreased to a peak, after which their expressions all returned to nondifferential levels from the shCtrl group (Fig. 5A).The shRNA-LGALS1 lentiviral vector also inhibited LGALS1 expression at the protein level (Fig. 5B).In addition, for further validation of the genes regulated by LGALS1, a rescue experiment was conducted.We observed that overexpression of LGALS1 significantly reduced the elevation in KCNJ2 and ELMO1 expression and minimized the decrease in HSPA6 expression induced by siLGALS1 alone (Fig. 5C).Taken together, these results implicate that LGALS1 can regulate the expressions of DEGs, such as KCNJ2, HSPA6 and ELMO1.

Identification of LGALS1-mediated AS events in A549 cells
Previous studies have demonstrated that AS is involved in lung cancer development. 27Thus, we further explored LGALS1-mediated AS events.A total of 69,385 AS events were identified, of which 18,587 were known AS, accounting for 26.79%, and 50,798 were novel AS, accounting for 73.21% (Fig. 6A, Supplementary Table 3).We found that both known and novel AS events contained 10 types (Fig. 6A, Supplementary Tables 4 and 5), of which A3SS and A5SS exhibited the largest share of the total AS events, corresponding to 26.03% and 24.95%, respectively (Fig. 6A, Supplementary Table 6).

Functional enrichment analysis of RASGs
To explore the biological function of LGALS1-regulated AS, RASGs were subjected to GO and KEGG enrichment analyses.The GO analysis revealed that RASGs were preferentially enriched in molecular functions related to transcriptional regulation, such as DNA-dependent transcription, regulation of transcription from RNA polymerase, DNA-dependent regulation of transcription, and gene expression (Fig. 7A).In addition, KEGG pathway analysis found that RASGs were mainly enriched in ErbB signaling and apoptosis pathways (Fig. 7B).These results demonstrated that LGALS1 may be involved in the regulation of cell apoptosis through regulating RASGs.

Validation of the LGALS1-mediated RASGs
Given that RASGs may be associated with NSCLC progression, we validated LGALS1-mediated RASGs using qPCR.Out of 9 detected AS events, the results of 3 A3SS AS events validated using RT-qPCR were consistent with the RNA-seq data, namely BCAP29, CSNKIE and MTFP1 (Fig. 8).Among them, the AS ratio of BCAP29 (M-W test, U = 13, p = 0.0142) and MTFP1 (M-W test, U = 11, p = 0.0078) were significantly increased in the siLGALS1 group compared to the siCtrl group, whereas CSNKIE (M-W test, U = 17, p = 0.0193) was found to be opposite.The RT-qPCR validation results of the remaining 6 AS events are shown in Supplementary Fig. 2.These results further confirm the LGALS1-related AS and RASE data.

Discussion
Despite extensive research regarding NSCLC, its mortality rate is still high, and more studies on NSCLC are warranted to understand the mechanism of the disease.In the present study, we have revealed insights into the transcriptomic landscape and AS dataset in NSCLC, driven by LGALS1 through RNA-seq.Many DEGs and RASGs in this study provide new insights and abundant potential candidate markers for NSCLC that could be targeted therapeutically.
As an RBP, LGALS1 is a β-galactoside-binding protein and is alternatively known as galectin-1. 28Many studies have demonstrated that intervention in the expression of LGALS1 plays a vital role in the progression of numerous human cancers.For instance, a recent study showed that LGALS1 is involved in NSCLC progression by interacting with NCAPG, and the knockdown of LGALS1 significantly suppresses proliferation, migration and invasion of A549 and H1299 cells. 10In oral cancer, LGALS1 knockdown suppressed cell proliferation by arresting cells at S phase and inhibited wound healing and migration. 29Furthermore, in certain brain tumors, silencing of LGALS1 downregulated the expression of genes involved in cell cycle progression, resulting in an accumulation of G2/M phase cells and eventually inhibiting tumor progression. 30There is also an extensive body of literature showing that LGALS1 participates in lung cancer progression.For example, LGALS1 is overexpressed in LUAD cells and induces their growth and invasive ability. 8The LGALS1 upregulated p38 MAPK, ERK and cyclo-oxygenase-2 expression to promote lung cancer progression. 9In the current study, we demonstrated that silencing of LGALS1 led to the aberrant expression of 225 genes in A549 cells and involved alterations in transcription and binding-related pathways, suggesting that LGALS1 may function in NSCLC.Therefore, further research is needed on RNA-binding properties and related biological functions of LGALS1.
Alternative splicing events are closely associated with the function of RBPs and cancer progression, 31 with previous studies showing that AS plays a crucial role in the development of lung cancer.For instance, the AS status of BIN1 with exon 12A inclusion (the BIN1+12A isoform) could recover the activity of tumor cells through the regulation of BIN1 in NSCLC. 32Furthermore, ESRP1 regulates the chemosensitivity of lung cancer cells through AS by inhibiting TGF-β/Smad signaling in SCLC patients. 33Moreover, VEGFxxxb family members encoded AS in VEGF-A, and VEGF165b/VEGF165 are positively correlated with lymph node metastasis in NSCLC. 34e current study showed that a total of 914 AS events were triggered after LGALS1 silencing, and the corresponding RASGs were implicated in ErbB signaling and apoptosis pathways.Encouragingly, some previous results support our observations.Wang et al. demonstrated that LGALS1 modulates vascular constriction by regulating AS of the cav1.2calcium channel. 14Therefore, LGALS1mediated AS events play an important role in NSCLC progression.
In addition, we discovered that LGALS1-related RASGs were significantly enriched in the apoptosis and ErbB signaling pathways.Apoptosis, or programmed cell death, plays a critical role in cancer occurrence, progression and drug resistance.Inhibition of LGALS1 in melanoma cells can restore the viability of apoptotic T lymphocytes 35 while targeting LGALS1 by shikonin, and its derivatives were shown to induce apoptosis and autophagy in colorectal carcinoma cells. 36Moreover, the ErbB signaling pathway is involved in cancer progression, for example, in cervical cancer. 37Finally, EGFR is a member of the ErbB receptor family, which regulates epithelial cell growth and survival, and its high expression or abnormal activation is correlated with drug resistance and even tumor progression in multiple cancers, such as breast cancer 38 and NSCLC. 39herefore, we suspect that LGALS1 may be involved in the NSCLC progression through the apoptosis and ErbB signaling pathways.

Limitations
There are some limitations of this study, one of which is that our results were not validated in animal models or clinical samples.Second, there are only 3 biological repeats for transcriptomic sequencing.Finally, no molecular experiments were performed to verify the regulatory relationship between LGALS1 on candidate target genes and AS events.Therefore, in the future, we will further explore the mechanism of LGALS1 involvement in NSCLC progression through a series of in vitro and in vivo experiments to address the above limitations.

Conclusions
In the current study, we characterized the transcriptomic landscape of A549 cells following LGALS1 silencing, identifying 225 DEGs that responded significantly to LGALS1 silencing.The expression patterns of ELMO1, KCNJ2 and HSPA6 were consistent between both RNA-seq and RT-qPCR experiments, and the expression of ELMO1 and KCNJ2 was upregulated, whereas HSPA6 expression was downregulated after LGALS1 silencing.In addition, LGALS1 silencing led to 914 RASEs consisting of 10 types of AS residing in 759 genes.These LGALS1-mediated RASGs were mainly enriched in apoptosis and ErbB signaling pathways.Our results show potential novel molecular mechanisms of RBP LGALS1-mediated AS in NSCLC development.

Fig. 1 .
Fig. 1. High expression of lectin galactoside-binding soluble 1 (LGALS1) indicates poor prognosis in non-small cell lung cancer (NSCLC).The pattern of overall survival (OS) (A), first progression (FP) (B) and post-progression survival (PPS) (C) of NSCLC patients with high and low LGALS1 expression levels HR -hazard ratio.

Fig. 3 .
Fig. 3. Functional clustering analysis of differentially expressed genes (DEGs).A. The only 2 Gene Ontology (GO) enrichment terms of upregulated DEGs between silenced lectin galactoside-binding soluble 1 (siLGALS1) and siCtrl group; B. The top 10 enriched Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways for upregulated DEGs between siLGALS1 and siCtrl group.The top 7 GO enrichment terms (C) and top 10 KEGG pathways (D) for downregulated DEGs between siLGALS1 and siCtrl group

Fig. 5 .
Fig. 5. Differentially expressed genes (DEGs) regulated by lectin galactoside-binding soluble 1 (LGALS1).A. Reverse transcription-quantitative polymerase chain reaction (RT-qPCR) was performed to detect mRNA expression of LGALS1, KCNJ2, HSPA6, and ELMO1 at different timepoints following LGALS1 knockdown using shRNA; B. Western blot was performed to detect the protein expression of LGALS1 at different timepoints after LGALS1 knockdown using shRNA; C. A rescue experiment was used to verify the regulatory role of LGALS1 on KCNJ2, HSPA6 and ELMO1 mRNA expression

Fig. 7 .
Fig. 7. Functional enrichment analysis of regulated alternative splicing genes (RASGs). A. The top 10 enriched Gene Ontology (GO) terms for RASGs; B. The top 10 enriched Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways for RASGs