Skip to main content

ORIGINAL RESEARCH article

Front. Cell Dev. Biol., 20 January 2022
Sec. Molecular and Cellular Oncology
This article is part of the Research Topic Single Cell Meets Metabolism and Cancer Biology View all 10 articles

Single-Cell RNA-Seq Reveals the Promoting Role of Ferroptosis Tendency During Lung Adenocarcinoma EMT Progression

Jiaxi Yao,&#x;Jiaxi Yao1,2Yuchong Zhang&#x;Yuchong Zhang1Mengling LiMengling Li3Zuyu SunZuyu Sun2Tao Liu
Tao Liu2*Mingfang Zhao
Mingfang Zhao1*Zhi Li
Zhi Li1*
  • 1Department of Medical Oncology, The First Hospital of China Medical University, Shenyang, China
  • 2Department of Urology, The First Hospital of China Medical University, Shenyang, China
  • 3Department of Clinical Epidemiology and Center of Evidence-Based Medicine, The First Hospital of China Medical University, Shenyang, China

Epithelial-mesenchymal transition (EMT) and ferroptosis are two important processes in biology. In tumor cells, they are intimately linked. We used single-cell RNA sequencing to investigate the regulatory connection between EMT and ferroptosis tendency in LUAD epithelial cells. We used Seurat to construct the expression matrix using the GEO dataset GSE131907 and extract epithelial cells. We found a positive correlation between the trends of EMT and ferroptosis tendency. Then we used SCENIC to analyze differentially activated transcription factors and constructed a molecular regulatory directed network by causal inference. Some ferroptosis markers (GPX4, SCP2, CAV1) were found to have strong regulatory effects on EMT. Cell communication networks were constructed by iTALK and implied that Ferro_High_EMT_High cells have a higher expression of SDC1, SDC4, and activation of LGALS9-HARVCR2 pathways. By deconvolution of bulk sequencing, the results of CIBERSORTx showed that the co-occurrence of ferroptosis tendency and EMT may lead to tumor metastasis and non-response to immunotherapy. Our findings showed there is a strong correlation between ferroptosis tendency and EMT. Ferroptosis may have a promotive effect on EMT. High propensities of ferroptosis and EMT may lead to poor prognosis and non-response to immunotherapy.

Introduction

Lung cancer is the most commonly diagnosed cancer and the leading cause of cancer death, which includes small cell lung cancer (SCLC) and non-small cell lung cancer (NSCLC). Lung adenocarcinoma (LUAD) is the most prevalent NSCLC. Currently, the main treatment for advanced LUAD including targeted therapy, radiotherapy, chemotherapy and immunotherapy. Despite the application of new therapies in recent years, advanced LUAD remains prone to metastasis, resulting in a poor prognosis with a 5-year survival rate of less than 20%. Therefore, further research on the mechanisms of lung cancer metastasis is urgently needed.

One of the most familiar factors in the metastasis of tumors is epithelial-mesenchymal transition (EMT). During EMT, epithelial cells lose their original polarity and intercellular adhesion and acquire invasive properties associated with the mesenchyme (Chen et al., 2021). Furthermore, it can contribute to tumor progression and has been recognized as a driver of tumor spread in many types of tumors (Puram et al., 2017). Current research suggests that EMT transcription factors like SNAI1, TWIST1 and ZEB1 may contribute to the promotion of tumor metastasis and drug resistance by EMT (Chen et al., 2021). EMT is associated with resistance to chemotherapy and immunotherapy, it is therefore considered to play an important role in the treatment of tumors (Panchy et al., 2019; Bulle and Lim, 2020). It has been previously published that in NSCLC, EMT may lead to exhaustion of CD8+ T cells, resulting in immunosuppression (Soundararajan et al., 2019). EMT also interacts with a variety of biological processes, leading to metastasis and poor prognosis. As a result, extensive research has been carried out on the relationship between EMT and other biological processes in recent years.

Ferroptosis is a newly identified iron-dependent form of programmed cell death caused by unrestricted lipid peroxidation and plasma membrane rupture (Dixon et al., 2012). The accumulation of iron, production of free radical and toxic lipid peroxides are the main features of ferroptosis, which was verified to be related to a variety of molecular mechanisms (Wu et al., 2021a). Cells underwent ferroptosis may bring an abundance of oxidized lipids, which can affect the phagocytosis of dendritic cells and the presentation of antigens, and then cause immune evasion of tumor cells (Demuynck et al., 2021). Previous research has shown EMT can increase cellular sensitivity to ferroptosis. On the one hand, it has been known that a selective vulnerability to ferroptosis is related to highly mesenchymal-like cell state (Viswanathan et al., 2017). Protein LYRIC can promote EMT, while at the same time, promote ferroptosis by inhibiting GPX4 (Bi et al., 2019). On the other hand, ferroptosis tendency also contributes to EMT (Du et al., 2021). Erastin, which is an inducer of ferroptosis, can act on lung epithelial cells and causes de-epithelialize (Sun et al., 2021). PTGS2 is another ferroptosis promotor, the silencing of which can inhibit the proliferation and migration of pancreatic cancer cells (Xu et al., 2021a). Combined with the conclusions in the previous literatures, EMT and ferroptosis tendency may form a positive feedback loop within a certain range. However, the existing studies were conducted based on cell line experiments that lack of complex tumor microenvironment, which may lead to results that are too artificial to address clinical questions. In vivo environment, the relationship between ferroptosis tendency and EMT has been little studied and is not yet clear. Hence, there is an urgent need for a more precise analysis.

Single-cell RNA sequencing (scRNA-seq) can partly solve the problems above. scRNA-seq uses samples obtained from human body for analysis and can represent the results of in-vivo experiments. It is a high-resolution tool that overcomes the limitations of traditional bulk sequencing, allowing us to measure tumor heterogeneity at the level of individual cell and is valuable for the study of the tumor microenvironment (Xu et al., 2021b; Liu et al., 2021). Previous studies of vitro experiments have shown that ferroptosis is involved in the EMT process, but the exact regulatory mechanism is unclear in vivo. Therefore, we used scRNA-seq to assess the relationship between ferroptosis tendency and EMT in LUAD.

In data analysis after scRNA-seq, we get many correlation conclusions. Various statistical models, machine learning, and deep learning models are used to get various correlations through analysis, but correlation cannot represent the objectivity and authenticity of causality. Correlation indicates two variables are related, which show an increase or decrease trend (Altman and Krzywinski, 2015). Causality is where the cause is partly responsible for the effect, and the effect is partly dependent on the cause. The main difference between causal inference and correlation inference is that the former analyzes the response of the effect variable when the cause is changed (Pearl, 2009; Vemuri, 2015). In addition, at present, regulation is only verified at single molecular level, probably because for basic experiments, the overall regulatory network is difficult to construct. To clarify the dependence of different biological processes, causal inference is increasingly used in biological data such as genome, proteome, and scRNA-seq data to explain the regulatory network between molecules (Sella et al., 2017).

Here, we used scRNA-seq data analysis and found that ferroptosis tendency in tumor cells is closely related to EMT. At the same time, we used SCENIC and causal inference to construct a directed network for key TFs and markers, and found that the tendency of ferroptosis has a positive regulatory effect on EMT, which plays an important role in tumor proliferation, metastasis and ineffective immunotherapy. This result will provide important support for the clinical treatment decision of LUAD.

Results

EMT and Ferroptosis in mLN are Positively Correlated

Processed scRNA-seq data were screened with strict QC criteria in the origin article, and we use the original meta. data file to annotate all cells (Kim et al., 2020). 36,071 ECs from different multiple histologic sections were obtained for our follow-up analysis (Figure 1A). After performing gene filtering, normalization, principle component analysis, we first set the resolution at 0.7, resulting in the 37 cell clusters, and visualized the original samples utilizing t-SNE plots (Figures 1B,C).

FIGURE 1
www.frontiersin.org

FIGURE 1. Overall flowchart and t-SNE dimensionality reduction clustering. (A) Flowchart of the work. (B,C) t-SNE plots showed the clustering results and the different origins of epithelial cells, including nLung, tLung, tL/B, mBrain, mLN.

