The transcriptome dataset GSE16515 comprising 52 patients with pancreatic cancer was obtained from the Gene Expression Omnibus (GEO) dataset (https://www.ncbi.nlm.nih.gov/gds/), The Gene Expression Omnibus Repository‌(GEO2R) tool was employed to pinpoint genes that are differentially expressed (DEGs). Using the GEIPA 2 website (http://gepia2.cancer-pku.cn/), we obtained differentially expressed genes of the The Cancer Genome Atlas - Pancreatic Adenocarcinoma (TCGA-PAAD) dataset. The GSE154778 dataset was used as the source of the single-cell dataset. Genes linked to NETosis were retrieved from the Gene Card (https://www.genecards.org/). The training set is the TCGA-PAAD queue, while the verification set is the external data set GSE57495.

Venn map was used to obtain NETosis-related DEGs. Univariate Cox analysis was employed to select NETosis-related genes linked to prognosis. The final gene signature was established using least absolute shrinkage and selection operator (LASSO) analysis [14].

The predict function determines each sample’s risk score based on the multivariate Cox model, and the median value is used to categorize the patients in the modeling and verification queue into groups with elevated risk and groups with reduced risk. To assess the disparities in the survival rates between the two cohorts, survival status was assessed using the Kaplan–Meier (KM) curve, and the receiver operating characteristic (ROC) curve was used to judge the prognostic model’s efficacy [15].

We investigated the association between infiltrating immune function and patient risk scores and utilized GSEABase and Gene Set Variation Analysis (GSVA) packages to determine the abundance of sixteen distinct varieties of infiltrating immune cells in both the high-risk and low-risk groups. The sensitivity of 196 medications was assessed for each patient risk group via the “oncoPredict” [16] package, a tool used for drug sensitivity prediction. The “calcPhenotype” function, which was trained using two datasets, GDSC2_Expr and GDSC2_Res, was used to forecast drug sensitivity based on patient data.

“Seurat”, an R package, was used for processing single-cell RNA sequencing (scRNA-seq) data, reducing dimensionality, performing normalization, and checking data quality. Extreme gene expression levels (nfeature RNA $>$200 and nfeature RNA 5000) were eliminated. The t-distributed stochastic neighbor embedding (t-SNE) approach was applied to standardize and decrease the dimensions and cluster the expression data. In accordance with the CellMarker database, the cells were labeled and categorized.

The human normal pancreatic duct epithelial cell line named HPDE6-C7 (STCC11108), along with the human pancreatic cancer cell lines PANC-1 (STCC11102) and SW1990 (FH0784), were procured from Servicebio and Fuheng biology. Vivacell was the supplier from which fetal calf serum and Dulbecco’s Modified Eagle’s Medium (C3113-0500) were procured. LGALS3 human pre-designed siRNA (HY-RS07604, siRNA sequences are presented in Table 1) and Cell Counting Kit-8 (HY-K0301) were purchased from MedChemExpress. Lipofectamine2000 (11668-019) was purchased from Invitrogen. The Total RNA Extraction Kit (DP430) was purchased from Tiangen. SweScript All-in-One RT SuperMix for qPCR (G3337) and 2xUniversal Blue SYBR Green qPCR Master (G3326) are sourced from Servicebio. PCR primers were synthesized by Qingke Biologics, Inc.

PANC-1 and SW1990 cells were cultured in DMEM + 10%FBS + 1%PS complete medium, 5% CO₂, and 37 ° constant temperature. PANC-1 and SW1990 cells that were in the log phase of growth were seeded into 6-well plates at a cell concentration of 1 $\times$ 10⁵. When the cell confluence was about 50%, they were replaced with a serum-free medium in which Lipofectamine2000 and Sh-LGALS3-siRNA plasmid were fully mixed at a ratio of 1:1. After the cells were treated for 6 hours, the serum-free medium was discarded and replaced with complete medium. The siRNA sequence with the highest interference effectiveness was chosen for further research after the total ribonucleic acid (RNA) was collected during 48 hours of continuous culture and the transfection efficiency was confirmed by reverse transcription quantitative real-time polymerase chain reaction (RT-qPCR). All cell lines were validated by short tandem repeat (STR) profiling and tested negative for mycoplasma.

Cells of the PANC-1 and SW1990 lines that were in the logarithmic growth stage were seeded into 96-well plates at a cell concentration of 5 $\times$ 10⁴ cells per milliliter, and each well was inoculated with 100 µL. After the cells adhere to the wall, each group of cells is provided with 5 multiple holes. Following a 48-hour treatment period, 10 µL Cell Counting Kit-8 (CCK-8) reagent was introduced into each well containing the complex, and then cultured in an incubator at 37 ° for 3 hours. After that, the 96-well plate was taken out, and an enzyme-labeled apparatus was employed to determine the value of the optical density (OD) at a wavelength of 450 nanometers and statistically analyzed. Cell survival rate = [(experimental hole absorbance-blank hole absorbance)/(control hole absorbance-blank hole absorbance)] $\times$ 100%.

The logarithmic growth period PANC-1 and SW1990 cell lines were digested and diluted to create 5 $\times$ 10⁶ cells per milliliter, which were then evenly distributed in the 6-well plates in 100 µL. At approximately 90% cell convergence, scratch with 200 µL gun head, cleaned with PBS, filled with culture medium, and imaged at 0, 24, and 48 hours, respectively.

The total RNA was isolated from HPDE6-C7, PANC-1, and SW1990 cells. Ultra-micro spectrophotometer (NanoDrop2000, Thermo Fisher Scientific, USA) was employed to measure the concentration of the total RNA. Subsequently, the total RNA was reverse-transcribed into complementary deoxyribonucleic acid (cDNA) using the PrimeScript™ RT reagent Kit with gDNA Eraser. Following the synthesis of cDNA, amplification was carried out, and the expression levels were measured with TB Green Premix Ex Taq II. The relative expression of the core genes was computed using the 2^{-Δ⁢Δ⁢Ct} method. The primer sequences are presented in Table 2.

The dataset was subjected to analysis using R project (version 4.4.2, Lucent Technologies, Murray Hill, New Jersey, USA), specifically version 4.3.1. For survival analysis, Kaplan–Meier survival analysis along with log-rank tests were employed. When dealing with variables, if a variable followed a normal distribution, the Student’s t-test was selected for analysis. In cases where the variable did not conform to a normal distribution, the Wilcoxon test was utilized instead. To build the ROC curve, the “timeROC” package was applied. Across all investigations, statistical significance was defined as a p-value less than 0.05.

Through a comparison of the gene expression patterns of 26 tumor tissues and 26 normal tissues within the GSE16515 dataset, 2682 DEGs were identified (log2[fold change] $>$1, Q 0.05). A total of 1969 genes were upregulated in the tumor group, whereas 711 genes were downregulated (Fig. 1A). The GEPIA database yielded a total of 9222 DEGs in PC (log2[fold change] $>$1, Q 0.05), of which 481 were downregulated and 8741 were upregulated in tumor tissues (Fig. 1B). We obtained 17 NETosis-related DEGs (Fig. 1C) by constructing Venn diagrams to further investigate the relationships between these DEGs and NETosis-related genes. Following an initial screening of genes associated with survival via univariate cox regression analysis, eleven genes (C5, CDC42, CXCL8, FN1, LGALS3, PLAUR, S100A11, DUOX2, SLC44A2, STAT1, and TNFSF10) were pinpointed that met the criterion of p 0.05 (Fig. 1D). LASSO regression analysis was conducted using the best $\lambda$ value construction of six marker gene features (FN1, STAT1, TNFSF10, SLC44A2, LGALS3, S100A11; Fig. 1E,F).

Fig. 1.

Identification of neutrophil extracellular trap related genes signature. (A) Volcano plot of differentially expressed genes (DEGs) in GSE16515 dataset. (B) DEGs in The Cancer Genome Atlas - Pancreatic Adenocarcinoma (TCGA-PAAD) dataset. (C) Venn of neutrophil extracellular trap (NETosis) related DEGs. (D) Result of univariate Cox regression analysis, and 11 genes showed p 0.05. (E) The least absolute shrinkage and selection operator (LASSO) regression of the 11 genes associated with overall survival (OS). (F) LASSO regression parameter selection is tuned using cross-validation.

A predictive gene model was established using the TCGA-PAAD dataset. The queue for the modeling was divided into groups with high and low risks based on the median of the risk score (Fig. 2A). A notable disparity in overall survival (OS) was evident between the low-risk and high-risk groups (p = 0.0159, Fig. 2B). Time-dependent ROC analysis was employed to evaluate the sensitivity and specificity of the prognostic model. The outcomes indicated that the area under the ROC curve (AUC) was 0.736 for 2-year survival, 0.792 for 3-year survival, and 0.872 for 4-year survival (Fig. 2C). Validation queues were derived from the GSE57495 dataset. Fig. 2D presents the dissemination of high- and low-risk classifications. Additionally, according to the Kaplan–Meier analysis, a significant difference in survival rates was detected between the high-risk and low-risk groups (p = 0.0263, Fig. 2E). The ROC curve for the validation dataset revealed that our model performed well in terms of prediction (0.698 for 2-year survival, 0.835 for 3-year survival, and 0.847 for 4-year survival; Fig. 2F).

Fig. 2.

Verification of the prognostic model associated with neutrophil extracellular trap. (A) Distribution of modeling queue based on the risk score. (B) The OS for modeling queue in the subunit is represented by Kaplan-Meier (KM) curve. (C) Receiver operating characteristic (ROC) curves of modeling queue. (D) Distribution of verification queue. (E) The OS for verification queue in the subunit is represented by KM curve. (F) ROC curves of verification queue. AUC, area under the ROC curve.

We employed single-sample immune infiltration analysis (ssGSEA) to contrast the concentration scores of sixteen distinct types of immune cells and the activity degrees of thirteen immune-associated pathways between the high- and low-risk cohorts. In the TCGA dataset, it was shown that the high-risk cohort generally exhibited a greater abundance of activated dendritic cells (aDCs), macrophages, Th1 cells, and Treg compared to the low-risk cohort (p 0.05, Fig. 3A). When conducting an analysis of the activities of immune-related pathways, it was found that in comparison to the low-risk cohort, the high-risk cohort exhibited a significantly higher level of enrichment for several pathways. These pathways include APC_co_inhibition, APC_co_stimulation, CCR, Check-point, HLA, Inflammation-promoting, MHC class I, Parainflammation, and Type IFN Response. The difference was statistically significant (p 0.05), as shown in Fig. 3B. We also examined the sensitivity of subgroups to frequently used chemotherapeutic medications for PC, namely, 5-fluorouracil, gemcitabine, and epirubicin. The group at high risk demonstrated a greater degree of sensitivity to all three pharmaceutical agents compared to the group at low risk. (p 0.01, Fig. 3C).

Fig. 3.

Immunoinfiltration and drug susceptibility of subgroups. (A) In the TCGA cohort, a comparison was made of the enrichment scores of 16 distinct types of immune cells between subgroups. (B) In the TCGA cohort, a comparison was carried out regarding the enrichment scores of 13 immune-associated pathways between subgroups. (C) The sensitivity of the 5-fluorouracil, gemcitabine, and epirubicin in various risk groups. ns: p $\mathrm{>}$ 0.05; *: p 0.05; **: p 0.01; ***: p 0.001; ****: p 0.0001.

Eight patients were included in the scRNA-seq dataset (Fig. 4A). A total of 49,056 genes and 6852 cells remained after removing some data with extremely low expression. Dimensionality reduction analysis was conducted using the t-SNE. Six cell clusters were identified (Fig. 4B). The bubble diagram’s results demonstrate that each cell cluster’s marker genes may be effectively distinguished (Fig. 4C). We used a t-SNE diagram to show the expression level of DEGs linked to NETosis in each cell cluster. The results demonstrated that epithelial cells exhibited high levels of LGALS3 and S100A11 (Fig. 4D). Since pancreatic cancer starts in epithelial cells, we speculate that these two genes have a crucial role in controlling the development and growth of tumor cells.

Fig. 4.

Single cell analysis results of pancreatic carcinoma. (A) Quality control results for eight patients. (B) Six clusters were identified based on the dimensionality reduction of t-distributed stochastic neighbor embedding (t-SNE). (C) Bubble diagram showing each cell cluster’s marker gene expression distribution. (D) The expression of NETosis-related DEGs in various cell clusters.

We verified the expression of LGALS3 and S100A11 in HPDE6-C7 and PANC-1 cells by RT-qPCR. The results demonstrated that pancreatic cancer cells had greater levels of LGALS3, and S100A11 expression than normal pancreatic duct epithelial cells (p 0.05, Fig. 5A). Findings from the survival analysis indicated that only pancreatic adenocarcinoma (PAAD) patients who exhibited elevated expression of the LGALS3 gene experienced poorer overall survival rates (p 0.05, Fig. 5B). This result emphasizes the significance of LGALS3 in pancreatic cancer, leading us to identify it as a critical gene and conduct an in-depth study into its mechanism. In order to explore the influence of LGALS3 on the functionality of pancreatic cancer cells, we inhibited the expression of LGALS3 in PANC-1 and SW1990 cells. After Sh-LGALS3-siRNA was transformed into PANC-1 and SW1990 cells, the total RNA was extracted and the interference efficiency was verified by RT-qPCR. The results demonstrated that the expression of LGALS3 was significantly inhibited by the Sh-LGALS3-siRNA sequence (p 0.001, Fig. 5C,D).

Fig. 5.

Verification of LGALS3 expression level. (A) Reverse Transcription Quantitative Real-Time Polymerase Chain Reaction verification of expression levels of LGALS3 and S100A11. (B) Survival analysis of S100A11 (p = 0.019) and LGALS3 (p = 0.079). (C) Results of interference efficiency of PANC-1 cells. (D) Results of interference efficiency of SW1990 cells. *: p 0.05; ***: p 0.001; ****: p 0.0001.

To further elucidate the mechanism by which LGALS3 modulates pancreatic cancer cells, we assessed the proliferative activity and migratory capacity of two types of pancreatic cancer cells subsequent to transfection with siRNA. The outcomes of the Cell Counting Kit-8 (CCK8) assay indicated that the proliferative capacity of PANC-1 and SW1990 cells in the Sh-LGALS3-1, Sh-LGALS3-2, and Sh-LGALS3-3 groups was notably reduced compared to that in the Sh-LGALS3-NC group (p 0.05, Fig. 6A,B). The results of the migration test demonstrated that interference with LGALS3 markedly reduced the migrating capacity of PANC-1 and SW1990 cells (p 0.05, Fig. 6C,D). Consequently, we verified that LGALS3 is capable of triggering the dissemination and metastasis of pancreatic cancer by modulating the proliferation and movement of cancer cells.

Fig. 6.

Interference with LGALS3 expression prevented PANC-1 and SW1990 cells from proliferating and migrating. (A) Cell Counting Kit-8 (CCK8) results of PANC-1 cells (p 0.05). (B) CCK8 results of SW1990 cells (p 0.05). (C) Results of migration of PANC-1 cells (p 0.05), scale length: 100 µm. (D) Results of migration of SW1990 cells (p 0.05), scale length: 100 µm. *: p 0.05; ***: p 0.001; ****: p 0.0001.