One scRNA-sequencing dataset (GSE117570) from the GEO database and two bulk-tissue RNA-sequencing datasets (TCGA_LUAD dataset from TCGA and GSE72094 dataset from GEO) were included in this study. Tumor bulk-tissue RNA-sequencing data and the corresponding clinical information of 594 patients with LUAD were obtained from TCGA database (https://portal.gdc.cancer.gov/) and used as the training cohort. The GSE72094 dataset containing RNA-sequencing and clinical data of 442 patients with LUAD were also collected as the validation dataset (25).

GSE117570 contains single-cell sequencing data of dissociated tumors or adjacent normal tissue from four NSCLC patients (three lung adenocarcinoma and one squamous cell lung cancer) (24). After exclusion of data from squamous cell lung cancer and its corresponding normal tissue, sequencing data from a total of 3,183 cells derived from three adenocarcinomas and their adjacent normal tissues were subjected to quality check and further study, which were all carried out by the Seurat package in R environment (R 3.5.1) (26). The percentage of mitochondrial genes was calculated by the PercentageFeatureSet function. All sequencing data had a reading depth of 10× genomics based on Illumina NextSeq 500. The relationship between sequencing depth and mitochondrial gene sequences, as well as total intracellular sequences, was evaluated by correlation analysis. HTseq counts of the remaining cells were normalized using a linear regression model with the LogNormalize method. During quality control, we further excluded a total of 2,078 low-quality cells with genes detected in <3 cells, <50 total detected genes, or ≥5% of mitochondria-expressed genes. The remaining cells were normalized using a linear regression model with the LogNormalize method. The top highly variant 1,500 genes were also identified by variance analysis.

Cell annotation by dimensionality reduction and trajectory analysis with pseudotime PCA was performed to identify significantly available dimensions with a p-value < 0.05 (27). Cluster classification analysis with dimensionality reduction for the top 15 principal components (PCs) was carried out using the t-distributed stochastic neighbor embedding (tSNE) algorithm to obtain the major clusters (28). Marker genes of each cluster were identified by differential expression analysis conducted using the limma package in R environment (29). Genes with an adjusted p value < 0.05 and | log2[fold change (FC)] | > 0.5 were considered to be significant. Annotation of cell clusters was carried out by the singleR package according to the composition patterns of the marker genes and were then manually verified and corrected with the CellMarker database (30, 31).

We then constructed the single-cell pseudotime trajectory using the Monocle 2 algorithm based on the scRNA-seq data derived from lung adenocarcinoma (32). We can interpret the differentiation trajectory of every single cell, as they were projected onto this space and ordered into a trajectory with branch points; cells within the same branch were considered to have a similar differentiation state, and vice versa. We also identified differentially expressed genes in cells with distinct differentiation states with | log2[fold change (FC)] | > 0.5 and FDR <0.05, which were defined as state-specific marker genes and were also defined as DRGs in our study. DRGs of different branches were applied to Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses to identify the enriched pathways. All pathway analyses were implemented by the “clusterProfiler” package in R environment. Pathways with an adjusted p value < 0.05 were considered to be significantly enriched.

To employ the DRGs identified from the trajectory analysis for classification of LUAD patients based on the bulk transcriptome data, we performed unsupervised consensus clustering on the GSE72094 dataset using the above-identified DRGs. Unsupervised consensus clustering was implemented by ConsensusClusterPlus’ package in R, using an algorithm based on k-means machine learning (33). The optimal number of clusters was decided based on the relative change in the area under the cumulative distribution function (CDF) curves of the consensus score and consensus heatmap. The clinical relevance of the clustering was evaluated in terms of its correlation with clinical features such as overall survival, age, sex, tumor stage, and oncogene mutation.

Tumor purity and the presence of infiltrating stromal/immune cells in tumor tissues were inferred based on the expression data using an ESTAMATE algorithm. Single sample gene set enrichment analysis (ssGSEA) and ESTAMATE signatures were integrated in the ESTAMATE algorithm, whose outputs include stromal score (represents the presence of stroma in tumor tissue), immune score (represents the infiltration of immune cells in tumor tissue), and estimate score (represents tumor purity) (34).

The composition of different immune cells was inferred by CIBERSORT, which is a web portal (http://cibersort.stanford.edu/) that performs cell type enrichment analysis for different immune cell types based on the bulk-tissue RNA-seq data (35). The relative infiltrating fractions of the 22 immune cell types in the tumor microenvironment were estimated by CIBERSORT.

Utilizing TCGA dataset, we conducted weighted gene co-expression network analysis (WGCNA) to explore the co-expression network of all the DRGs, which was implemented by the WGCNA package in R software. The distance between each DRG was calculated using the Pearson correlation coefficient, and a weighted co-expression network was constructed. Genes with missing values were identified and discarded, and a cluster tree was constructed to test the outliers. Network topology analysis was then applied to choose the soft-thresholding power. The expression matrix was transformed into an adjacency matrix and, subsequently, into a topological matrix (TOM). The hierarchical clustering tree of genes was shaped with the corresponding dissimilarity calculated. Based on the TOM, we used the average-linkage hierarchical clustering method to cluster genes. Modules with similar patterns through the Dynamic Tree Cut were identified and merged. We followed the standard of a mixed dynamic shear tree and set the minimum number of genes in each DRG network module to 20. Gene Significance (GS) and Module Membership (MM) were the quantifications of the relationship of genes to the trait and modules. Genes in the modules demonstrating the most remarkable association with clinical traits (such as age, sex, smoking history, stage, survival status, and survival time) were identified for further analysis. Genes with GS >0.5 and MM >0.5 were selected as candidate hub genes.

The candidate hub genes identified in WGCNA analysis were further applied to univariate survival analysis. Genes with a Univariate Cox value of p < 0.05 were further screened with the least absolute shrinkage and selection operator (Lasso) Cox model, which was implemented by package glmnet in R environment (36). Genes screened by Lasso Cox regression analysis were selected as hub genes and used for the construction of the risk prediction model, which was denoted as RiskScore in our study. Lambda derived from Lasso analysis was used as the coefficient of each gene in the calculation of RiskScore with the following formula: RiskScore =

where t represents the “t”th gene and n represents the total number of genes. The GSE72094 dataset was used as a testing cohort for external validation of RiskScore. The predictive power of the nomogram in the training cohort and validation cohort was evaluated and presented using receiver operating characteristic (ROC) curve analysis.

TCGA_LUAD dataset was applied as a training set for the development of a prognostic nomogram. Overall survival (OS) was defined as the time from the date of diagnosis or surgical resection to the date of death or most recent follow-up. Clinical parameters such as age, TNM stage, and RiskScore were included in univariate and multivariate Cox regression analyses using the backward stepwise Cox proportional hazard model. The coefficients for each parameter derived from multivariate analysis were used for the construction of a prognostic nomogram. The survival statuses of patients at 1, 3, and 5 years were used as the endpoint parameters for the development of the nomogram (37). The discrimination and calibration of the nomogram for the endpoint index were measured by the concordance index (C-index = 0.727) and by calibration plot comparing the expected and observed survival probabilities, respectively. The GSE72094 dataset was used as a testing cohort for external validation. The predictive power of the nomogram in the training cohort and validation cohort was evaluated and presented using the ROC curve.

Human lung cancer cell lines (A549) were purchased from the American Type Culture Collection (ATCC, USA). Cells were cultured with Dulbecco’s Modified Eagle’s Medium/Nutrient Mixture F-12 (DMEM/F12) (DMEM/F-12; Invitrogen, Carlsbad, CA, USA) supplemented with 10% fetal bovine serum (FBS; Gibco, Grand Island, NY, USA) and 4 mM of L-glutamine. Cells were cultured and incubated in a 5% CO₂ atmosphere at 37°C.

A549 cells (5 × 10⁴ cells per well) were seeded in 24-well plates and cultured overnight. Cells were transfected with 100 nM non-targeting siRNA, si-WFDC2 (Ruibo, Liaocheng, Shandong, China) using the Lipofectamine 3000 Transfection Reagent (Thermo Fisher, Waltham, MA, USA) following the manufacturer’s instructions. The negative control group was treated only with a transfection reagent. After 48 h, the transfection efficiency was confirmed by qRT-PCR.

The influence of WFDC2 on the proliferation of A549 cells was further evaluated using the EdU Apollo 567 In Vitro Kit (Solarbio, Beijing, China). Briefly, A549 cells transfected with or without si-WFDC2 were placed in 96-well plates at the density of 1 × 10³ cells. EdU solution (50 µM) was added to each well for 2 h. Then, cells were fixed with 4% paraformaldehyde for 10 min at room temperature. Cells were treated with glycine (2 mg/ml) and Apollo (1×), respectively. The nuclei were stained with 4′,6-diamidino-2-phenylindole (DAPI, Sigma, St. Louis, MO, USA). Positive staining was observed and captured under a fluorescence microscope (Olympus, Tokyo, Japan).

The proliferation of A549 cells transfected with or without si-WFDC2 cells were analyzed using the Cell Counting Kit-8 assay. Cells were seeded at 1 × 10³ cells per well into 96-well plates and cultured overnight at 37°C and 5% CO₂. After 24 h, CCK‐8 (Dojindo, Kumamoto, Japan) was added to each well and cultured for 2 h. The OD values at 450 nm were detected using the Multiskan Go Spectrophotometer (Thermo Fisher Scientific, USA) at days 1, 2, 3, and 4, respectively.

A549 cells transfected with or without si-WFDC2 (1 × 10³) were first blended into top agar (1.5 ml) for colony formation assay. Then the mixture was added onto base agar in each well. Three weeks post‐seeding, colonies were stained with 0.5% Crystal Violet for 15 min. To count the colonies, a single-lens reflex camera (Nikon) was used.

Cell apoptosis was evaluated utilizing FITC Annexin V Apoptosis Kit (BD Biosciences, San Jose, CA, USA) by flow cytometry assay. A549 cells (5 × 10⁵) transfected with or without si-WFDC2 were incubated with 1× binding buffer (100 μl) supplemented with PI (5 μl) and FITC Annexin V (5 μl) at room temperature and then analyzed using a flow cytometer (BD Biosciences).

The cell cycle was analyzed by a Cell Cycle Detection Kit (Keygen Biotech, Nanjing, China). Briefly, A549 cells transfected with or without si-WFDC2 were collected and fixed in 70% cold ethanol overnight. After that, cells were treated with 100 μl RNase A (100 mg/ml) for 30 min at 37°C and stained with 400 μl PI (50 mg/ml) at 4°C for 30 min. DNA content was measured using flow cytometry (Beckman Coulter FC500); the number of cells in the G1, S, and G2 phases were quantified by cell-cycle analysis software (FlowJo v10.0, Palo Alto, CA, USA).

The migration and invasion of A549 cells were evaluated using Transwell assays. A549 cells transfected with or without si-WFDC2 were pre-starved for 12 h, and 5 × 10⁴ cells were seeded in the upper chambers of the Transwell inserts coated with or without Matrigel^® (BD Biosciences, Bedford, MA, USA) solution, and the medium (DMEM supplemented with 10% FBS) was added into the lower chambers. After a 24-h incubation, the migrated and invaded cells were fixed with formalin for 10 min and stained using 0.1% crystal violet for 20 min at room temperature. Finally, the cell numbers were counted in five randomly selected fields under the microscope (Olympus, Japan).

Total RNAs were extracted using TRIzol^® Reagent (Invitrogen; Thermo Fisher Scientific, Inc.) according to the manufacturer’s instructions. The RNAs were transcribed into cDNAs with First Strand cDNA Synthesis Kit (Thermo Scientific, USA). Specific cDNAs were amplified with iTaq™ Universal SYBR^® Green (Bio-Rad, Hercules, CA, USA) utilizing Eco Real-Time PCR System (Illumina, San Diego, CA, USA). The results were analyzed using the 2−ΔΔCT relative quantitative method, with GAPDH as an internal control.

Between-group differences with respect to continuous variables were assessed using two-sample t test or Wilcoxon test. Between-group differences with respect to categorical variables were assessed using the χ² test, CMH-χ² test, or Fisher’s exact test, as appropriate. Pearson correlation analysis was conducted to assess the correlation between two continuous variables. Continuous variables were transferred to dichotomous variables with the median value as the cutoff point before inclusion in survival analysis. Log-rank test was used to evaluate the survival relevance of different parameters, which were presented with Kaplan–Meier plots. Logistic regression univariate analysis was conducted to explore the prognostic significance of each DRG. ROC curve analysis was performed to assess the performance of RiskScore in predicting the postsurgery survival probability of LUAD at 1, 3, and 5 years, respectively. All statistical analyses and data presentations were performed in R language 3.5.1 (http://www. r-project.org). p values less than 0.05 were considered indicative of statistical significance.

Single-cell sequencing data from dataset GSE117570 were analyzed. Standard quality control and normalization were performed for a total of 5,189 cells derived from the tumor mass of three cases of lung adenocarcinoma and 3,183 cells derived from adjacent normal tissues. A total of 2,078 low-quality cells were excluded from further analysis ( Figure 1A ). Significant positive correlations between sequencing depth and mitochondrial gene sequences (R = 0.08) and between sequencing depth and total intracellular sequences (R = 0.91) were detected ( Figure 1B ). A total of 3,162 genes detected in the dataset were included in the study. Variance analysis was performed for all the included genes, among which 1,500 genes were found to be highly variable (p < 0.05 and log(fold change) >1) across cells, while 1,662 genes showed no significant variation ( Figure 1C ). To classify the cells derived from tumor mass into different subsets, principal component analysis (PCA) was performed based on scRNA sequencing data. Although a total of 15 principal components (PCs) were identified ( Figure 1E ), PCA could not clearly separate the cells ( Figure 1D ). The top 20 highly expressed genes that were highly correlated with each of the first eight components are shown in Figure 1F .

To further identify cell clusters with distinct features, the t-distributed stochastic neighbor embedding (tSNE) algorithm was used to narrow down the cells into 10 distinct clusters. We also performed differential expression analysis to identify differentially expressed genes in different clusters, which were defined as the marker genes of the clusters. A total of 1,187 marker genes were identified; the top 20 marker genes for each cluster are presented in the heatmap ( Supplementary Figure 1 ).

The 10 clusters identified previously were annotated using singleR and CellMarker. Based on the expression pattern of the marker genes, the cellular subsets that were recognized during the annotation included T cells (cluster 0, containing 566 cells), monocytes (clusters 1 and 2, containing 334 cells and 265 cells, respectively), B cells (clusters 3 and 6, containing 130 cells and 86 cells, respectively), macrophages (cluster 4, containing 89 cells), epithelial cells (clusters 5 and 8, containing 87 cells and 51 cells, respectively), tissue stem cells (cluster 7, containing 57 cells), and endothelial cells (cluster 9, containing 30 cells) ( Figure 2A and Supplementary Table 6 ).

We further performed trajectory analysis to identify cells with a distinct differentiated pattern. All the cells derived from the tumor mass were projected into a three-dimensional graph with one root (II) and four branches (branches I, III, IV, and V) ( Figures 2B, C ). Branch-specific feature genes were identified as genes that were incrementally upregulated or downregulated from root to each branch ( Supplementary Tables 1–4 ). Feature genes that were upregulated in each branch were defined as differentiation-related genes (DRFGs) and were applied to GO and KEGG pathway analyses; all the significantly enriched pathways for different branches are shown in Supplementary Figures 2A–D . Pathways that were significantly enriched in each branch were mainly associated with immune regulation, such as immune cell infiltration or activation, antigen presentation, and cytokine production. The DRFGs in branches I/V were involved in neutrophil degranulation and activation, immune receptor activity, and inflammatory disease occurrence ( Supplementary Figures 2A, D ); those in branch II were also associated with a positive regulation of cytokine production and antigen processing and presentation ( Supplementary Figure 2B ), and those in branch III were closely linked to response to molecule of bacterial origin and Staphylococcus aureus infection pathways ( Supplementary Figure 2C ).

To classify LUAD patients based on the expression pattern of DRGs identified from scRNA sequencing analysis, machine learning-based unsupervised consensus clustering was carried out using the bulk transcriptome data of 442 LUAD tumors from the GSE72094 dataset. The number of clusters (k value) was determined based on the area under the CDF curve and consensus heatmap. When the k value was 4, the heatmap clearly classified the tumor into four clusters with distinct expression patterns ( Supplementary Figure 3A ), and the relative change in the area under the CDF curve was minimal with the increase in the k value ( Figures 3 A, B ). Thus, all LUAD patients from the GSE72094 dataset were classified into four groups, and their correlation with clinicopathological features was evaluated. The Kaplan–Meier survival plot was used to present the survival difference among patients of different groups; the results showed that patients of C3 had the worst overall survival compared with the other three groups (p = 0.003) ( Figure 3C ). The distribution of different clinicopathological features such as age, sex, smoking history, tumor stage, and mutation status of oncogenes or tumor-suppressor genes (KRAS, EGFR, TP53, STK11) in different groups was compared and presented as percentage bar plot ( Figures 3D–H , Supplementary Figures 3B–D ). The distribution of all these clinical features was comparable across the four groups except for the mutation status of TP53 and STK11. TP53 mutation occurred more frequently among patients of C3 ( Figure 3G ). The mutation of STK11 was significantly more frequent in C3 and C4 and extremely rare in C1 and C2 ( Figure 3H ).

We further evaluated the difference in the tumor microenvironment of the four different groups. As shown in Figure 4A , tumors from different groups significantly differed from one another in terms of immune score, stromal score, and tumor purity as inferred from the expression data. Of note, tumors from C3 demonstrated the lowest immune score and stromal score but the highest score for tumor purity. Also, tumors of the four groups differed from each other with respect to the composition of different immune cells ( Figure 4B , Supplementary Table 7 ). Specifically, tumors of C3 had a significantly lower level of adaptive immune population (such as B cells and T cells) but significantly more abundant innate immune cells (such as NK cells and macrophages) compared with other groups ( Supplementary Figure 3E ). We further evaluated the difference in immune checkpoint expression among tumors of the four groups. As shown in Figure 4C , the expression pattern of most of the known immune checkpoints was significantly different among the four groups.

Considering the significant difference in immune characteristics and overall survival among patients of the four groups, we anticipated that the immune microenvironment plays an important role in dictating tumor progression and clinical outcome. We further evaluated the prognostic relevance of different immune populations as well as immune-associated genes in LUAD patients using the GSE72094 dataset. As shown in Supplementary Figure 1 , lower levels of macrophage M0, activated mast cells, and plasma cells and higher levels of resting mast cells and CD4 memory resting cells were associated with significantly prolonged overall survival (p < 0.05). Immune-associated genes whose overexpression predicted favorable survival included CD40LG and PTPRC. Some other immune-associated genes such as YTHDF1, LDHA, PVR, and TNFSF4 correlated with unfavorable prognosis when highly expressed.

To better classify DRGs based on their co-expression pattern, WGCNA was performed for DRGs that are specific to each branch to construct co-expression patterns using TCGA-LUAD dataset. A total of four modules, MEturquoise, MEblue, MEbrown, and MEgrey, were identified, containing 326, 43, 29, and 14 DRGs, respectively ( Supplementary Table 5 ). The correlation between the expression of each module and clinical traits such as age, sex, smoking history, stage, survival status, and survival time was evaluated. As shown in Figure 5A , two modules (MEturquoise and Eblue) demonstrated the most remarkable association with clinical traits. Both MEturquoise and MEblue showed a significant correlation with deceased survival status, older age, no smoking history, and early tumor stage ( Figure 5A ). DRGs from modules MEturquoise and MEblue were selected as potential hub DRGs for the construction of the DRG-based prognostic model. All the DRGs from modules MEturquoise and MEblue were applied to differential expression analysis to identify their differential expression status between cells from the tumor mass and normal tissue ( Figure 5B , Supplementary Figure 5 ). A total of 15 genes were highly expressed in cells derived from tumor mass, and 85 genes were highly expressed in cells from normal tissue; the rest of the DRGs demonstrated no differential expression ( Figure 5B , Supplementary Figure 5 ).

Univariate Cox-regression survival analysis was performed to identify potential hub DRGs with significant prognostic relevance, which was further applied to multivariate Cox-regression survival analysis to evaluate their independent prognostic significance ( Figure 5B ). All survival analyses were based on TCGA-LUAD dataset. Only genes with independent prognostic significance were applied to the development of a prognostic model. The lambda value from Lasso analysis was carried out to construct the prognostic model.

A total of nine DRGs (TIMP1, VCAN, HLA-DPA1, HLA-DQA1, KRT8, MS4A7, SFTPB, WFDC2, FABP5) were finally recruited for the development of the prognostic model (Risk score) using TCGA-LUAD dataset. Pairwise comparison was carried out to evaluate the differential expression status of all the hub genes between tumor tissues and adjacent normal tissues using the expression data derived from 19 paired samples from our center (38). As shown in Supplementary Figure 2 , FABP5 and MS4A7 demonstrated a significantly higher expression in normal tissue (p < 0.05), while the expression levels of TIMP1, KRT8, and WFDC2 were significantly higher in tumor tissue (p < 0.05); the four remaining genes VCAN, HLA-DPA1, HLA-DQA1, and SFTPB demonstrated no significant differential expression between tumor and normal tissue (p > 0.05). The expression level of the hub genes in different cell types as annotated based on scRNA sequencing data was evaluated ( Figures 5C, D ). SFTPB, WFDC2, and KRT8 were highly expressed in epithelial cells; MS4A7, HLA-DPA1, and HLA-DQA1 were mainly expressed in immune cells such as macrophages and monocytes; and TIMP1 and VCAN were predominantly expressed by tissue stem cells ( Figures 5C, D ).

A prognostic model (RiskScore) integrating the expression level of the nine hub DRGs and their lambda value from Lasso analysis as coefficients was formulated as riskScore = SFTPB* -0.04820129 + WFDC2* -0.08632731 + HLA-DPA1* 0.22186232 + TIMP1* 0.26765595 + MS4A7*0.25567211 + HLA-DQA1*-0.15251071 + VCAN*0.12144478 + KRT8*0.26015183 + FABP5*0.22527365. Patients from TCGA cohort were categorized as high-risk group and low-risk group using the median value of RiskScore as the cutoff point ( Figure 6A ). The log-rank test showed that patients with low RiskScore had significantly better overall survival in TCGA-LUAD dataset (p < 0.001) ( Figure 6C ). To demonstrate the predicting power of RiskScore in an independent cohort, we evaluated the survival relevance of RiskScore in the GSE72094 dataset as a validation cohort. We obtained the risk score for each patients using the RiskScore model, based on which patients were categorized as high-risk group and low-risk group with the median value of risk score as the cutoff point. The high-risk group demonstrated significantly worse overall survival as compared to the low-risk group (p = 0.001) ( Figures 6B, D ). Time-dependent ROC curve analysis revealed impressive predictive power of the prognostic model in TCGA dataset [area under the curve (AUC) for prediction of 1-, 3-, and 5-year overall survival rates: 0.702, 0.651, and 0.632, respectively) ( Figure 6E ). The predictive ability of the prognostic model was also evaluated in an independent cohort (GSE72094), which showed comparable efficiency in predicting the overall survival rates at years 1, 3, and 5 years (AUC: 0.691, 0.610, and 0.679, respectively) ( Figure 6F ).

Univariate and multivariate Cox survival analyses were conducted to compare the prognostic value of RiskScore with that of the other clinicopathological variables such as age, sex, smoking history, and TNM stage ( Figures 6A, B ). Only TNM stage and RiskScore showed significant prognostic relevance in both univariate and multivariate survival analyses (p < 0.001) ( Figures 7A, B ); the hazard ratios (95% confidential intervals) in multivariate analysis were 1.54 (1.33–1.783) and 1.649 (1.466–1.855), respectively.

A prognostic nomogram integrating RiskScore and clinicopathological variables such as age and TNM stage was developed based on TCGA-LUAD dataset to better predict the prognosis of LUAD patients ( Figure 7C ). Survival statuses at 1, 3, and 5 years were applied as parameters of clinical outcome. The calibration plots showed excellent agreement between the OS predictions and the actual observations of the 0.5-, 1-, and 3-year survival rates in TCGA cohort ( Figure 7D ).

WFDC2, which was mainly expressed in epithelial cells and served as an immunity-related gene, was most significantly differentially expressed in our own RNA-seq database of 19 LUAD patients ( Supplementary Figure 2 , Supplementary Table 1 ). Therefore, WFDC2 was chosen for further verification in vitro. To clarify the functions of WFDC2 in lung cancer, the expressions of WFDC2 were silenced in A549 cells ( Figure 8A ). WFDC2 knockdown significantly reduced EdU-positive cells ( Figure 8B ). The results of CCK-8 assay also showed that the proliferation of A549 cells was decreased when WFDC2 was silenced ( Figure 8C ). Meanwhile, WFDC2 depletion repressed the colony numbers of A549 cells ( Figure 8D ). These results indicated that knockdown of WFDC2 inhibited the proliferation of A549 cells. Then we conducted flow cytometry to examine the effects of WFDC2 on the apoptosis of A549 cells. Annexin V-FITC and PI staining revealed that the number of apoptotic cells was obviously increased when WFDC2 was silenced in A549 cells ( Figure 8E ). For cell-cycle analysis, we observed that silence of WFDC2 significantly increased the fraction of cells in the S phase while it decreased the fraction of cells in the G1 phase, which indicated that the cell proliferation of A549 cells was reduced ( Figure 8F ). Finally, Transwell assay indicated that WFDC2 inhibition also decreased the numbers of migrated ( Figure 8G ) and invaded A549 ( Figure 8H ) cells which indicated that WFDC2 inhibition attenuated the migration and invasion of A549 cells (Fig.2C). All these results indicated that WFDC2 acted as a pro-oncogenic regulator in lung cancer.