EMT is a cell reprogramming process, which ECs acquire a mesenchymal phenotype. A series of changes activated after EMT help tumor cells to invade surrounding tissues and metastasize (Dongre and Weinberg, 2019; Pastushenko and Blanpain, 2019). Since the tumor tissues come from different metastatic sites, we put the research center on EMT, using 23 EMT up-regulated and 3 EMT down-regulated gene sets obtained from MSigDB to perform ssGSEA on all ECs (Subramanian et al., 2005). The analysis revealed that the expression of EMT-related gene set was down-regulated in normal tissues (green box), while globally up-regulated in tumor derived tissues, especially in mLN (purple box) (Figure 2A). This result suggests that there are some EMT-generating cells in the tumor tissue, which may cause the primary tumor metastasis of lung cancer. Considering EMT was strongly associated with ferroptosis tendency, we further explored various forms of cellular death including ferroptosis using the ssGSEA method (Figure 2B). Compared with other forms of death, we found that ferroptosis has specific activation in different tissues. In mLN, the ferroptosis pathway in most ECs is activated, which was the same as the activation of the EMT pathway (Figure 2B). It is worth noting that some cells in mBrain and other sample tissues are also highly sensitive to ferroptosis. However, unlike mLN, in other tissues, only a part of cells exhibited a high sensitivity to ferroptosis. This may have some relation to the tumor microenvironment. Additionally, the heterogeneity in other tissues such as mBrain may also be a reason for this inconsistency. To further confirm the relationship between ferroptosis tendency and EMT, we divided all ECs into four groups based on the median score: Ferro_High_EMT_High, Ferro_High_EMT_Low, Ferro_Low_EMT_High, Ferro_Low_EMT_Low (Figure 2C). Through the pie chart, we can clearly see that the proportion of “Ferro_High_EMT_High” cells (double-high group) in mLN was as high as 49%, which suggested that there was a strong correlation between EMT and ferroptosis in mLN.

FIGURE 2
www.frontiersin.org

FIGURE 2. EMT and ferroptosis were significantly active in mLN. (A) Heatmap showed the expression of EMT-related genes in cells of different origins. ssGSEA score of EMT was low in normal tissues (green box) and high in mLN (purple box). (B) t-SNE plot showed the score of various forms of cellular death. Box plot showed ferroptosis tendency scores of cells from different tissue sources. mLN had the highest score (red box). (C) All ECs was divided into four groups based on the median score: Ferro_High_EMT_High, Ferro_High_EMT_Low, Ferro_Low_EMT_High, Ferro_Low_EMT_Low and are summarized in the pie chart.

Marker Expression and Enrichment Analysis Suggested a Strong Correlation Between Ferroptosis Tendency and EMT in all Epithelial Cells

Among the proportions of cells from the four original samples, we found that the sum of the proportions of the double-high group and the double-low group accounted for about 70% (Figure 2C). This result proved that ferroptosis and EMT might occur simultaneously in ECs. In order to prove that this result was not unique in mLN, we divided all ECs into “high” and “low” groups according to the median ssGSEA score of ferroptosis, and all ECs were grouped and analyzed according to original samples (Figure 3A). We found that the EMT score was significantly increased in the group with high ferroptosis group, no matter in the primary tumor (tLung, tL/B) or the metastatic site (mLN, mBrain). This result proves that EMT displays homotropic cooperativity with ferroptosis tendency are not only in metastatic lymph nodes, but are closely related in all lung cancer cells, which is a common phenomenon.

FIGURE 3
www.frontiersin.org

FIGURE 3. In tLung, tL/B, mLN, mBrain, EMT all exhibited homotropic changes with ferroptosis. (A) Heatmap showed the EMT score was significantly increased in the group with high ferroptosis tendency score group. (B) Dotplot also indicated a strong association between the expression of EMT markers and ferroptosis markers in four original sites. Darker dots suggest higher average expressions.

