Heterogeneity and Cancer-Related Features in Lymphangioleiomyomatosis Cells and Tissue

This article is featured in Highlights of This Issue, p. 1793

Lymphangioleiomyomatosis (LAM), a rare multisystem disease that almost exclusively affects women, is characterized by cystic lung destruction, which can lead to respiratory failure in severe cases (1, 2). Damage to the lung parenchyma is caused by proliferation of spindle-shaped and epithelioid LAM cells, which express alpha-smooth muscle actin (αSMA). Extrapulmonary LAM manifestations frequently include renal angiomyolipomas (AML) and lymphatic abnormalities (1, 2). Empirically, LAM is defined as a low-grade metastasizing neoplasm whose cell (or cells) of origin is uncertain (3–5). Recent analyses of single-cell gene expression profiles indicate that LAM cells may originate in the uterus (6) and/or lung mesenchyme (7). However, the similarity with cancer molecular portraits and the extent of disease heterogeneity are not yet fully understood.

LAM can occur sporadically or in the presence of tuberous sclerosis complex (TSC), an autosomal-dominant multisystem disorder caused by mutations in the tumor suppressor genes TSC1 and TSC2 (8). In sporadic LAM (S-LAM), somatic inactivation of TSC2 or much less commonly TSC1 in an unknown cell type(s) appears to be sufficient for disease development (9). Mutations of TSC1 or TSC2 cause abnormal activation of the mTOR, which is the basis for the current, FDA-approved, standard-of-care use of an mTOR inhibitor for lung and kidney disease (10). Rapamycin has greatly improved outcomes for women with LAM, but some patients continue to lose lung function, albeit at a lower rate, while receiving treatment (11, 12). Disease heterogeneity might contribute to the differences in clinical benefit of rapamycin (13).

mTOR kinase activity is abnormally enhanced in many cancer types, and generally linked to stem cell-like features, which are also present in LAM cells (14–17). Indeed, LAM presents several other fundamental hallmarks of cancer, including continued cell proliferation and resistance to cell death, expression of factors promoting tissue invasion and metastasis, and immune evasion (4, 18). In addition, LAM presents features of steroid-sensitive cancers: LAM lung lesions commonly appear heterogeneous, but diseased cells are frequently positive for the expression of estrogen and progesterone receptors (ER and PR, respectively; refs. 19, 20); LAM disease is generally more progressive in premenopausal women and may be exacerbated by pregnancy or estrogen-based medications (1, 2). Estrogens have been shown to boost the metastatic potential in a LAM cell model (21), and progesterone and estradiol synergistically interact in this process (22). The similarities between LAM and cancer might extend to shared disease susceptibility (23, 24) and, overall, suggest that comparison with cancer molecular portraits can provide further insight into LAM biology. Using expression of gene and protein LAM signatures, and by analyzing LAM single-cell and tissue gene expression profiles, we investigated biological features of LAM and their similarities with cancer. This integrative study proposes causes and markers of LAM cell and tissue heterogeneity, reveals potential transcriptional regulators, and defines molecular and cellular relationships with cancer.

