To explore the EMT-related genes, which are associated with the oncogenesis and development of LUAD, transcriptome profiles of LUAD patients were obtained from The Cancer Genome Atlas (TCGA) database (https://portal.gdc.cancer.gov/). The “limma” package of R-language Bioconductor was used to identify the differentially expressed genes (DEGs) of LUAD patients under the conditions of false discovery rate <0.05 and |log(FC)| ≥ 4. The EMT-related genes were respectively extracted from the Genecard database (https://www.genecards.org/), the OMIM database (https://omim.org/), the EMTome (http://www.emtome.org/), the Epithelial-Mesenchymal Transition Gene Database (http://dbemt.bioinfo-minzhao.org/index.html), and the NCBI (https://www.ncbi.nlm.nih.gov/). Then the overlapping genes between the EMT-related genes and the LUAD associated genes (EMT-LUAD related genes) were identified with the Venn diagram and used for the following analysis (Wu et al., 2019; Li et al., 2020).

The ‘survival’ package in R-language was used to analyze the correlation between the EMT-LUAD related genes and the survival of LUAD patients.

First, to identify the EMT-LUAD related genes with significant prognostic value (p < .05), the univariate Cox proportional hazards regression analysis was used, and the correlation between the EMT-LUAD related genes and the survival of LUAD patients was analyzed.

Next, to identify EMT-LUAD related genes with independent significant prognostic value, a multivariate Cox proportional hazards regression model was constructed. The EMT-LUAD related genes with independent significant prognostic value was further analyzed (the mode of stepwise search is “both”). As a result, the risk scores for all EMT-LUAD related genes with independent significant prognostic value of each included individual were obtained [predict.coxph (type = “risk”) of “survival” package in R language], and the median risk score was used as the cut off value to divide the LUAD patients into two subgroups, including a high-risk score group and a low-risk score group.

Then, the expression heatmap of EMT-LUAD related genes with independent significant prognostic value in the high- and low-risk groups was plotted in R language using the package “pheatmap”. Meanwhile, the risk score curves and the patient survival scatter plots of all included individuals, the survival curves of the high- and low-risk groups, and the time-dependent ROC curves of this model were carried out in R-language with the package “survival” and “survminer”.

Last, the association between the EMT-LUAD related genes with independent significant prognostic value and the clinical characteristics were analyzed in R-language with the package “survival”, including the [age (≤60 versus >60), gender (male versus female), stage (stage I and II versus stage III and IV), T (T1 and 2 versus T3 and 4) and N(N0 versus N1–3)] (Fisher and Lin, 1999).

The TIMER, CIBERSORT, CIBERSORT-ABS, and QUANTISEQ algorithms were used to explore the cellular immune response between the low-risk group and the high-risk group mentioned above.

Meanwhile, the differences of the cellular immune response between different algorithms were visualized with the Heatmap.

Furthermore, ssGSEA was used to quantify the subgroups of the tumor-infiltrating immune cells, the immune check-points, as well as the immune functions between the low-risk group and the high-risk group (Tang et al., 2021).

The expression differences of m⁶A-related genes, including ALKBH5, HNRNPC, METTL14, METTL3, YTHDF2, IGF2BP1, IGF2BP3, METTL16, RBMX, IGF2BP2, and HNRNPA2BA, were analyzed between the low-risk and the high-risk groups mentioned above in LUAD patients (Huang et al., 2020).

Drug targets of the active ingredients in Chonglou (dioscin (Tao et al., 2018; Ding et al., 2020; Wu et al., 2020), gracillin (Liu et al., 2021; Yang J. et al., 2021; Yang L. et al., 2021), polyphyllin I (Chang et al., 2017; Tian et al., 2020; Lai et al., 2021; Shen et al., 2021), polyphyllin II (Niu et al., 2020), polyphyllin III (Zhou et al., 2021), polyphyllin VI (Liu et al., 2020; Teng et al., 2020), polyphyllin VII (He et al., 2020; Zhao et al., 2021) and polyphyllin H) were collected by using the online tools (Su et al., 2019), including PubChem (https://pubchem.ncbi.nlm.nih.gov/), Swiss Target Prediction (http://www.swisstargetprediction.ch/index.php) and TargetNet (http://targetnet.scbdd.com/home/index/).

The overlapping gene set between the drug targets of Chonglou active ingredients with the EMT-LUAD related genes were enquired in the STRING database (version 11.5) (EMT-LUAD-Chonglou related genes), which was used for the following protein-protein interactions (PPI) network, hub genes and enrichment analyses.

The protein-protein interaction (PPI) network map of the EMT-LUAD-Chonglou related genes was further obtained by using the online tool STRING database (version 11.5). Cytoscape software was then applied to analyze the degree values, and the top five genes were further identified as hub genes with the cytoHubba plugin (Li et al., 2021).

To explore the gene biological process of the EMT-LUAD-Chonglou related genes identified above, the package ggplot2 in R software was used to perform the gene ontology (GO), and the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses (Li et al., 2019; Liang et al., 2020).

Furthermore, a graphical representation of the drug-disease-target-GO function-pathway was constructed by using the Cytoscape software to create the GO and KEGG pathway analyses (Xu et al., 2021).

To carry out the molecular docking between the active ingredients of Chonglou with the five hub genes, the protein structures were obtained from the PDB database (https://www.rcsb.org/), and the compound structures were downloaded from the PubChem database (https://pubchem.ncbi.nlm.nih.gov/).

After setting the key parameters, such as the corresponding protein, the hydrogenation and the Gasteiger charge for merging nonpolar hydrogen atoms (Trott and Olson, 2010), the AutoDock Tools (version 1.5.7) were then used for the molecular docking to obtain the groups presenting the lowest protein binding energy with a rigid macromolecule as the protein structure and a genetic algorithm as the algorithm.

PyMOL (version 2.5) software was used to visualize the docking results.

The PC-9 cell lines present in this study were obtained from Nation Collection of Authenticated Cell Cultures. All cells were cultured in a 5% CO2 humidified incubator with a growth medium including FBS (10%), streptomycin (100 μg/ml) and penicillin (100 U/ml). Human lung adenocarcinoma cells were cultured at 37°C with DMEM.

After reviewing the literature we identified Dioscin and polyphyllin III as the same compound with the same CAS number, but with different spatial structures. Finally, We determined the drug concentrations of these seven active ingredients of Chonglou. Among them, the concentration of polyphyllin I, Polyphyllin II, Polyphyllin VII and Polyphyllin H was 1 μM, the concentration of Gracillin was 2 μM, the concentration of Polyphyllin VI was 32 μM, and the concentration of Polyphyllin III was 4 μM.

PC-9 cells (4×10⁴/100μl/well) plated in 96-well plate were placed in 96-well plates and treated with 2 μM of Gracillin, 1 μM of Polyphyllin I, 1 μM of Polyphyllin II, 4 μM of Polyphyllin III, 32 μM of Polyphyllin VI, 1 μM of Polyphyllin VII and 1 μM of Polyphyllin H. Cell viability was then measured using the Cell Counting Kit-8 (CCK-8) to determine cell viability at 24, 48, and 72 h.

Cell viability of human lung adenocarcinoma (PC-9) cells (3×10³/100 μl/well) plated in specific xCELLigence-E plates was sequentially and in real time monitored for 72 h using the xCELLigence-RealTime Cell Analyzer (Roche Applied Science, Mannheim, Germany).

PC-9 cells were cultured for 24 h in the upper chamber of a 24-well plate (3×10⁴/well/200 µl of serum-free culture medium) and the lower chamber was chemotactic culture medium (10% FBS, 300 µl). Cells were migrated to the other side of the chamber, washed lightly twice with PBS, fixed in 4% paraformaldehyde for 30 min, coloured with 1% crystal violet solution for 10 min as well as photographed under the microscope.

The xCELLigence RTCA DP system was applied to assess PC-9 cell invasion and migration. Briefly, for the specific CIM-Plate in the xCELLigence RTCA DP system, the lower chamber contained medium with 10% FBS (165 μl/well). 40 μl/serum-free medium were used in the upper chamber to determine background. Cells in serum-free medium (3 × 10⁴/100 μl/well) were then placed in the upper chamber. Cell migration capacity was monitored continuously for 24 h. For the cell invasion capacity assay, with the exception of applying the upper chamber pre-coated with Matrigel and monitored continuously for 48 h, the same conditions as for the migration assay were used.

After a comprehensive search from the public database, 4920 EMT-related genes were collected from Genecard, OMIM, and EMTome databases (Figure 1A, left circle), meanwhile, 804 DEGs between the LUAD and the normal control tissues were identified from TCGA database under the conditions of false discovery rate <0.05 and |log(FC)| ≥ 4 (Figure 1A, right circle). A total of 166 overlapping genes were recognized as the EMT-LUAD related genes (Figure1A, overlapping part), including 7 down-regulated and 159 up-regulated genes in the LUAD tissues versus the normal control tissues from LUAD patients (Figure 1B).

As a statistical method applicable to exploring the prognostic value of a single variable, the univariate Cox proportional hazards regression analysis was used to identify the potential association between the 166 EMT-LUAD related genes obtained above with the overall survival (OS) of LUAD patients. A total of 16 genes (AGER, DEPDC1, FAM83A, HMGA2, HOXA13, KRT16, SERPINB5, S100P, TUBB3, IGF2BP1, IGF2BP3, NTSR1, KRT6A, FAM83A-AS1, HJURP, and TMPRSS11E) were identified to be significantly associated with the OS of LUAD patients (Table 1 and Figure 2, p < 0.05), indicating the important prognostic values of these genes in LUAD patients.

Since the existing limitation of univariate Cox proportional hazards regression analysis, it is not applicable to identify whether a particular factor is an independent prognostic factor. Therefore, the multivariate Cox proportional hazards regression model was constructed, including correlation coefficient and hazard ratio.

As a descriptive statistic, the correlation coefficient is a bivariate statistic to summarize the association between two variables, which uses a number between −1 and 1 to represent the strength and direction of an association between variables, such as to describe how a specific predictor factor affects the survival of patients. The absolute value reflects the correlation strength, and the positive or negative value reflect a positive or negative correlation (Schober et al., 2018). The hazard ratio reflects how a predictor factor affects the survival (Trkulja and Hrabač, 2019). For each of the sixteen EMT-LUAD related genes with important prognostic value in LUAD patients identified by the above univariate Cox proportional hazards regression analysis, the correlation coefficient and the hazard ratio was calculated. Eight genes, including IGF2BP1, HOXA13, NTSR1, KRT6A, TMPRSS11E, AGER, FAM83A, and FAM83A-AS1, were identified to be independently and significantly associated with the prognosis of LUAD patients (Table 2). Among them, AGER and FAM83A-AS1 were favorable prognostic factors (hazard ratio <1), and others (FAM83A, HOXA13, IGF2BP1, NTSR1, KRT6A, and TMPRSS11E) were harmful prognostic factors (hazard ratio >1).

To further investigate the association between the expression levels of these eight genes with the prognosis of LUAD patients, the risk scores for these eight genes were obtained to calculate the median risk score, which was used as the cut off value to divide the 477 LUAD patients into a high-risk score group (n = 237) and a low-risk score group (n = 240). Then different clinical characteristics of LUAD patients were analyzed based on the high- and low-risk scores.

The expression differences of the eight EMT-LUAD related genes with independent significant prognostic value between the high-group and the low-risk group were visualized by a Heatmap (Figure 3A).

The curve of risk score and the survival status of the patients were displayed (Figures 3A–C), which showed a difference between the high-risk group and the low-risk group in the survival rate, and patients in the high-risk group had a shorter overall survival than patients in the low-risk group (about 25% vs. 50% of the 5-year survival rate, p < 0.001) (Figure 3D). The AUC values of the ROC curves for the 1, 3, and 5-year overall survival were .729, .7, and .679, respectively, indicating a favorable predictive accuracy of the prognostic signature (Figure 3E).

Furthermore, the clinical characteristic analysis of the eight genes showed that the expression levels of FAM83A, FAM83A-AS1, and TMPRSS11E were significantly related (p < .05) to the number and scope of lymph node metastasis and the spread in the lymph node area, which showed that theses the expression of these three genes in the N1-3 status is higher than those in the N0 status. Moreover, FAM83A and FAM83A-AS1 were up-regulated in the stage III-IV patients than in the stage I-II patients (Table 3 and Figure 4), while AGER, HOXA13, IGF2BP1, IGF2BP1, and KRT6A were not significantly related to the development of LUAD.)

Since the immune system plays an important role in the oncogenesis and development of LUAD (Anichini et al., 2020), it is helpful for the clinical practice to explore the immunological regulatory effect of the EMT-LUAD related genes with independent prognostic values. Therefore, the correlation analysis between the immune cell infiltration with the risk scores of the eight EMT-LUAD related genes with independent prognostic values in the LUAD patients were further performed. And the immune responses calculated by using the TIMER, CIBERSORT, CIBERSORT-ABS, and QUANTISEQ algorithms between the high-risk group and the low-risk group were visualized with Heatmap, which revealed the correlation between the immune cell infiltration and the EMT-LUAD related genes with independent prognostic values. As shown in Figure 5A, immune cells were identified to be characteristically enriched in the low group or the high-risk groups, such as B cells calculated by TIMER, CIBERSORT and CIBERSORT-ABS, monocytes calculated by CIBERSORT and CIBERSORT-ABS, mast cells activated calculated by CIBERSORT and CIBERSORT-ABS and T cell CD4⁺ memory resting calculated by CIBERSORT-ABS being significantly enriched in the low-risk group, while neutrophils calculated by CIBERSORT and QUANTISEQ, T cell CD4⁺ memory activated calculated by CIBERSORT and CIBERSORT-ABS, mast cell resting calculated by CIBERSORT and CIBERSORT-ABS, macrophage M0 calculated by CIBERSORT and T cell CD4⁺ (non-regulatory) calculated by QUANTISEQ being significantly enriched in the high-risk group, indicating the close association between the immune infiltrations and the risk scores of the eight EMT-LUAD related genes with independent prognostic values in the LUAD patients.

Furthermore, the correlation analysis between the immune cell subsets/functions and the risk scores of the eight EMT-LUAD related genes with independent prognostic values in the LUAD patient was performed, which showed that APC_co_stimulation, T_cell_co-stimulation, and type II INF response was significantly lower in the high-than in the low-risk group (Figure 5B).

Since the checkpoint inhibitor is one of the most applicable immunotherapy strategy (Waldman et al., 2020), to identify the association between the immune checkpoints with the risks of the LUAD patients has important implications for the development of immunotherapy. Our results displayed a significant difference in the expression of CD274, CD176, LGALS9, TNFSF15, TNFSF14, and TNFSF25 between the high- and the low-risk groups. Among them, the CD274 and CD276 were downregulated, while the LGALS9, TNFSF15, TNFSF14, and TNFSF25 were upregulated in the high-risk group than in the low-risk group (Figure 5C).

M⁶A-related genes are important regulators in cancer development, and affect the pathological process of various cancers (Huang et al., 2020). Further comparative analysis of m⁶A-related genes between the high-risk group and the low-risk group mentioned above revealed significant differences in HNRNPC, METTL14, METTL3, YTHDF2, IGF2BP1, IGF2BP2, IGF2BP3, METTL16, and HNRNPA2B1. Among them, the expressions of HNRNPC, IGF2BP1, IGF2BP3, IGF2BP2, and HNRNPA2B1 were upregulated in the high-risk group than in the low-risk group, while the others (METTL14, METTL3, YTHDF2, and METTL16) were downregulated in the high-risk group than in the low-risk group (Figure 6).

Since the EMT-related genes play important roles in LUAD and the active ingredient of Chonglou was widely used in the treatment of various cancers, it is necessary to explore the potential drug targets of the Chonglou active ingredients from the EMT-LUAD related genes (EMT-LUAD-Chonglou related genes). Therefore, Swiss Target Prediction and TargetNet database were used for data extraction and 561 drug targets of Chonglou active ingredients were obtained (Figure 7A, left circle). Then the EMT-LUAD-Chonglou related genes were further identified by Venn diagram from the 561 drug targets of Chonglou active ingredients and the 166 EMT-LUAD related genes (Figure 7A, right circle), and 12 EMT-LUAD-Chonglou related genes were identified (Figure 7A, the overlapping part).

PPI network is used to identify the proteins based on their interaction with each other to participate in the life processes such as the biological signaling, the gene expression regulation, the energy and substance metabolism, and the cell cycle regulation (Cafarelli et al., 2017). By STRING analysis, we identified that the PPI network consisting of 12 intersection targets of EMT-LUAD-Chonglou related genes (Figure 7B). Then, the Cytoscape was further used to identify the five core targets by degree value, and F2, MMP1, MMP3, MMP10, and MMP13 were identified (Figure 7C), which suggests that these five core targets play crucial roles in the LUAD development by regulation of other genes.

In order to explore the biological pathways playing important roles in Chonglou against LUAD, the 12 intersective genes obtained by STRING analysis mentioned above were submitted for GO and KEGG enrichment analysis. The results indicated that Chonglou active ingredients mainly affected the BPs, including extracellular matrix disassembly, collagen metabolic process, collagen catabolic process, extracellular matrix organization, extracellular structure organization, and response to amyloid-beta. As for the MFs, Chonglou active ingredients mainly affected the endopeptidase activity, metalloendopeptidase activity, metallopeptidase activity, serine-type endopeptidase activity, serine-type peptidase activity, and serine hydrolase activity and neuropeptide receptor activity. Furthermore, in the KEGG pathway analysis four pathways were identified, including neuroactive ligand-receptor interaction, IL-7 signaling pathway, coronavirus disease-COVID-19, and rheumatoid arthritis (Figure 7E). Moreover, the bioinformatics and network pharmacology were constructed, the results highlighted that extracellular matrix and anti-inflammatory were the key target/pathway of EMT-LUAD-Chonglou related genes (Figure 8).

Molecular docking models can predict the binding mode and binding affinity between the molecules to account for the differences in compound activities. To explore the potential bindings between the Chonglou active ingredients and the five core targets mentioned above, the molecular docking analysis was performed. The results revealed that the Chonglou active ingredients had strong binding abilities with the core targets, which showed that the binding energies between them were all less than −5 kcal mol⁻¹ (Table 4). The docking results were visualized and the positions of the small molecules to form the hydrogen bonds with amino acid residues were marked on the map (Figure 9). Taken together, these bond interactions indicated that the Chonglou active ingredients had strong binding activities with the five core proteins.

The active ingredient of Chonglou inhibited the cell viability and cell proliferation of lung adenocarcinoma cells. In the CCK-8 assay, the active ingredient of Chonglou significantly reduced the viability of lung adenocarcinoma cells intime-dependent manner, and subsequent observation of cell morphology under crystalline sub-staining revealed that the active ingredient of Chonglou significantly reduced the number of cell compared to untreated control cells. Furthermore, the viability of PC-9 cells was measured by xCELLigence RTCA DP system significantly inhibited the cell viability of lung adenocarcinoma cells in a time-dependent manner (Figure 10).

The Transwell migration examination and Real time migration and invasion assay were performed to investigate whether the eight active components of Chonglou affect the metastasis of human lung adenocarcinoma cells. Firstly, it was found that Gracillin, polyphyllin I, polyphyllin II, polyphyllin Ⅲ, polyphyllin Ⅵ, polyphyllin Ⅶ and polyphyllin H inhibited the migration and invasion of tumor cells compared with the control group, and it was found that the inhibition of PC-9 cell migration ability was time-dependent (0–24 h) (Figure 11).