Next, we detected the markers of EMT and ferroptosis in the four original samples (Figure 3B; Supplementary Table S1). The ferroptosis suppressor (CAV1, GPX4, ISCU, ATG5) and ferroptosis driver (HIF1A, KEAP1, PRKAA1, PRKAA2, SCP2, SLC38A1, VEGFA) were obtained through FerrDB database (http://www.zhounan.org/ferrdb/) (Zhou and Bao, 2020). EMT markers (downregulated: CDH1; upregulated: CDH2, ZEB1, SANI1, TWIST1, VIM) were selected according to previous literature (Thiery et al., 2009; Zhang et al., 2021a). We used dotplot to display the above markers expressions in different original samples (Figure 3B). We found a strong correlation between the expression of EMT markers and ferroptosis markers in metastasis and primary tumors. Compared with the primary tumor, in the metastasis group, the mesenchymal markers (CDH2, VIM, TWIST1, SNAI1, ZEB1) and ferroptosis drivers increased, and the epithelial marker CDH1 and ferroptosis suppressor decreased. This result once again confirmed from the level of molecular markers that the relationship between EMT and ferroptosis tendency played a vital role in tumor metastasis process.

Causal Inference Network Suggested Ferroptosis Tendency Promoted the Occurrence of EMT

In order to explore the transcriptional regulatory network inside the double-high group, we then used SCENIC to infer transcription factors (TFs) regulatory information underlying each EC phenotype. SCENIC analysis revealed that some TFs had obvious differential activation and deactivation in the double-high group. We then compared the double-high group with other groups, and detected 5 up-regulated and 18 downregulated TFs (Figure 4A). In the double-high group, we can clearly see that the HOX family genes (HOXA1, HOXB5, HOXC9) is highly expressed. HOX family genes regulate the proliferation, migration and invasion of tumor cells by regulating the transcription of target genes, and participate in the occurrence and development of tumors and drug resistance (Zhang et al., 2003; Bhatlekar et al., 2014; Zhang et al., 2018). In addition, CBX3 is also closely related to metastasis and poor prognosis (Zhong et al., 2019). Further, decreased expression of down-regulated TFs, such as (GATA3, TRPS1, NFKB1, FOXA2, KLF9, etc.), also indicates the influence of the double-high group on poor prognosis and metastasis (Tang et al., 2011; Chou et al., 2013; Gong et al., 2018; Low et al., 2020; Wang et al., 2021).

FIGURE 4
www.frontiersin.org

FIGURE 4. Ferroptosis tendency may promote the occurrence of EMT. (A) SCENIC analysis revealed transcription factors regulatory information in the four groups. 5 up-regulated (red) and 18 down-regulated (blue) transcription factors were detected in the double-high group. (B) The cause-effect relationships among transcription factors, ferroptosis and EMT markers. (C,D) The causal inference network shows regulatory relationships demonstrated around CDH2 and GPX4. GPX4 can negatively regulate the expression of CDH2, while CDH1 can up-regulate the expression of GPX4.

In order to further analyze the regulatory relationship among transcription factors, ferroptosis and EMT markers, we first constructed an undirected network graph through String (version 11.5) and visualized using Cytoscape (version 3.8.2) software (Supplementary Figure S1A). By deleting disconnected nodes in the network, we used connected genes as variables for causal inference (Supplementary Table S2). Figure 4B provided an overview of directed network graph to describe the cause-effect relationships between selected genes. It was apparent from Figures 4C,D that GPX4 as an ferroptosis suppressor could negatively regulate the expression of CDH2, which means that promoting ferroptosis tendency will promote the EMT of malignant epithelial cells. In addition, CDH1 as an epithelial marker can also up-regulate the expression of GPX4 and promote the tendency of ferroptosis. Some other regulatory relationships (SCP2, HOXB5, CAV1) also confirm that ferroptosis can positively regulate the occurrence of EMT (Supplementary Figures S1B–S1D). According to previous literature, during the EMT process, the high transcription level of SNAI1, TWIST1 and ZEB1 will increase the sensitivity of cells to ferroptosis (Chen et al., 2021). Taken together, we believe that EMT of epithelial cells will increase the sensitivity of ferroptosis, and the occurrence of ferroptosis can enhance the transformation of epithelial cells to mesenchymal cells, which form a positive feedback mechanism loop (Xu et al., 2021a; Sun et al., 2021). To validate the clinical relevance of molecules in the gene regulatory network, we performed a validation of the ferroptosis marker in the TCGA-LUAD database. GPX4, as a suppressor of ferroptosis, has an HR of 0.74, p < 0.05. While PRKAA2 as a driver of ferroptosis has an HR of 1.37, p < 0.05 (Supplementary Figures S1E,F). In other words, some markers of ferroptosis can lead to poor prognosis in lung adenocarcinoma, which confirms that the promotion of EMT by ferroptosis leads to poor prognosis.

A Certain Level of the Intracellular Reactive Oxygen Species (ROS) is Generated in Double-High Group

In order to further analyze the differences between the double-high group and other groups, we first screened DEGs through the Findmarker function in Seurat (Supplementary Table S3). Next, we used ClusterProfiler 4.0 for GO enrichment analysis (Figure 5A). The enrichment analysis results showed that some ROS-related biology progresses and molecular functions were activated (“neutrophil degranulation,” “response to reactive oxygen species,” “calium-dependent protein binding,” “antioxidant activity,” etc.) (Dou et al., 2015; Clemens et al., 2017; Sakurada et al., 2019; Yang et al., 2019). We used ssGSEA to score the enrichment of two ROS-related genesets (“HOUSTIS_ROS,” “MODULE_310”), which also shown that the accumulation of ROS in the double-high group was higher (Figure 5B).

FIGURE 5
www.frontiersin.org

FIGURE 5. A higher accumulation of ROS occurred in the double-high group. (A) GO enrichment analysis showed some activated ROS-related biology progresses and molecular functions (“neutrophil degranulation,” “response to reactive oxygen species,” “calium-dependent protein binding,” “antioxidant activity”). (B) ssGSEA indicated higher enrichment scores of two ROS-related genesets in the double-high group.

Several reports have shown that ROS cause changes in downstream molecular such as Src, FAK and 15-PGDH, leading to the activation of EMT in multiple tumors (Wang et al., 2018; Cheung et al., 2020; Peng et al., 2020). Interestingly, the process of ferroptosis was often accompanied by an increase in ROS level and a decrease in the antioxidant system (Zheng and Conrad, 2020; Yao et al., 2021a). As the site for intracellular oxidative respiration, mitochondrion is a major source of ROS. By promoting mitochondrial respiration, the ferroptosis promoter erastin can increase the production of ROS and thus promote ferroptosis by enhancing mitochondrial respiration (Yao et al., 2021b). Combined with the results of the above causal inference, we believe that ferroptosis leads to the low expression of GPX4, which promotes the accumulation of ROS and reduces the inhibition of CDH2 expression, thereby promoting the occurrence of EMT. According to previous literature, we believe that there is a critical value for this promoting. When the strength of promoting ferroptosis is higher than that of EMT, malignant epithelial cells are more likely to undergo programmed death in situ instead of EMT.

Cells Communications Revealed Proliferative Activity and Non-Response to Immunotherapy of Double-High Group

To describe the molecular association between epithelial cells and lymphocytes, we first constructed cell communication networks based on the interaction of potential receptor-ligand pairs using iTALK. Depending on the function of receptor-ligand pairs, we evaluate cellular interactions from three modules: “Growth Factor,” “Cytokine,” and “Immune Checkpoint” (Figures 6A–C).

FIGURE 6
www.frontiersin.org

FIGURE 6. The double-high group showed higher proliferative activity and non-response to immunotherapy. (A) Intracellular and intercellular communications is more active in the double-high group in growth factor module of cell communication networks. The expression of VEGFA, HBEGF, TGFB1, CD9 is higher than other groups. (B) Cytokine module of cell communication networks showed higher expressions of SDC1 and SDC4 in the double-high group. (C) In immune checkpoint module of cell communication networks, the interaction of LGALS9-HARVCR2 between tumor cells and immune cells is more activated in the double-high group. (D,E) The cell cycle and E2F targets scores of the double-high group were also significantly higher than those of the other groups. (F,G) Survival curves showed that the double-high group was linked to a poorer overall survival, and is even worse than the “Ferro_Low_EMT_High” group. (H,I) Both two LUAD immunotherapy cohorts (GSE126044, GSE165252) showed that the non-response group has more double-high group cells than the response group.

In growth factor module, it is clear that there are significant intracellular and intercellular communications in the double-high group (Figure 6A). This result suggests the stimulating effect of autocrine or paracrine growth factors on double-high group cells. Compared to other groups, the double-high group has more VEGFA, HBEGF, TGFB1, CD9 expression. We found that double-high group cells potentially stimulated tumor growth and metastasis through angiogenesis signaling molecules (VEGFA-ITGB1, VEGFA-EGFR, VEGFA-ITGAV) (Maae et al., 2012; Chou et al., 2013; Kim et al., 2017; Sullivan et al., 2019). A specific activation of TGFB1-CXCR4 of the double-high group revealed a high degree of EMT (Xu et al., 2009). The high expression of CD9 also suggests that cancer cells will be more aggressive and rapidly form large tumors (Wang et al., 2019a). Because of the activation of multiple growth factor pathways, we further analyzed the cell cycle (Figures 6D,E). We found that the proliferation of the double-high group also increased significantly. In cytokine module, we found that SDC1 and SDC4 had higher expressions in the double-high group than in other groups (Figure 6B). The high expression of SDC1 and SDC4 indicates proliferation, metastasis and drug resistance of tumor cells, which will result in a poor prognosis (Yao et al., 2019; Kim et al., 2021). Pearson correlation was applied to analyse the correlation of ssGSEA sore of the malignant epithelial cells (Supplementary Figure S2A).

In the immune checkpoint module (Figure 6C), we found that an increase in ferroptosis tendency score would lead to a specific interaction of LGALS9-HARVCR2 between tumor cells and immune cells, with more activation in the double-high group than in the single-high (Sun et al., 2020; Zhang et al., 2021b). A large number of studies have shown that the activation of LGALS9-HARVCR2 in various tumor tissues is related to the development of tumors and negatively regulation of immune signaling. The activation of the LGALS9-HAVCR2 signal can co-express with PD-1 in tumor-infiltrating immune cells, and cooperate to mediate effector T cell exhaustion and dysfunction (Bi et al., 2021; Dixon et al., 2021). This result suggests that double-high group cells may cause non-response to immunotherapy.

Evaluate the Impact of Double-High Group on Survival and Immunotherapy Efficacy

To investigate the clinical role of the cell classification identified in this study, we used CIBERSORTx to evaluate the scores of each cell type in the TCGA LUAD cohort (Supplementary Table S4). “Ferro_High_EMT_Low” group suggested a good prognosis, while the “Ferro_Low_EMT_Low” group (double-low group) had no significant difference in OS (Supplementary Figure S2B). High EMT was related to poor overall survival (OS) (Figures 6F,G). However, it was striking that the cells in the double-high group resulted in worse OS than the “Ferro_Low_EMT_High” group. It was worth noting that the double-high group also had a higher enrichment score in patients with metastases (Supplementary Figures S2C,S2D). The above results once again confirmed that ferroptosis may have a positive regulation on EMT.

Cell communication results indicated that there was a close relationship between immune cell and double-high group cells. Therefore, we chose LUAD immunotherapy cohorts (GSE126044, GSE165252) to verify the therapeutic effect. Consistent with the results of cell communication, more double-high group cells were seen in the immunotherapy non-response group (Figures 6H,I). Our results predicted that the double-high group cells promote tumor metastasis and lead to a poor prognosis, and cause non-response to immunotherapy.

Discussion

Ferroptosis and EMT are closely related as two important biological processes, but the relationship between the two needs to be further confirmed. Therefore, this article used scRNA-seq technology to analyze the correlation between EMT and ferroptosis tendency from the single-cell level. Then, we used causal inference network to objectively interpret the promotion effect of ferroptosis activation on EMT in more details.

Some literature has proved EMT-related signaling contributes to ferroptosis in vitro experiments. ZEB1 is one of the major transcription factors that regulate EMT by binding to the E-box in E-cadherin (Larsen et al., 2016). In the mesenchymal state of cells, ZEB1 correlates with the sensitivity of GPX4 inhibition and thus increases the susceptibility to ferroptosis (Chen et al., 2021). On the contrary, E-cadherin can mediate cell-to-cell contact and inhibit EMT while protecting cells from ferroptosis (Wu et al., 2019). Susceptibility to ferroptosis can be restored when expression of EMT-promoting transcription factors such as SNAI1, TWIST1 and ZEB1 is elevated. On the contrary, ferroptosis tendency can enhance the development of EMT. When PTSG2, as a driver of ferroptosis, was blocked, the expression of N-cadherin was suppressed and the expression of E-cadherin was promoted, which suggested that inhibition of ferroptosis can lead to inhibition of EMT in parallel (Xu et al., 2021a). Not coincidentally, the ferroptosis promoter Erastin reduced the expression of E-cadherin, again suggesting that ferroptosis tendency may be engaged in the development of EMT (Sun et al., 2021). However, the above results were obtained from cell line experiments and were not validated in vivo experiments. Our study used single-cell sequencing technology, with samples from human tissues, which is closer to the clinical patient’s condition and has a greater reference value for clinical treatment. At the same time, we confirmed that ferroptosis tendency has a strong effect on EMT through more objective and true causal inferences. We used difference analysis to identify DEGs between the double-high cells and other cells, and enrichment analysis suggested that ROS played an important role in the mutual promoting of EMT and ferroptosis (Figures 5A,B).

In the study of transcriptional regulation, we used the technique of SCENIC and causal inference. We presented GPX4 plays a negative regulatory role on CDH2. There are also some other meaningful regulatory relationships in the transcriptional regulatory network, which also reveals that ferroptosis tendency can promote EMT (Supplementary Figure S1B). SCP2, as a ferroptosis driver, was found to promote other ferroptosis drivers (PRKAA1, PRKAA2) and EMT marker (CDH2), while inhibiting ferroptosis suppressors such as GPX4 and CAV1. Transcription factors play a key role in the regulation of biological behavior, and HOXB5 is one of them. In lung cancer, it was observed in vitro experiment that proliferation and metastasis of NSCLC cell lines were significantly inhibited when HOXB5 was knocked down (Zhang et al., 2018). Similar results can be found in HNSCC. It has been suggested that HOXB5 regulates EGFR transcription by binding to the EGFR promoter region and therefore it could act as an oncogenic driver in HNSCC (Lee et al., 2020). We found that HOXB5 also promoted the activation of EMT and ferroptosis in the regulatory network (Supplementary Figure S1C). However, there are currently few reports of a link between ferroptosis and HOXB5, and we believe that this may be a direction that can be explored in the future. CAV1 promotes the expression of GPX4 and CDH1, and inhibits CDH2, which implies that it suppresses ferroptosis and EMT (Supplementary Figure S1D). Previous studies have also reported that decrease of CAV1 can promote EMT, more strongly validating the reliability of our results (Zhang et al., 2016). We also verified the correlation between each molecule in bulk sequencing, which is consistent with our results (Supplementary Figure S3).

Immunotherapy has been an important treatment for lung cancer in recent years, but patients’ response to it has been variable and the reasons are not sufficiently clear. By deconvolution of single cell data, we found an immunosuppressive effect of double-high group cells in immunotherapy. In our study, we have seen that LGALS9-HARVCR2 was significantly expressed in both the double-high group and the Ferro_High_EMT_Low group, suggesting that elevation of ferroptosis may play a role in the suppression of immunotherapy. We can also observe that the expression of LGALS9 was higher in the double-high group than in Ferro_High_EMT_Low group, demonstrating a stronger inhibition of immunotherapy in the double-high group. This may provide some clues to the poor efficacy of immunotherapy in some patients and can be a reliable predictor of the risk of immunotherapy. Some of these internal molecular modulations may provide a good direction for future clinical decisions on immunotherapy. Another phenomenon associated with the immune microenvironment is that in tumorigenesis, ferroptosis plays a dual role in promoting and inhibiting tumors, which is attributed to the release of damage-associated molecular patterns (DAMPs). Immunogenic cell death occurs as a result of the release of DAMPs, thereby stimulating anti-tumor immunity (Galluzzi et al., 2017). At the same time, however, the inflammatory response is promoted, thus supporting the growth of tumors (Tang et al., 2012; Chen et al., 2021).

Notably, the current single-cell methods for constructing gene regulatory networks (GRN) using causal inference are time-series-based expression and scRNA-seq expression. MIIC is a scRNA-seq expression profile-based causal inference method, which is based on mutual information for causal inference. This method is essentially based on curve fitting and may have statistical error with the real results. Although the EMT marker SNAI1 received inhibition from SCP2, which is different from the finding that ferroptosis have a promotive effect on EMT, we consider this situation as an uncertain association that does not affect the final results. MIIC, as a scRNA-seq expression profile-based method for constructing gene regulatory networks, can achieve a high accuracy and still shows better results in our GRN.

In summary, we demonstrated the contribution of ferroptosis tendency to EMT through a variety of advanced technologies, and we found that the mutual promotion between them may lead to tumor metastasis, growth and non-response to immunotherapy.

Materials and Methods

Samples and Epithelial Cells Extraction

The data used for this research were downloaded from online public database, which present no ethical issues. We obtained processed single-cell LUAD data from GEO data sets (GSE131907) (Kim et al., 2020). The clinical information of the patients (Supplementary Table S5) and metadata (Supplementary Table S6) come from the original study.

To evaluate EMT enrich score of each epithelial cell, we extracted all epithelial cells (ECs) for subsequent analysis. Because the number of ECs derived from PE is small, we extract normal lung tissues (nLung), early-stage lung cancers (tLung), advanved-stage lung cancers (tL/B), metastatic lymph nodes (mLN) and brain metastases (mBrain) for subsequent analysis (Figure 1A).

Single Cell RNA-Seq Data Quality Control and Data Processing

According to the original articles, malignant tumor cells were separated from non-malignant cells in tumor tissues and three quality measures were applied on raw matrix for each cell: mitochondrial genes (≤20%, unique molecular identifiers (UMIs), and gene count (ranging from 100 to 150,000 and 200 to 10,000). A total of 2000 highly variable genes identified by the FindVariableFeatures function in the Seurat (version 4.0.3) package were used for principal component analysis (PCA)-based dimensionality reduction with RunPCA (Hao et al., 2021). t-distributed stochastic neighbor embedding (t-SNE) was utilized to visualize single-cell clustering using principal components (PCs) 1 to 30. Dotplot function was used to display the expression of EMT and ferroptosis markers in different samples.

Transcription Factor Regulatory Analysis and PPI Network

Single-cell regulatory network inference and clustering (SCENIC) is a new computational method to map genes regulatory networks and identify stable cell states by evaluating the activity of each cell from scRNA-seq data (Aibar et al., 2017). SCENIC was performed on 2000 cells sampled from all ECs, and the regulons were calculated based on transcription factors (TFs) or their target genes. Only regulons significantly upregulated or downregulated were involved in further analysis.

We constructed a protein-protein interaction network for a total of 37 genes (5 upregulated and 14 downregulated TFs, 12 ferroptosis markers, 6 EMT markers) using String version 11.5 (https://string-db.org) (Szklarczyk et al., 2019), and visualized with the force-directed layout algorithm by Cytoscape software (version 3.8.2) (Shannon et al., 2003).

Causal Inference Network

Multivariate Information-based Inductive Causation (MIIC) algorithm relies on a novel information-theoretic method that combines constraint-based learning approach and maximum likelihood framework (Sella et al., 2017; Affeldt et al., 2016; Verny et al., 2017). MIIC online server (https://miic.curie.fr) aims to learn the most appropriate causal, non-causal or mixed graphical model from the available data (Sella et al., 2017). In brief, starting from a fully connected graph, MIIC iteratively removes dispensable edges, by uncovering significant information contributions from indirect paths. This amounts to progressively uncover the best supported conditional independencies by iteratively ‘taking off’ the most significant indirect contributions of positive conditional 3-point information, from every 2-point (mutual) information.

In order to infer the causal relationship between transcription factors from SCENIC and ferroptosis and EMT markers, we extracted TPM data of malignant tumor cells for further analysis. We used huge. npn function of the huge R package to implement the Gaussianization to help relax the assumption of normality (Zhao et al., 2012). We prepared our dataset formatted as a table with variable names (gene names) specified as column names. Category order data was uploaded with a fourth column named “group” added to regroup. When running MIIC on the web interface, we set the following parameters:

Neff: It is Effective number of independent samples, which is set to −1.

Seed: Value used to initialize the pseudorandom number generator in the random sampling phase which is performed only if an effective number of samples is set.

Complexity: Complexity criterion to take into account finite size effects from N or Neff samples for the network reconstruction. It is set to NML.

Orientation: It is set to yes, which means the orientation of v-structures should be enabled.

Propagation: It is set to yes, which means orientations should be propagated downstream of v-structure.

Latent: It is set to no, which means it is not allowed to detect the effects of unobserved (latent) com-mon causes on the relationships between observed variables, represented by bidirected edges.

Gene ontology (GO) Analysis and Single-Sample Gene set Enrichment Analysis (ssGSEA)

Differentially expressed genes (DEGs) of subgroups were recognized by the Findmarker function in Seurat. ClusterProfiler (version 4.0.0) was applied to analyze differences between two groups (Wu et al., 2021b). Based on the GO database, functional enrichment analysis of marker genes was performed on each cell to explore their potential biological functions. ssGSEA was applied to calculate the biological function scores of each cell using the R package GSVA (1.38.2) with the “ssgsea” option for the method argument (Hänzelmann et al., 2013). According to the median scores of EMT and Ferroptosis, we divided all ECs into four groups: Ferro_High_EMT_High, Ferro_High_EMT_Low, Ferro_Low_EMT_High, Ferro_Low_EMT_Low.

Cell–Cell Communication Analysis With iTALK

We mapped the receptor-ligand pairs using iTALK (version 0.1.0) between cell subsets we identified and immune cells to identify cell–cell communications (Wang et al., 2019b). iTALK annotates ligand-receptors into 4 categories: cytokines, growth factors, immune checkpoints and others. For this analysis, we applied two filtering steps of cells inclusion criteria (Chen et al., 2021): Excluded normal lung cells derived from nLung and included all malignant ECs derived from tumors tissues (Puram et al., 2017), Included immune cells identified in the original article: “NK,” “Cytotoxic CD8+ T,” “Exhausted CD8+ T,” “Treg,” “CD4+ Th.” Considering the better graphics of the chord diagram, we choose to display the top 25 receptor-ligand pairs.

Deconvolution of scRNA-Seq Data and Correlation to Public Datasets

We extracted the expression matrices of the four groups of cells from the scRNA-seq data. In order to use the cell types in our scRNA-seq to generate the CIBERSORTx feature matrices, we ran the “Create Signature Matrix” module (https://cibersortx.stanford.edu/runcibersortx.php) (Newman et al., 2019). We used the generated feature matrix to perform CIBERSORTx deconvolution on the bulk RNA-seq cohorts. Among the bulk RNA-seq cohorts, TCGA_LUAD was downloaded from GDC Data Portal (https://portal.gdc.cancer.gov/projects/TCGA-LUAD) to observe the overall survival (OS) and TNM stage, while GSE126044 and GSE165252 were included to observe the therapeutic effect of immunotherapy (Cho et al., 2020; van den Ende et al., 2021).

Data Availability Statement

The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding authors.

Author Contributions

ZL, MZ, and TL conceived and managed the study. JY performed data analysis and visualization. YZ organized literatures and wrote the manuscript. ML performed data analysis. ZS collected materials. All authors reviewed the manuscript and consented for publication.

Funding

This study was supported by The Natural Science Foundation of Liaoning Province of China (2021-MS-179); National Natural Science Foundation of China (No.81302023); Scientific Research Foundation of Education Department in Liaoning Province (Grant LJKZ0739); Science and Technology Plan Project of Shenyang (Grant 21-173-9-30).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fcell.2021.822315/full#supplementary-material

Supplementary Figure S1 | Some meaningful regulatory relationships. (A) undirected network graph visualized using Cytoscape. (B–D) The cause-effect relationships between some genes. (B) SCP2 promoted PRKAA1, PRKAA2 and CDH2, while inhibiting GPX4 and CAV1. (C) HOXB5 promoted CDH2 and PRKAA1, PRKAA2, and is meanwhile promoted by HIF1A, KEAP1 and ATG5. (D) CAV1 promoted the expression of GPX4 and CDH1, and inhibits CDH2. (E,F) Kaplan-Meier survival curves of GPX4 and PRKAA2.

Supplementary Figure S2 | (A) The correlation of ssGSEA sore of the malignant epithelial cells. (B) The scores of each cell type in the TCGA LUAD cohort. “Ferro_Low_EMT_Low” group (double-low group) had no significant difference in OS and “Ferro_High_EMT_Low” group suggested a good prognosis. (C,D) The double-high group had a higher enrichment score in patients with metastases. (C) Patients with N1 & N2 stages have significantly higher enrichment scores than those with N0 stage. (D) Patients in stage M1 received significantly higher scores than those in stage M0.

Supplementary Figure S3 | Analysis of molecular correlations in TCGA dataset.

References

Affeldt, S., Verny, L., and Isambert, H. (2016). 3off2: A Network Reconstruction Algorithm Based on 2-point and 3-point Information Statistics. BMC Bioinformatics 17 (Suppl. 2), 12. doi:10.1186/s12859-015-0856-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Aibar, S., González-Blas, C. B., Moerman, T., Huynh-Thu, V. A., Imrichova, H., Hulselmans, G., et al. (2017). SCENIC: Single-Cell Regulatory Network Inference and Clustering. Nat. Methods 14 (11), 1083–1086. doi:10.1038/nmeth.4463

PubMed Abstract | CrossRef Full Text | Google Scholar

Altman, N., and Krzywinski, M. (2015). Association, Correlation and Causation. Nat. Methods 12 (10), 899–900. doi:10.1038/nmeth.3587

PubMed Abstract | CrossRef Full Text | Google Scholar

Bhatlekar, S., Fields, J. Z., and Boman, B. M. (2014). HOX Genes and Their Role in the Development of Human Cancers. J. Mol. Med. 92 (8), 811–823. doi:10.1007/s00109-014-1181-y

CrossRef Full Text | Google Scholar

Bi, J., Yang, S., Li, L., Dai, Q., Borcherding, N., Wagner, B. A., et al. (2019). Metadherin Enhances Vulnerability of Cancer Cells to Ferroptosis. Cell Death Dis. 10 (10), 682. doi:10.1038/s41419-019-1897-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Bi, K., He, M. X., Bakouny, Z., Kanodia, A., Napolitano, S., Wu, J., et al. (2021). Tumor and Immune Reprogramming during Immunotherapy in Advanced Renal Cell Carcinoma. Cancer Cell 39 (5), 649–661. doi:10.1016/j.ccell.2021.02.015

PubMed Abstract | CrossRef Full Text | Google Scholar

Bulle, A., and Lim, K.-H. (2020). Beyond Just a Tight Fortress: Contribution of Stroma to Epithelial-Mesenchymal Transition in Pancreatic Cancer. Sig Transduct Target. Ther. 5 (1), 249. doi:10.1038/s41392-020-00341-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, X., Kang, R., Kroemer, G., and Tang, D. (2021). Broadening Horizons: the Role of Ferroptosis in Cancer. Nat. Rev. Clin. Oncol. 18 (5), 280–296. doi:10.1038/s41571-020-00462-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Cheung, E. C., DeNicola, G. M., Nixon, C., Blyth, K., Labuschagne, C. F., Tuveson, D. A., et al. (2020). Dynamic ROS Control by TIGAR Regulates the Initiation and Progression of Pancreatic Cancer. Cancer Cell 37 (2), 168–182. doi:10.1016/j.ccell.2019.12.012

PubMed Abstract | CrossRef Full Text | Google Scholar

Cho, J.-W., Hong, M. H., Ha, S.-J., Kim, Y.-J., Cho, B. C., Lee, I., et al. (2020). Genome-wide Identification of Differentially Methylated Promoters and Enhancers Associated with Response to Anti-PD-1 Therapy in Non-small Cell Lung Cancer. Exp. Mol. Med. 52 (9), 1550–1563. doi:10.1038/s12276-020-00493-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Chou, J., Lin, J. H., Brenot, A., Kim, J.-w., Provot, S., and Werb, Z. (2013). GATA3 Suppresses Metastasis and Modulates the Tumour Microenvironment by Regulating microRNA-29b Expression. Nat. Cel Biol. 15 (2), 201–213. doi:10.1038/ncb2672

CrossRef Full Text | Google Scholar

Clemens, R. A., Chong, J., Grimes, D., Hu, Y., and Lowell, C. A. (2017). STIM1 and STIM2 Cooperatively Regulate Mouse Neutrophil Store-Operated Calcium Entry and Cytokine Production. Blood 130 (13), 1565–1577. doi:10.1182/blood-2016-11-751230

PubMed Abstract | CrossRef Full Text | Google Scholar

Demuynck, R., Efimova, I., Naessens, F., and Krysko, D. V. (2021). Immunogenic Ferroptosis and where to Find it?. J. Immunother. Cancer 9 (12), e003430. doi:10.1136/jitc-2021-003430

CrossRef Full Text | Google Scholar

Dixon, K. O., Tabaka, M., Schramm, M. A., Xiao, S., Tang, R., Dionne, D., et al. (2021). TIM-3 Restrains Anti-tumour Immunity by Regulating Inflammasome Activation. Nature 595 (7865), 101–106. doi:10.1038/s41586-021-03626-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Dixon, S. J., Lemberg, K. M., Lamprecht, M. R., Skouta, R., Zaitsev, E. M., Gleason, C. E., et al. (2012). Ferroptosis: an Iron-dependent Form of Nonapoptotic Cell Death. Cell 149 (5), 1060–1072. doi:10.1016/j.cell.2012.03.042

PubMed Abstract | CrossRef Full Text | Google Scholar

Dongre, A., and Weinberg, R. A. (2019). New Insights into the Mechanisms of Epithelial-Mesenchymal Transition and Implications for Cancer. Nat. Rev. Mol. Cel Biol. 20 (2), 69–84. doi:10.1038/s41580-018-0080-4

CrossRef Full Text | Google Scholar

Dou, Q. Q., Rengaramchandran, A., Selvan, S. T., Paulmurugan, R., and Zhang, Y. (2015). Core - Shell Upconversion Nanoparticle - Semiconductor Heterostructures for Photodynamic Therapy. Sci. Rep. 5, 8252. doi:10.1038/srep08252

PubMed Abstract | CrossRef Full Text | Google Scholar

Du, Y., Miao, W., Jiang, X., Cao, J., Wang, B., Wang, Y., et al. (2021). The Epithelial to Mesenchymal Transition Related Gene Calumenin Is an Adverse Prognostic Factor of Bladder Cancer Correlated with Tumor Microenvironment Remodeling, Gene Mutation, and Ferroptosis. Front. Oncol. 11, 683951. doi:10.3389/fonc.2021.683951

PubMed Abstract | CrossRef Full Text | Google Scholar

Galluzzi, L., Buqué, A., Kepp, O., Zitvogel, L., and Kroemer, G. (2017). Immunogenic Cell Death in Cancer and Infectious Disease. Nat. Rev. Immunol. 17 (2), 97–111. doi:10.1038/nri.2016.107

PubMed Abstract | CrossRef Full Text | Google Scholar

Gong, X., Liu, W., Wu, L., Ma, Z., Wang, Y., Yu, S., et al. (2018). Transcriptional Repressor GATA Binding 1-mediated Repression of SRY-Box 2 Expression Suppresses Cancer Stem Cell Functions and Tumor Initiation. J. Biol. Chem. 293 (48), 18646–18654. doi:10.1074/jbc.RA118.003983

CrossRef Full Text | Google Scholar

Hänzelmann, S., Castelo, R., and Guinney, J. (2013). GSVA: Gene Set Variation Analysis for Microarray and RNA-Seq Data. BMC Bioinformatics 14, 7. doi:10.1186/1471-2105-14-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Hao, Y., Hao, S., Andersen-Nissen, E., Mauck, W. M., Zheng, S., Butler, A., et al. (2021). Integrated Analysis of Multimodal Single-Cell Data. Cell 184 (13), 3573–3587. doi:10.1016/j.cell.2021.04.048

PubMed Abstract | CrossRef Full Text | Google Scholar

Kim, M., Jang, K., Miller, P., Picon-Ruiz, M., Yeasky, T. M., El-Ashry, D., et al. (2017). VEGFA Links Self-Renewal and Metastasis by Inducing Sox2 to Repress miR-452, Driving Slug. Oncogene 36 (36), 5199–5211. doi:10.1038/onc.2017.4

PubMed Abstract | CrossRef Full Text | Google Scholar

Kim, N., Kim, H. K., Lee, K., Hong, Y., Cho, J. H., Choi, J. W., et al. (2020). Single-cell RNA Sequencing Demonstrates the Molecular and Cellular Reprogramming of Metastatic Lung Adenocarcinoma. Nat. Commun. 11 (1), 2285. doi:10.1038/s41467-020-16164-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Kim, S., Han, Y., Kim, S. I., Lee, J., Jo, H., Wang, W., et al. (2021). Computational Modeling of Malignant Ascites Reveals CCL5-SDC4 Interaction in the Immune Microenvironment of Ovarian Cancer. Mol. Carcinogenesis 60 (5), 297–312. doi:10.1002/mc.23289

CrossRef Full Text | Google Scholar

Larsen, J. E., Nathan, V., Osborne, J. K., Farrow, R. K., Deb, D., Sullivan, J. P., et al. (2016). ZEB1 Drives Epithelial-To-Mesenchymal Transition in Lung Cancer. J. Clin. Invest. 126 (9), 3219–3235. doi:10.1172/jci76725

CrossRef Full Text | Google Scholar

Lee, K., Chang, J. W., Oh, C., Liu, L., Jung, S.-N., Won, H.-R., et al. (2020). HOXB5 Acts as an Oncogenic Driver in Head and Neck Squamous Cell Carcinoma via EGFR/Akt/Wnt/β-catenin Signaling axis. Eur. J. Surg. Oncol. 46 (6), 1066–1073. doi:10.1016/j.ejso.2019.12.009

CrossRef Full Text | Google Scholar

Liu, S., Thennavan, A., Garay, J. P., Marron, J. S., and Perou, C. M. (2021). MultiK: an Automated Tool to Determine Optimal Cluster Numbers in Single-Cell RNA Sequencing Data. Genome Biol. 22 (1), 232. doi:10.1186/s13059-021-02445-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Low, J. T., Christie, M., Ernst, M., Dumoutier, L., Preaudet, A., Ni, Y., et al. (2020). Loss of NFKB1 Results in Expression of Tumor Necrosis Factor and Activation of Signal Transducer and Activator of Transcription 1 to Promote Gastric Tumorigenesis in Mice. Gastroenterology 159 (4), 1444–1458. doi:10.1053/j.gastro.2020.06.039

PubMed Abstract | CrossRef Full Text | Google Scholar

Maae, E., Andersen, R. F., Steffensen, K. D., Jakobsen, E. H., Brandslund, I., Sørensen, F. B., et al. (2012). Prognostic Impact of VEGFA Germline Polymorphisms in Patients with HER2-Positive Primary Breast Cancer. Anticancer Res. 32 (9), 3619–3627.

PubMed Abstract | Google Scholar

Newman, A. M., Steen, C. B., Liu, C. L., Gentles, A. J., Chaudhuri, A. A., Scherer, F., et al. (2019). Determining Cell Type Abundance and Expression from Bulk Tissues with Digital Cytometry. Nat. Biotechnol. 37 (7), 773–782. doi:10.1038/s41587-019-0114-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Panchy, N., Azeredo-Tseng, C., Luo, M., Randall, N., and Hong, T. (2019). Integrative Transcriptomic Analysis Reveals a Multiphasic Epithelial-Mesenchymal Spectrum in Cancer and Non-tumorigenic Cells. Front. Oncol. 9, 1479. doi:10.3389/fonc.2019.01479

PubMed Abstract | CrossRef Full Text | Google Scholar

Pastushenko, I., and Blanpain, C. (2019). EMT Transition States during Tumor Progression and Metastasis. Trends Cel Biol. 29 (3), 212–226. doi:10.1016/j.tcb.2018.12.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Pearl, J. (2009). Causal Inference in Statistics: An Overview. Statist. Surv. 3 (none), 96–146. doi:10.1214/09-SS057

CrossRef Full Text | Google Scholar

Peng, Y.-S., Syu, J.-P., Wang, S.-D., Pan, P.-C., and Kung, H.-N. (2020). BSA-bounded P-Cresyl Sulfate Potentiates the Malignancy of Bladder Carcinoma by Triggering Cell Migration and EMT through the ROS/Src/FAK Signaling Pathway. Cell Biol Toxicol. 36 (4), 287–300. doi:10.1007/s10565-019-09509-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Puram, S. V., Tirosh, I., Parikh, A. S., Patel, A. P., Yizhak, K., Gillespie, S., et al. (2017). Single-Cell Transcriptomic Analysis of Primary and Metastatic Tumor Ecosystems in Head and Neck Cancer. Cell 171 (7), 1611–1624. doi:10.1016/j.cell.2017.10.044

PubMed Abstract | CrossRef Full Text | Google Scholar

Sakurada, R., Odagiri, K., Hakamata, A., Kamiya, C., Wei, J., and Watanabe, H. (2019). Calcium Release from Endoplasmic Reticulum Involves Calmodulin-Mediated NADPH Oxidase-Derived Reactive Oxygen Species Production in Endothelial Cells. Ijms 20 (7), 1644. doi:10.3390/ijms20071644

PubMed Abstract | CrossRef Full Text | Google Scholar

Sella, N., Verny, L., Uguzzoni, G., Affeldt, S., and Isambert, H. (2017). MIIC Online: a Web Server to Reconstruct Causal or Non-causal Networks from Non-perturbative Data. Bioinformatics 34 (13), 2311–2313. doi:10.1093/bioinformatics/btx844

PubMed Abstract | CrossRef Full Text | Google Scholar

Shannon, P., Markiel, A., Ozier, O., Baliga, N. S., Wang, J. T., Ramage, D., et al. (2003). Cytoscape: a Software Environment for Integrated Models of Biomolecular Interaction Networks. Genome Res. 13 (11), 2498–2504. doi:10.1101/gr.1239303

PubMed Abstract | CrossRef Full Text | Google Scholar

Soundararajan, R., Fradette, J., Konen, J., Moulder, S., Zhang, X., Gibbons, D., et al. (2019). Targeting the Interplay between Epithelial-To-Mesenchymal-Transition and the Immune System for Effective Immunotherapy. Cancers 11 (5), 714. doi:10.3390/cancers11050714

PubMed Abstract | CrossRef Full Text | Google Scholar

Subramanian, A., Tamayo, P., Mootha, V. K., Mukherjee, S., Ebert, B. L., Gillette, M. A., et al. (2005). Gene Set Enrichment Analysis: A Knowledge-Based Approach for Interpreting Genome-wide Expression Profiles. Proc. Natl. Acad. Sci. 102 (43), 15545–15550. doi:10.1073/pnas.0506580102

PubMed Abstract | CrossRef Full Text | Google Scholar

Sullivan, I., Riera, P., Andrés, M., Altés, A., Majem, M., Blanco, R., et al. (2019). Prognostic Effect of VEGF Gene Variants in Metastatic Non-small-cell Lung Cancer Patients. Angiogenesis 22 (3), 433–440. doi:10.1007/s10456-019-09668-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Sun, J., Zhang, Z., Bao, S., Yan, C., Hou, P., Wu, N., et al. (2020). Identification of Tumor Immune Infiltration-Associated lncRNAs for Improving Prognosis and Immunotherapy Response of Patients with Non-small Cell Lung Cancer. J. Immunother. Cancer 8 (1), e000110. doi:10.1136/jitc-2019-000110

CrossRef Full Text | Google Scholar

Sun, L., Dong, H., Zhang, W., Wang, N., Ni, N., Bai, X., et al. (2021). Lipid Peroxidation, GSH Depletion, and SLC7A11 Inhibition Are Common Causes of EMT and Ferroptosis in A549 Cells, but Different in Specific Mechanisms. DNA Cel Biol. 40 (2), 172–183. doi:10.1089/dna.2020.5730

CrossRef Full Text | Google Scholar

Szklarczyk, D., Gable, A. L., Lyon, D., Junge, A., Wyder, S., Huerta-Cepas, J., et al. (2019). STRING V11: Protein-Protein Association Networks with Increased Coverage, Supporting Functional Discovery in Genome-wide Experimental Datasets. Nucleic Acids Res. 47 (D1), D607–D613. doi:10.1093/nar/gky1131

PubMed Abstract | CrossRef Full Text | Google Scholar

Tang, D., Kang, R., Coyne, C. B., Zeh, H. J., and Lotze, M. T. (2012). PAMPs and DAMPs: Signal 0s that spur Autophagy and Immunity. Immunol. Rev. 249 (1), 158–175. doi:10.1111/j.1600-065X.2012.01146.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Tang, Y., Shu, G., Yuan, X., Jing, N., and Song, J. (2011). FOXA2 Functions as a Suppressor of Tumor Metastasis by Inhibition of Epithelial-To-Mesenchymal Transition in Human Lung Cancers. Cell Res. 21 (2), 316–326. doi:10.1038/cr.2010.126

PubMed Abstract | CrossRef Full Text | Google Scholar

Thiery, J. P., Acloque, H., Huang, R. Y. J., and Nieto, M. A. (2009). Epithelial-mesenchymal Transitions in Development and Disease. Cell 139 (5), 871–890. doi:10.1016/j.cell.2009.11.007

PubMed Abstract | CrossRef Full Text | Google Scholar

van den Ende, T., de Clercq, N. C., van Berge Henegouwen, M. I., Gisbertz, S. S., Geijsen, E. D., Verhoeven, R. H. A., et al. (2021). Neoadjuvant Chemoradiotherapy Combined with Atezolizumab for Resectable Esophageal Adenocarcinoma: A Single-Arm Phase II Feasibility Trial (PERFECT). Clin. Cancer Res. 27 (12), 3351–3359. doi:10.1158/1078-0432.Ccr-20-4443

PubMed Abstract | CrossRef Full Text | Google Scholar

Vemuri, V. K. (2015). Counterfactuals and Causal Inference: Methods and Principles for Social Research, by Stephen L. Morgan and Christopher Winship. J. Inf. Tech. Case Appl. Res. 17 (1), 61–64. doi:10.1080/15228053.2015.1014752

CrossRef Full Text | Google Scholar

Verny, L., Sella, N., Affeldt, S., Singh, P. P., and Isambert, H. (2017). Learning Causal Networks with Latent Variables from Multivariate Information in Genomic Data. Plos Comput. Biol. 13 (10), e1005662. doi:10.1371/journal.pcbi.1005662

PubMed Abstract | CrossRef Full Text | Google Scholar

Viswanathan, V. S., Ryan, M. J., Dhruv, H. D., Gill, S., Eichhoff, O. M., Seashore-Ludlow, B., et al. (2017). Dependency of a Therapy-Resistant State of Cancer Cells on a Lipid Peroxidase Pathway. Nature 547 (7664), 453–457. doi:10.1038/nature23007

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, K., Liu, S., Dou, Z., Zhang, S., and Yang, X. (2021). Loss of Krüppel‐like Factor 9 Facilitates Stemness in Ovarian Cancer Ascites‐derived Multicellular Spheroids via Notch1/slug Signaling. Cancer Sci. 112, 4220–4233. doi:10.1111/cas.15100

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, V. M.-Y., Ferreira, R. M. M., Almagro, J., Evan, T., Legrave, N., Zaw Thin, M., et al. (2019). CD9 Identifies Pancreatic Cancer Stem Cells and Modulates Glutamine Metabolism to Fuel Tumour Growth. Nat. Cel Biol. 21 (11), 1425–1435. doi:10.1038/s41556-019-0407-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, W., Hu, Y., Wang, X., Wang, Q., and Deng, H. (2018). ROS-mediated 15-Hydroxyprostaglandin Dehydrogenase Degradation via Cysteine Oxidation Promotes NAD+-Mediated Epithelial-Mesenchymal Transition. Cel Chem. Biol. 25 (3), 255–261. doi:10.1016/j.chembiol.2017.12.008

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, Y., Wang, R., Zhang, S., Song, S., Jiang, C., Han, G., et al. (2019). iTALK: an R Package to Characterize and Illustrate Intercellular Communication:507871. doi:10.1101/507871

CrossRef Full Text | Google Scholar

Wu, H., Wang, Y., Tong, L., Yan, H., and Sun, Z. (2021). Global Research Trends of Ferroptosis: A Rapidly Evolving Field with Enormous Potential. Front. Cel Dev. Biol. 9, 646311. doi:10.3389/fcell.2021.646311

CrossRef Full Text | Google Scholar

Wu, J., Minikes, A. M., Gao, M., Bian, H., Li, Y., Stockwell, B. R., et al. (2019). Intercellular Interaction Dictates Cancer Cell Ferroptosis via NF2-YAP Signalling. Nature 572 (7769), 402–406. doi:10.1038/s41586-019-1426-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Wu, T., Hu, E., Xu, S., Chen, M., Guo, P., Dai, Z., et al. (2021). clusterProfiler 4.0: A Universal Enrichment Tool for Interpreting Omics Data. The Innovation 2 (3), 100141. doi:10.1016/j.xinn.2021.100141

PubMed Abstract | CrossRef Full Text | Google Scholar

Xu, F., Zhang, Z., Zhao, Y., Zhou, Y., Pei, H., and Bai, L. (2021). Bioinformatic Mining and Validation of the Effects of Ferroptosis Regulators on the Prognosis and Progression of Pancreatic Adenocarcinoma. Gene 795, 145804. doi:10.1016/j.gene.2021.145804

PubMed Abstract | CrossRef Full Text | Google Scholar

Xu, J., Lamouille, S., and Derynck, R. (2009). TGF-β-induced Epithelial to Mesenchymal Transition. Cel Res. 19 (2), 156–172. doi:10.1038/cr.2009.5

CrossRef Full Text | Google Scholar

Xu, J., Liao, K., Yang, X., Wu, C., Wu, W., and Han, S. (2021). Using Single-Cell Sequencing Technology to Detect Circulating Tumor Cells in Solid Tumors. Mol. Cancer 20 (1), 104. doi:10.1186/s12943-021-01392-w

PubMed Abstract | CrossRef Full Text | Google Scholar

Yang, M., Akbar, U., and Mohan, C. (2019). Curcumin in Autoimmune and Rheumatic Diseases. Nutrients 11 (5), 1004. doi:10.3390/nu11051004

PubMed Abstract | CrossRef Full Text | Google Scholar

Yao, W., Rose, J. L., Wang, W., Seth, S., Jiang, H., Taguchi, A., et al. (2019). Syndecan 1 Is a Critical Mediator of Macropinocytosis in Pancreatic Cancer. Nature 568 (7752), 410–414. doi:10.1038/s41586-019-1062-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Yao, X., Li, W., Fang, D., Xiao, C., Wu, X., Li, M., et al. (2021). Emerging Roles of Energy Metabolism in Ferroptosis Regulation of Tumor Cells. Adv. Sci. 8 (22), 2100997. doi:10.1002/advs.202100997

CrossRef Full Text | Google Scholar

Yao, Y., Chen, Z., Zhang, H., Chen, C., Zeng, M., Yunis, J., et al. (2021). Selenium-GPX4 axis Protects Follicular Helper T Cells from Ferroptosis. Nat. Immunol. 22 (9), 1127–1139. doi:10.1038/s41590-021-00996-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, B., Li, N., and Zhang, H. (2018). Knockdown of Homeobox B5 (HOXB5) Inhibits Cell Proliferation, Migration, and Invasion in Non-small Cell Lung Cancer Cells through Inactivation of the Wnt/β-Catenin Pathway. Oncol. Res. 26 (1), 37–44. doi:10.3727/096504017X14900530835262

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, K., Yang, G., Wu, W., Zhang, J., Xia, X., Jiang, T., et al. (2016). Decreased Expression of Caveolin-1 and E-Cadherin Correlates with the Clinicopathologic Features of Gastric Cancer and the EMT Process. Pra 11 (2), 236–244. doi:10.2174/1574892811666160128151437

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, N., Ng, A. S., Cai, S., Li, Q., Yang, L., and Kerr, D. (2021). Novel Therapeutic Strategies: Targeting Epithelial-Mesenchymal Transition in Colorectal Cancer. Lancet Oncol. 22 (8), e358–e368. doi:10.1016/s1470-2045(21)00343-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, X., Zhu, T., Chen, Y., Mertani, H. C., Lee, K.-O., and Lobie, P. E. (2003). Human Growth Hormone-Regulated HOXA1 Is a Human Mammary Epithelial Oncogene. J. Biol. Chem. 278 (9), 7580–7590. doi:10.1074/jbc.M212050200

CrossRef Full Text | Google Scholar

Zhang, Z., Bao, S., Yan, C., Hou, P., Zhou, M., and Sun, J. (2021). Computational Principles and Practice for Decoding Immune Contexture in the Tumor Microenvironment. Brief Bioinform 22 (3). doi:10.1093/bib/bbaa075

CrossRef Full Text | Google Scholar

Zhao, T., Liu, H., Roeder, K., Lafferty, J., and Wasserman, L. (2012). The Huge Package for High-Dimensional Undirected Graph Estimation in R. J. Mach Learn. Res. 13, 1059–1062.

PubMed Abstract | Google Scholar

Zheng, J., and Conrad, M. (2020). The Metabolic Underpinnings of Ferroptosis. Cel Metab. 32 (6), 920–937. doi:10.1016/j.cmet.2020.10.011

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhong, X., Kan, A., Zhang, W., Zhou, J., Zhang, H., Chen, J., et al. (2019). CBX3/HP1γ Promotes Tumor Proliferation and Predicts Poor Survival in Hepatocellular Carcinoma. Aging 11 (15), 5483–5497. doi:10.18632/aging.102132

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhou, N., and Bao, J. (2020). FerrDb: a Manually Curated Resource for Regulators and Markers of Ferroptosis and Ferroptosis-Disease Associations, 2020, 2020. doi:10.1093/database/baaa021

PubMed Abstract | CrossRef Full Text | Google Scholar

Keywords: lung adenocarcinoma, single-cell sequencing, ferroptosis, epithelial-mesenchymal transition, causal inference

Citation: Yao J, Zhang Y, Li M, Sun Z, Liu T, Zhao M and Li Z (2022) Single-Cell RNA-Seq Reveals the Promoting Role of Ferroptosis Tendency During Lung Adenocarcinoma EMT Progression. Front. Cell Dev. Biol. 9:822315. doi: 10.3389/fcell.2021.822315

Received: 25 November 2021; Accepted: 30 December 2021;
Published: 20 January 2022.

Edited by:

Shengli Li, Shanghai Jiao Tong University, China

Reviewed by:

Guangjian Liu, Guangdong Provincial People’s Hospital, China
Desi Shang, Harbin Medical University, China

Copyright © 2022 Yao, Zhang, Li, Sun, Liu, Zhao and Li. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

*Correspondence: Tao Liu, bGl1dGFvZGFuaWVsQDE2My5jb20=; Mingfang Zhao, emhhb21mNjE4QDEyNi5jb20=; Zhi Li, emxpQGNtdS5lZHUuY24=

These authors have contributed equally to this work

Disclaimer: All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article or claim that may be made by its manufacturer is not guaranteed or endorsed by the publisher.