Proteins expressed in LAM lung lesions were compiled from the literature using PubMed searches (up to July 2020) with the MeSH terms “lymphangioleiomyomatosis” and “lung,” “tissue,”, “protein,” “expression,” and/or “IHC,” and through curation of the relevant publications. Table 1 lists the compiled proteins making up the “LAM protein” (LAM^p) signature; the publication sources are detailed in Supplementary Table S1A. The study of LAM^p was accompanied by parallel analyses of a LAM^core signature, which emerged from the analysis of LAM-diseased cells using single-cell RNA sequencing (scRNA-seq) assays in four LAM lung tissue samples (6). The signature expression levels or scores were derived from the combined expression analysis of the corresponding gene constituents using the single-sample gene set expression analysis (ssGSEA) algorithm (25), calculated with the Gene Set Variation Analysis (GSVA) application (https://github.com/rcastelo/GSVA; ref. 26). The same methodology was applied to other gene sets analyzed in this study. A negative signature score for a given cell or tissue means that the gene set has a lower level of expression than the same gene set with a positive value in a different sample, or a different gene set with a positive value (25). The threshold used to define high and low scores of LAM^core and LAM^p in LAM cells was based on the corresponding average signature scores. The GSEA (27) tool was applied using standard parameters and curated C2 gene sets.

Gene expression data from 14 LAM lung bulk tissue samples and control-defined cell lines were obtained from the Gene Expression Omnibus (GEO) reference GSE12027 (28). This report used two cell lines, melanoma (Malme-3M) and pulmonary-artery smooth muscle, as controls or references because LAM cells are phenotypically similar to smooth muscle cells, and contain the premelanosomal structures found in melanocytes and melanoma cells (29). The CIBERSORTx algorithm (30) was applied to bulk gene expression profiles to infer the proportions of immune cell subsets. The Cancer Genome Atlas (TCGA) normal and primary tumor bulk tissue RNA-seq data (RNA-seq v2 with expectation maximization quantification, RSEM) were obtained from the Genomic Data Commons Data Portal (https://portal.gdc.cancer.gov). “Normal” refers to histologically normal tissue and tissue adjacent to surgically removed tumors, as defined by TCGA. Metastases and recurrent tumors were excluded from the study. The cancer types are named in accordance with TCGA study abbreviations. For all TCGA studies, gene expression values were log-transformed, and genes with less than 75% representation in a sample were filtered out.

The scRNA-seq LAM lung data were analyzed by implementing the process described in the original study (6), including preprocessing of LAM1–4 samples, and using the Seurat v3 package (31). Data were downloaded from the GEO reference GSE135851. Briefly, only cells with ≥ 500 expressed [unique molecular identifier (UMI) > 0] genes and < 10% UMIs mapping to mitochondrial genes were considered, except for the analysis of tissue sample LAM4, which included cells with ≥ 300 expressed (UMI > 0) genes and < 10% UMIs mapping to mitochondrial genes. Expression of each gene in a given cell was determined by computing the natural logarithm of (10,000 UMI_gene/UMI_cell +1). Highly variable genes (n = 2,000) were used for dimension reduction, applying runPCA (PCA, n = 200) and the Uniform Manifold Approximation and Projection (UMAP) method (32). The previously defined LAM cell cluster was identified on the basis of the expression of specific genes (6). We identified LAM cells (n = 133) in this dataset by applying the same analytic procedure as originally reported, excluding lung mesenchymal cells, and by confirming the expression of defined LAM genes in a precise cell cluster (Supplementary Fig. S1). This analytic process did not enable us to rule out the possible existence of additional diseased cells, but we limited our analyses to the defined cluster for consistency with the original analysis (6).

Cancer cell line data, from a recent study of cell heterogeneity (33), were downloaded from the Broad Institute's single-cell portal (https://singlecell.broadinstitute.org/single_cell). The cancer cell lines of breast, head and neck, kidney, lung, ovarian, pancreatic, skin, stomach, and uterine tissue were selected (n = 23,546 cells; Supplementary Table S1B), and genes that did not pass the original quality control were excluded from subsequent analyses. For each cell line, a Seurat object was created and only genes detected in ≥ 2 cells were considered, with ≥500 UMIs and <5% UMIs mapping to mitochondrial genes. After merging objects, gene expression was determined by computing the natural logarithm of (10,000 UMI_gene/UMI_cell +1). The LAM and cancer cell profiles were integrated using the FindIntegrationAnchors and IntegrateData functions. Highly variable genes (n = 2,000) were also used for dimensional reduction by applying runPCA and UMAP (PCA, n = 200).

Enrichment at transcription factor–binding sites was determined with oPOSSUM v3.0 (http://opossum.cisreg.ca), using default parameters and JASPAR core profiles, and promoter genomic sequences set to −2 kb from the transcription start site. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses were performed in g:Profiler in R (https://biit.cs.ut.ee/gprofiler/page/r) using the default parameters for an unranked gene list.

Formalin-fixed paraffin-embedded LAM lung tissue from seven patients was analyzed in this study by IHC assays. These samples corresponded to lung biopsies or organ transplants of individuals with a clinical LAM diagnosis based on outlined guidelines (1), and collected by Spanish LAM reference hospitals (University Hospital Vall d'Hebron; University Hospital La Princesa; University Hospital Clínica Puerta del Hierro; and University Hospital of Bellvitge), coordinated by the Spanish LAM Association (AELAM). LAM diagnosis included in all cases transbronchial lung tissue biopsy and subsequent pathologic evaluation, high-resolution CT scanning, and pulmonary function tests. All patients provided written, informed consent and the biomarker study was approved by the ethics committee of IDIBELL. The studies conformed to the principles set out in the World Medical Association Declaration of Helsinki and the Department of Health and Human Services Belmont Report.

Assays were performed on serial paraffin sections using an EnVision kit (Dako). Antigens were retrieved using citrate pH 6.0 (anti-LUM, anti-LCN2, and anti-RUNX1) or EDTA pH 9.0 (anti-IRF1 and anti-YB1) buffer. Endogenous peroxidase was blocked by pre-incubation in a solution of 3% H₂O₂, performed in 1× phosphate-buffered saline with 10% goat serum. Slides were incubated overnight at 4°C with primary antibody diluted in blocking solution. The antibodies were anti-IRF1 (C-20, Santa Cruz Biotechnology; dilution 1:100), anti-LCN2 (MAB1757, R&D Systems; dilution 1:50), anti-LUM (sc-166871, Santa Cruz Biotechnology; dilution 1:50), anti-RUNX1 (ab23980, Abcam; dilution 1:200), and anti-YB1 (ab76149, Abcam; dilution 1:50). Secondary peroxidase-conjugated antibodies (Envision+ system-HRP, Dako) were used. Positive controls are shown in Supplementary Fig. S2. Sections were counterstained with hematoxylin and examined under a Nikon Eclipse 80i microscope.

LAM blood samples were collected during the 2017 and 2018 annual AELAM patient conferences, so the period between undertaking pulmonary function tests and sample acquisition varied, making it unfeasible to assess LCN2 levels by pulmonary function. The data collected consisted of age at diagnosis, age at sample extraction, diagnosis of AML, pneumothorax, TSC, and therapy used at the time of sample extraction. All patients provided written informed consent, and the study was approved by the ethics committees of IDIBELL and the Instituto de Investigación Sanitaria La Princesa, Hospital de Henares, Spain. Control samples were obtained from healthy women. VEGF-D levels were measured using a commercially available ELISA kit according to the manufacturer's protocol (R&D Systems). Emphysema plasma samples were collected from patients attending the Interstitial Lung Disease Unit of the University Hospital of Bellvitge. These were adult individuals (31–59 years old) with a smoking history and dyspnea on exertion, but whose FEV₁/FVC ratio was less than 0.70, and who had not been diagnosed with chronic obstructive pulmonary disease. The plasma samples of patients with Langerhans cell histiocytosis and Sjögren syndrome were collected by the ILD Center of Excellence, St. Antonius Hospital Biobank, Nieuwegein, the Netherlands. The study was approved by the St. Antonius Hospital ethics committee (reference R05–08A), and all participants provided written informed consent. Detailed information about the individuals included in the plasma study is provided in Supplementary Table S1C. The studies conformed to the principles set out in the WMA Declaration of Helsinki and the Department of Health and Human Services Belmont Report.

LCN2 was quantified using the human LCN2/NGAL Quantikine ELISA Kit (R&D Systems) in accordance with the manufacturer's instructions. The samples were diluted 1:100. Plasma samples of patients with cerebrospinal fluid leak were used as positive controls (34). The limits of quantification and detection were 52 and 23 pg/mL, respectively.

The correlation between signature scores was computed using Spearman rank (r_s) correlation coefficient. The two-tailed Mann–Whitney U test was used to assess signature scores and gene expression differences between pairs of groups. The log-fold change of expression differences between two defined cell groups was used to rank genes for GSEA (27). Genes differentially expressed between cell subsets were identified using the FindMarkers function of the Seurat v3 package with default parameters. The bimodality of the LAM^core score distribution in LAM cells was determined using the bimodality coefficient and Hartigan dip statistic (35). All data and reagents in this study are available from the authors upon request.

LAM lung lesions often display heterogeneity at the molecular and cellular levels (36). Even diagnostic markers—typically αSMA and the premelanosomal protein gp100, detected using the HMB-45 mAb—may be heterogeneous among diseased cells (1, 2, 37, 38). Because the activities of genes and proteins are precisely coordinated to execute cellular functions (39), analysis of a set of biomarkers corresponding to proteins expressed in LAM lung lesions could provide further insight into disease biology (40). To define this LAM-associated protein set, we compiled evidence from the literature preceding single-cell LAM studies (6, 7); this curation provided a list of 67 proteins shown to be expressed in lung LAM lesions (Table 1), but whose levels relative to other conditions, tissue and/or cells are generally unclear. This 67-protein set defined the LAM^p signature used in subsequent analyses. Using gene expression data of bulk LAM lung tissue nodules (28), 24 (35%) of the LAM^p genes showed significant overexpression relative to two control cell lines (Supplementary Fig. S3). This percentage of overexpressed LAM^p genes in bulk LAM tissue was found to be higher than would be expected at random considering all the genes examined (one-tailed proportion test, P = 0.029).

To complement the study of LAM^p, we used a transcriptional signature of LAM single cells, named LAM^core, which was obtained by the landmark analysis of scRNA-seq profiles in four LAM lung tissues (6). Of the 777 genes in LAM^core, 14 were identified in LAM^p (hypergeometric test for overlap P < 10^–6; Table 1). Then, using the scRNA-seq data from the LAM^core study (6), we identified 21% of LAM^p genes (14 of 66 available for analysis) as being significantly (FDR < 5%) overexpressed in LAM cells relative to non-LAM cells (Supplementary Table S2A). This proportion was again found to be higher than would be expected at random (hypergeometric test for overlap P = 5 × 10^–15). Therefore, as expected, the literature-curated set of proteins expressed in LAM lung lesions (i.e., LAM^p) and the single-cell LAM cell transcriptome signature (i.e., LAM^core) significantly overlap with respect to their gene identities and differences in expression relative to other lung cell types.

To further assess the LAM^p and LAM^core signature scores relative to the disease, we compared their expression values (excluding shared genes; Table 1) between LAM and non-LAM scRNA-seq profiles (6). We carried out non-parametric, ssGSEA (25) to compute signature values in individual cell transcriptomes. This yielded an enrichment score denoting the degree to which the genes in a given set were upregulated or downregulated in a coordinated manner within a cell. The LAM^p score was computed using 56 (84%) of the literature-curated genes: 11 were not measured or were entirely unexpressed in the dataset (AGTR2, CCR2, CXCR1, DCT, FLT4, ITGB3, MMP9, REN, TNFRSF11B, TNFSF11, and VEGFD). The LAM^core was scored on the basis of all the originally defined genes in this signature (n = 777; ref. 6). As expected, both signatures showed a significantly higher level of expression in LAM than in non-LAM cells (Fig. 1A). Unexpectedly, however, the signatures were found to be negatively correlated in LAM cells (Spearman correlation coefficient (r_s) = −0.35, P = 5 × 10^–5; Fig. 1B, left), but strongly positively correlated in non-LAM cells (r_s = 0.42, P < 2 × 10^–16; Fig. 1B, right). It is of particular note that LAM^core was positively correlated (r_s = 0.81, P < 2 × 10^–16) with a gene set characteristic of LAM cells identified in an independent scRNA-seq study (7).

The differential correlation of LAM^p and LAM^core signatures might indicate that although the two sets partially overlap with respect to their constituents, and are coexpressed in non-LAM cells, they recognize partially different diseased-cell phenotypes. Following up on this idea, the LAM^core signature score was found to be bimodally distributed in LAM cells (Fig. 1C). The bimodality allowed us to define two cell subsets (Fig. 1C), which might be indicative of disease heterogeneity. Evaluation of the LAM^p-LAM^core relationship in these subsets revealed a negative correlation only in LAM cells with high LAM^core values (Fig. 1D). The LAM^core was also bimodally distributed in non-LAM cells (Fig. 1C). All the lung mesenchymal cells, which we expected to be similar to LAM cells (6), were located within the highest 10% of LAM^core values in this analysis. In turn, the LAM^p score distribution did not show a significant difference from unimodality in LAM cells and non-LAM cells (Fig. 1E).

Following on from the observation of two potential LAM^core-based subsets of diseased cells, for each measured gene expression, we computed the log-fold change between the LAM cells with high and low LAM^core values (Supplementary Table S2B), and used these differences to evaluate associations with a curated collection of gene sets representing defined biological states in health and disease (27). This analysis demonstrated the phenotypic variability among LAM cells: the LAM^core score was found to be positively correlated (FDR < 5%) with sets linked to the extracellular matrix organization (“NABA matrisome”), hypoxia under certain conditions (“Elvidge MCF7 hypoxia upregulated”), cell stemness (“Boquest CD31-negative stromal stem cell upregulated,” “Wong adult tissue cell stem module,” and “Lee neural stem cell upregulated”), and cancer cell invasion (“Anastassiou multicancer invasiveness signature”), among other similar associations (Fig. 1F; Supplementary Table S2C). In turn, the LAM^core score in LAM cells was found to be negatively correlated with sets linked to amino acid metabolism (Reactome metabolism of amino acids and derivatives pathway), translation (Reactome eukaryotic translation elongation pathway), and sensitivity to rapamycin (“Bilanges rapamycin sensitive via TSC1/2;” ref. 41), among other similar associations (Fig. 1F; Supplementary Table S2D). An analogous analysis of LAM^p scores in LAM cells did not reveal any correlations that were the opposite of those detected with LAM^core with an FDR < 5%, but several sets were uniquely and significantly correlated with LAM^core or LAM^p (Supplementary Table S2C–S2F): positive correlations with LAM^core that were not identified (nominal P > 0.10) with LAM^p included sets linked to the onco-miR-21 (“Gabriely miR-21 targets in glioma”) and OCT4 (“Benporath OCT4 targets in human embryonic stem cells”) targets, among others (Fig. 1G). On the other hand, uniquely positive correlations with LAM^p included sets linked to the cell cycle (“Whitfield genes expressed in G₂–M phase”), hypoxia (“Harris hypoxia induced”), and oxidative phosphorylation (WikiPathway oxidative phosphorylation), among others (Fig. 1H). Therefore, LAM cells with high or low values of LAM^core and LAM^p signatures are differentially associated with biological processes and functions.

The predicted subsets of LAM diseased cells were further evaluated by IHC studies of the gene products with the largest difference in expression between LAM^core-low and LAM^core-high cells: the keratan sulfate proteoglycan lumican (LUM), predicted to be overexpressed in LAM^core-high cells; and the transcription factor Y-Box binding protein 1 (YB1), predicted to be overexpressed in LAM^core-low cells (Supplementary Table S2B). Both proteins were identified as being expressed in LAM lung lesions, but LUM expression was mostly limited to diseased cells with a spindle cell-like phenotype, while the level of YB1 expression was lower in this type than in epithelioid cells, when detected within the lesions and in the epithelial layer (Fig. 1I).

Given the key role of estrogens and progesterone in LAM biology and disease progression (42), we further evaluated differences between LAM cells with high and low LAM^core relative to the expression of signatures of ER target genes (43) and progesterone-induced genes in ER-positive breast cancer (44). The LAM^core-high, but not the LAM^core-low LAM cells tended to be positively associated with these two signatures (Supplementary Fig. S4). Indeed, among the curated gene sets, a positive association with LAM^core corresponded to genes that were downregulated in primary ovarian cells exposed to progesterone (45) (“Wilcox response to progesterone downregulated”; Fig. 1D). Violin plots further illustrate the overexpression of ER/PR-related signatures in LAM cells with high LAM^core scores (Fig. 1J). Collectively, these data depict two subsets of LAM cells that differ in gene expression profiles of key biological states, which may in turn be connected to cancer features.

To evaluate associations between the defined LAM^p and LAM^core signatures, and cancer molecular portraits, we examined their values (excluding shared genes) in gene expression data of normal-adjacent and primary tumor tissue of 15 solid cancer types studied in TCGA (46). These cancers comprised tissue-of-origin categories proposed for LAM, in particular those of breast [breast invasive carcinoma (BRCA)], kidney [kidney renal clear cell carcinoma (KIRC)], lung [lung adenocarcinoma and lung squamous cell carcinoma; (LUAD) and (LUSC), respectively], and uterus [uterine corpus endometrial carcinoma (UCEC)]. In addition to the positive LAM^p-LAM^core correlation in non-LAM cells (Fig. 1B, right), the scores of these signatures also proved to be positively correlated in most normal tissue settings, and in all tumor types (Supplementary Fig. S5A). LAM^p and LAM^core were mostly underexpressed in tumors relative to normal tissue, with the exception of glioblastoma (GBM), head and neck squamous cell carcinoma (HNSC), pancreatic adenocarcinoma (PAAD), and stomach adenocarcinoma (STAD; Fig. 2). The LAM^core signature did not show any significant differences between normal and tumor tissues in GBM, HNSC, and PAAD, whereas LAM^p showed the opposite trend in these types and in STAD, that is, significantly higher levels of expression in tumors than in normal tissue counterparts (Fig. 2). Intriguingly, these cancer types (i.e., GBM, HNSC, PAAD, and STAD) are characterized by high stromal content and/or mesenchymal cell features, which may broaden the concept that the signatures capture some largely overlapping phenotypes. Nonetheless, both signatures generally showed higher levels of expression in normal tissue, which is consistent with LAM being a low-grade neoplasm (5).

Unsupervised hierarchical clustering of TCGA gene expression data using the LAM^core gene list commonly identified relatively high scores in normal-adjacent tissue of urothelial bladder carcinoma (BLCA), HNSC, KIRC, and STAD studies; in tumors, high signature scores were also widely identified in BLCA, KIRC, and STAD (Supplementary Fig. S5B). An analogous analysis of LAM^p mainly identified high scores in normal-adjacent tissue of BRCA, but tumor analyses yielded similar results to those found in LAM^core, largely identifying high scores in BLCA, HNSC, KIRC, and STAD, as well as in PAAD (Supplementary Fig. S5B). These cancer connections broaden the range of observations described above regarding tumor stromal content and cancer cell phenotype.

The availability of single-cell transcriptome profiles of LAM cells from different patients (6), and of cancer cell lines representing a variety of neoplasms (33), offers the possibility of determining disease similarities more precisely. For this purpose, the data from both studies were processed and integrated for cell clustering analysis (Materials and Methods). The cancer cell lines that were clustered closer to LAM cells originated from breast, gastric, lung, and urothelial cancer (Fig. 3A). These similarities involved six cancer cell lines with at least two cells identified: BT549, corresponding to triple-negative breast cancer with a mesenchymal phenotype (47); LMSU, corresponding to a gastric lymph node metastasis; BFTC909, corresponding to a transitional cell or urothelial carcinoma; and three models of pleural mesothelioma, represented by the ACCMESO1, NCIH226, and NCIH2452 cell lines (Fig. 3B). The BT549 cell line mapped most closely to LAM cells and is a model that is sensitive to mTOR inhibition; targeting mTOR with everolimus in this setting reduces VEGF expression (48). The expression of ER in this cell line is not completely null and can be upregulated using chromatin remodelers (49). In addition, we identified three pleural mesothelioma cell models, which also typically show PI3K/mTOR overactivation and sensitivity to rapamycin (50).

To further assess the similarity between the identified cancer cell lines and LAM cells, we determined which genes were differentially expressed between these two groups and all other cancer cells included in the initial clustering analysis. Fifty genes were identified as commonly overexpressed (FDR < 5%; Supplementary Table S2G) in the six cancer cell lines and LAM cells. This number represented a higher proportion than would be expected by chance: 819 and 225 genes were overexpressed in the LAM cells and the six cancer cell lines, respectively; hypergeometric test P = 1 × 10^–24. The 50-gene set included two known markers of LAM, ACTA2, and MMP2 (both included in LAM^p and LAM^core; Table 1), and the autophagy receptor SQSTM1/p62; of the remaining genes, 18 were also included in LAM^core (Supplementary Table S2G). In addition, analysis of GO terms and pathway annotations in the 50-gene set exposed the strongest associations (FDR < 5%) with extracellular matrix biology (Supplementary Table S2H), which was consistent with the identification of “NABA matrisome” as the highest positive correlation with LAM^core in LAM cells (Fig. 1F). Therefore, single-cell transcriptome analyses connect LAM to four cancer settings, including those of breast cancer, gastric lymph node metastasis, mesothelioma, and urothelial carcinoma.

Having determined the similarities between the LAM^core and LAM^p signatures in single-cell and bulk tissue data, the representation of transcriptional factor-binding motifs in the corresponding gene promoters—covering a length 2 kb upstream from each transcription start site—was evaluated using JASPAR core binding profiles (51). Excluding shared genes, 11 transcription factors were identified as being overrepresented (Z-scores > 5, P = 5 × 10^–7) in both signatures: CEBPA, ELF5, FOXA1, IRF1, NFATC2, RUNX1, SOX5, SOX9, SPIB, STAT1, and TBP (Supplementary Table S2I). STAT1 is a mediator of LAM cell survival and proliferation (52, 53), and ELF5 and FOXA1 are coregulators of ER function (54). SOX9 supports mTOR signaling and lung metastasis in breast cancer (55), and is also expressed in LAM lung lesions (16), which means that it is a component of LAM^p (Table 1).

Of the remaining transcription factor predictions, RUNX1 regulates the fate of ER-positive mammary epithelial luminal cells (56) and its expression characterizes circulating tumor cells of endometrial cancers, in which it promotes invasion and metastasis (57). Intriguingly, the previous results identified LAM^core and LAM^p positive correlations (FDR < 5%) with a set corresponding to genes regulated by RUNX1-RUNX1T1 (ref. 58; Supplementary Table S2C). The transcription factor predictions also included IRF1, which could provide a link to chemokine expression in LAM lung lesions and fluids (28). Therefore, we evaluated the expression of these two factors in the LAM lung tissue of 7 patients using IHC assays. RUNX1 was strongly positive in LAM epithelioid cells in 4 of the patients; diseased cells with a spindle cell-like phenotype showed weak positivity for this factor in all cases (Fig. 4, top). In turn, LAM cells with a spindle cell-like phenotype showed clear IRF1 expression in 5 of the patients and, as expected, infiltrating immune cells were also found to be strongly positive (Fig. 4, bottom).

The above results broaden the concept of heterogeneity among LAM cells. Next, we aimed to assess heterogeneity among LAM lung tissues. For this study, the variance of gene expression measured by microarray assays in 14 LAM lung tissue samples (28) was computed, and the gene decile with the highest variance was defined. This gene subset showed significant overrepresentation of several GO terms corresponding to immune system function, including “Chemokine activity” and “Cytokine activity” (Fig. 5A; Supplementary Table S2J). Several other terms overrepresented in the LAM variable gene set were linked to cell differentiation, development, and DNA transcription (Supplementary Table S2J). In parallel, pathway-based analyses identified “Cytokine-cytokine receptor interaction” as the most overrepresented (Fig. 5A; Supplementary Table S2J), reinforcing the concept that genes with large expression variation among LAM lung tissue are mainly linked to the immune system.

Hierarchical unsupervised clustering with bootstrap resampling using the aforementioned gene decile identified a LAM subgroup comprising four of the 14 lung samples (Fig. 5B). The four-sample subgroup showed a lower level of expression of the signature corresponding to “Bilanges rapamycin sensitive via TSC1/2” (Fig. 5B; ref. 41). Differential gene expression analysis identified 37 overexpressed and six underexpressed genes (FDR < 5%) in this subgroup, and several of the highlighted genes were linked to immune system functions (LCN2, MAP2K3, OSMR, PTPN2, PTPRJ, and TRIM48; Fig. 5C). The total set of 43 genes presented significant overrepresentation of IRF1-binding sites (Z-score = 4.63), but an underrepresentation of RUNX1-binding sites (Z-score = −7.23). We then inferred the immune cell contents (30) in the 14 samples and found that the immune-associated subgroup had significantly higher levels of eosinophil numbers (Mann–Whitney test, P = 0.008; Supplementary Table S2K). Collectively, these data depict two subsets of LAM lung lesions that differ in immune system factors. Therefore, variability of immune cell content and/or of expression of defined factors by LAM and/or tissue microenvironment cells may contribute to disease heterogeneity. It is of particular note that LAM^core in LAM cells showed positive correlations with TNF-induced gene sets (Fig. 1G).

A review of the cancer-related literature featuring the overexpressed genes identified in the immune-associated LAM tissue subgroup highlighted lipocalin 2 (LCN2; Fig. 5C). LCN2 is overexpressed and associated with poor prognosis in several solid malignancies (59), and mediates breast cancer metastasis to lung (60). LCN2 is a multifunctional, innate-immunity protein that protects from inflammation in epithelial tissue (61). IHC analysis of LAM lung tissue (n = 7) revealed strong LCN2 expression in all cases, in epithelioid and spindle-like diseased cells (Fig. 6A). The degree of LCN2 positivity varied between lesions (Fig. 6A). Unexpectedly, however, the level of LCN2 expression in scRNA-seq LAM cells was found to be significantly lower than in non-LAM cells (log-fold change = −0.32, P = 0.001). This apparent inconsistency with the IHC results is reminiscent of the negative correlation between LAM^p and LAM^core, and could indicate that the identities of the diseased cells analyzed by IHC and scRNA-seq assays are not the same.

LCN2 is a plasma biomarker of several diseases linked to tissue inflammation, including acute and chronic kidney disease, chronic obstructive pulmonary disease, and cancers (62). We analyzed LCN2 levels in plasma from 78 patients with LAM, 21 healthy age-matched female controls, and 32 patients with LAM-related pulmonary diseases (12 patients with emphysema, 10 with pulmonary Langerhans cell histiocytosis, and 10 with Sjögren syndrome) using an ELISA. The level of plasma LCN2 was significantly higher in patients with LAM than in healthy women (two-tailed Mann–Whitney test P = 0.047), but there were no significant differences between patients with LAM and those with the other pulmonary diseases (Fig. 6B). Patients with LAM who were not taking rapamycin showed higher LCN2 plasma levels than did healthy women (two-tailed Mann–Whitney test P = 0.003; Fig. 6C). Four LAM cases appeared to be outliers, with relatively high LCN2 levels (Fig. 6B and C); if these cases were ignored, there was no significant difference between patients with LAM and healthy women, although the overabundance of LCN2 in untreated LAM relative to healthy women remained significant (two-tailed Mann–Whitney test P = 0.009). Measuring VEGF-D in LAM plasma has been clinically implemented for disease diagnosis and monitoring (11). Patients with a defined high plasma VEGF-D level (> 800 pg/mL) also showed overabundance of LCN2 relative to healthy women (two-tailed Mann–Whitney test P = 0.041; Fig. 6D), which further connects LCN2 to the pathology.

Using complementary gene/protein signatures, this study indicates sources of LAM cell and tissue heterogeneity, in parallel with the similarities and differences with respect to cancer molecular portraits. A literature-based signature of proteins expressed in LAM lung lesions (LAM^p) captures a diseased-cell phenotype that partially differs from that identified by RNA-seq profiling of LAM single diseased cells [LAM^core (6)]. These signatures are found to be specifically anti-correlated in LAM cells, in particular in diseased cells with high values of LAM^core. The signature distinctions might be influenced by their different origin (i.e., proteins expressed in LAM lung lesions, LAM^p; or genes expressed in LAM lung single cells, LAM^core). However, both signatures significantly overlap in terms of gene identities, are positively coexpressed in most normal and cancer settings, and are overexpressed in LAM single cells relative to non-LAM single cells. As inferred from the curated gene set associations with LAM^p and LAM^core scores, LAM cell subsets differ in terms of their hormone-mediated signaling, proliferation, and stemness features, which in turn may be linked to distinct metabolic activities (e.g., metabolism of amino acids and rapamycin-sensitivity genes). These observations are further evidence of the existence of a phenotypic gradient of diseased cells that may be associated with the established spindle-like and epithelioid forms (13). Paracrine interactions in the tissue microenvironment may also contribute to the gradient. However, TSC1 or TSC2 mutations in profiled single cells need to be identified to assign diseased-cell features definitively.

When compared with gene expression profiles in cancer studies, the LAM^p and LAM^core signatures show frequent associations with normal-adjacent tissue and/or with primary tumors predicted to have high stromal content and/or mesenchymal cell characteristics. These may be sources of bias when trying to assess LAM cell origin using bulk tissue data. Inference of the LAM origin may also be hindered because normal cells are not expected to show abnormal activation of mTOR, even though normal stem cells are frequently characterized by mTOR activation (63). On the other hand, cancer cells commonly harbor many other genetic alterations that cause complex biological changes, despite frequently showing mTOR overactivation (64). In spite of these limitations, we have compared LAM and cancer single-cell transcriptome profiles, the results demonstrating connections with the diseased cell phenotype (e.g., BT549, mesenchymal; LMSU, metastatic) and/or proposed disease origins. The identification of three pleural mesothelioma cell lines with their characteristic overactivation of PI3K/mTOR and sensitivity to mTOR inhibition (50) is particularly intriguing, especially considering the recent demonstration that Tsc2 loss in the lung mesenchymal lineage causes LAM-like disease in mice (7). In parallel, our study exposes RUNX1 as being a regulator of LAM cell signatures. This factor is highly expressed in the mesenchymal and epithelial compartments of the developing and postnatal lung (65), and regulates lung inflammation (66). Moreover, genetic variation in the RUNX1 locus, and in IRF1, has been associated with pulmonary function measures in the general population (refs. 67, 68; Supplementary Fig. S6). Therefore, LAM cells might originate in a lung cell population that is regulated by these factors during tissue development and/or inflammation-related insults.

By analyzing gene expression profiles in bulk LAM lung tissue, our study also indicates sources of disease heterogeneity. This heterogeneity appears to result mainly from the interplay with the immune system, as is supported by the expression of IRF1 and RUNX1 in LAM lung lesions. The activity of the immune system may influence LAM progression and application of immunotherapies (18, 69, 70). A 43-gene set is proposed to differentiate two types of lesions, and the expression of the multifunctional innate immunity protein LCN2 may contribute to this distinction, further linking the disease to inflammation (60, 61). LCN2 is a lipocalin that stabilizes matrix metalloproteinases, which are known to mediate lung tissue destruction in LAM (71–73). Moreover, RUNX1 negatively regulates LCN2 expression and prevents autophagy in skeletal muscle (74), whereas LCN2 expression is induced by IFNγ and TNFα (75); furthermore, LAM^core positively correlates with TNF-induced gene sets. The depicted diseased tissue, as well as single-cell differences might influence the response to rapamycin, but analyses of relatively large cohorts are required to test this hypothesis.

B. Saez reports personal fees from MSD, Janssen, and Chiesi outside the submitted work. C. Valenzuela reports personal fees from Boehringer Ingelheim, F. Hoffmann-La Roche, Ltd., and BMS outside the submitted work. M. Molina-Molina reports grants from Roche and Boehringer Ingelheim, grants and personal fees from Esteve-Teijin, and personal fees from Origo Pharma outside the submitted work. No disclosures were reported by the other authors.

R. Espín: Data curation, software, formal analysis, validation, investigation, visualization, methodology, writing–review and editing. A. Baiges: Data curation, formal analysis, validation, investigation, visualization, methodology, writing–review and editing. E. Blommaert: Formal analysis, investigation, methodology, writing–review and editing. C. Herranz: Data curation, formal analysis, writing–review and editing. A. Roman: Resources, data curation, writing–review and editing. B. Saez: Resources, data curation, writing–review and editing. J. Ancochea: Resources, data curation, writing–review and editing. C. Valenzuela: Resources, data curation, writing–review and editing. P. Ussetti: Resources, data curation, writing–review and editing. R. Laporta: Resources, data curation, writing–review and editing. J.A. Rodriguez-Portal: Resources, data curation, writing–review and editing. C.H.M. van Moorsel: Resources, data curation, writing–review and editing. J.J. van der Vis: Resources, data curation, writing–review and editing. M.J.R. Quanjel: Resources, data curation, writing–review and editing. A. Villar-Piqué: Resources, methodology, writing–review and editing. D. Diaz-Lucena: Resources, methodology, writing–review and editing. F. Llorens: Resources, methodology, writing–review and editing. A. Casanova: Resources, data curation, writing–review and editing. M. Molina-Molina: Resources, data curation, writing–review and editing. M. Plass: Software, investigation, methodology, writing–review and editing. F. Mateo: Resources, investigation, methodology, writing–review and editing. J. Moss: Data curation, investigation, methodology, writing–review and editing. M.A. Pujana: Conceptualization, resources, data curation, formal analysis, supervision, funding acquisition, validation, investigation, visualization, methodology, writing–original draft, project administration, writing–review and editing.

We thank the individuals who provided samples and the AELAM foundation for its continued support for LAM research. We also wish to thank Oscar M. Tirado for support with LUM immunohistochemistry. Our results are partly based upon data generated by TCGA Research Network (https://www.cancer.gov/tcga), and we express our gratitude to TCGA consortia and coordinators for producing the data and the clinical information used in our study. This research was partially supported by AELAM (ICO-IDIBELL agreement, to M.A. Pujana), The LAM Foundation Seed Grant 2019, to M.A. Pujana, Carlos III Institute of Health grant PI18/01029, to M.A. Pujana and ICI19/00047 to M. Molina-Molina [co-funded by European Regional Development Fund (ERDF), a way to build Europe], Generalitat de Catalunya SGR grant 2017-449, to M.A. Pujana, the CERCA Program for IDIBELL institutional support, and ZonMW-TopZorg grant 842002003, to C.H.M. van Moorsel. M. Plass was supported by a “Ramón y Cajal” contract of the Spanish Ministry of Science and Innovation (RYC2018-024564-I) and J. Moss was supported by the Intramural Research Program of NIH/NHLBI.

The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.