Integrated analysis of long non-coding RNAs and mRNA profiles reveals potential sex-dependent biomarkers of bevacizumab/erlotinib response in advanced lung cancer.

BACKGROUND: While lung cancer patient outcomes are well-recognized to vary as a function of patient sex, there has been insufficient research regarding the relationship between patient sex and EGFR(Epidermal growth factor receptor) response efficacy. The present study therefore sought to identify novel sex-related biomarkers of bevacizumab/erlotinib (BE) responses in non-small cell lung cancer (NSCLC) patients.
METHODS: The exon array data in the Gene Expression Omnibus (GEO) dataset were analyzed in order to identify patterns of mRNA and lncRNA expression associated with BE resistance in NSCLC. These differentially expressed (DE) lncRNAs and mRNAs were identified via DE Analysis Filtering. These DE mRNAs were then assessed for their potential functional roles via pathway enrichment analyses, with overlapping functions possibly associated with the BE resistance. The mRNAs in these overlapping groups were then assessed for their correlations with patient survival, and lncRNA-mRNA co-expression networks were generated for each patient subset. A protein-protein interaction (PPI) network was also generated based upon these DE mRNAs.
RESULTS: In females we identified 172 DE lncRNAs and 1766 DE mRNAs associated with BE responses, while in males we identified 78 DE lncRNAs and 485 DE mRNAs associated with such responses. Based on the overlap between these two datasets, we identified a total of 37 GO functions and 18 pathways associated with BE responses. Co-expression and PPI networks suggested that the key lncRNAs and mRNAs associated with these BE response mechanisms weredifferent in the male and female patients.
CONCLUSIONS: This work is the first to conduct a global profiling of the relationship between lncRNA and mRNA expression patterns, patient sex, and BE responses in individuals suffering from NSCLC. Together these results suggest that the integrative lncRNA-mRNA expression analyses may offer invaluable new therapeutic insights that can guide the tailored treatment of lung cancer in order to ensure optimal BE responses.

Introduction

Lung cancer is a leading driver of cancer-associated death globally [1], with most patients only being diagnosed when the disease is in its advanced stages and is no longer eligible to undergo curative surgical resection [2]. In non-small cell lung cancer (NSCLC), epidermal growth factor receptor (EGFR) mutations are frequently detected [3]. Erlotinib was developed as an inhibitor of EGFR to treat patients with NSCLC that have activating mutations in EGFR [4-6], but how effective this compound is in patients with wild type EGFR remains less clear. Some initial studies have suggested that firstline treatment with erlotinib or a similar kinase inhibitor in such patients may even be more harmful that conventional chemotherapy [7], although more recent biomarker analysis of large-scale trial found that a subset of EGFR wild-type patients do benefit from such treatment [8].

A leading therapeutic strategy for treating cancer revolves around the simultaneous targeting of multiple targets via a combination treatment approach in an effort to simultaneously suppress multiple signaling pathways within tumor cells [9, 10]. The combination of erlotinib with the bevacizumab, which is an antibody specific for vascular endothelial growth factor (VEGF) has been found to be more effective than treatment with erlotinib alone [11]. While EGFR activating mutations are the biomarkers most closely associated with patient responses to this combined bevacizumab + erlotinib (BE) treatment, certain patients without such mutations have also been found to benefit from combination therapy. Previous studies have investigated the efficacy and safety of erlotinib plus bevacizumab treatment in patients with non-small-cell lung cancer [12, 13]. However, the exact molecular mechanisms underlying BE responsiveness in these patients remain uncertain. It is therefore vital that additional biomarkers of BE responses be identified in order to more effectively identify those EGFR wild-type patients who are most likely to benefit from BE treatment in order to guide appropriate patient treatment strategies.

Solid tumor biomarkers and therapeutic targets have increasingly been identified using omics technologies [14], with whole-genome exon arrays being ideal for discovering such biomarkers with exon-level resolution [15]. These arrays permit for gene- and exon-level analysis of the expression of both mRNAs and long non-coding RNAs (lncRNAs). These lncRNAs are known to be capable of regulating a variety of activities including RNA splicing, stability, and shuttling between the nucleus and the cytoplasm [16]. The integrated analysis of both mRNA and lncRNA expression patterns in tumor cells has the potential to offer invaluable insight into the molecular basis for tumor pathogenesis, thereby potentially identifying prognostic or therapeutic targets to guide patient treatment efforts. To date, however, no studies have constructed a lncRNA-mRNA interaction network specifically associated with BE therapeutic responses in lung cancer.

There is clear evidence that men typically have poorer lung cancer outcomes than do women [17]. There is additionally evidence suggesting that there are sex-specific differences in immunological functionality and in patient responses to targeted molecular therapies. We therefore hypothesized that such sexual dimorphism may have an impact on the relative efficacy of BE responses in male and female lung cancer patients.

In the present study, we analyzed patient subgroups from among individuals participating in the EGFR wild-type group of the SAKK 19/05 trial in an effort to better understand what factors influence patient responses to BE therapy. In particular, we focused on the relationship between patient sex and genetic factors as they related to BE responses in this patient cohort. For this analysis, we analyzed published whole-exome array data available in the Gene Expression Omnibus (GEO) in order to identify those mRNAs and lncRNAs that were differentially expressed (DE) between patients that did and did not respond to BE therapy in male and female patient subgroups. When then constructed separate lncRNA-mRNA coexpression networks pertaining to BE responses in males and females separately, and analyzed the similarities and differences between the networks in these two groups. Together, the results of our analysis highlight that there are sex-dependent differences in patient responses to targeted molecular therapies [18-20], suggesting that such integrative lncRNA/mRNA-based analyses may help to better predict patient therapeutic responses, thereby allowing clinicians to select appropriate treatments for lung cancer patients in an individually targeted manner.

Materials and methods

Data sets

We obtained publicly available gene expression profile data (GSE37138) from GEO (http://www.ncbi.nlm.nih.gov/geo) for the present analysis. This dataset included baseline samples obtained from the whole bronchoscopic biopsy of 42 suffering from advanced non-squamous NSCLC. These patients had then undergone treatment with BE either until disease progression was observed or until toxicity necessitated the termination of this treatment. Based on whether or not they exhibited stable disease after 12 weeks of treatment, patients were categorized as either responders (n = 22; 12 male, 10 female) or non-responders (n = 20; 12 male, 8 female). Additional detail regarding the SAKK 19/05 trial is discussed in published articles [21]. Data and probe sequences were analyzed in accordance with the GPL5188 [HuEx-1_0-st] Affymetrix Human Exon 1.0 ST Array, with background signals being adjusted according to the Robust Multichip Average.

mRNA and lncRNA profiling

Using the provided exon array data, we sought to assess the mRNA and lncRNA expression profiles in these patient samples. To that end, we began by constructing an mRNA and lncRNA alignment database. The exon array enrichment uniquely mapped to genes and the lncRNAs were indeed inferred from genes. In order to identify mRNAs, those transcripts that were preceded by either NM or XM in the NCBI refseq database and those transcripts annotated to code for protein in the Ensembl database (http://www.genome.ucsc.edu/) were pooled together, with any repeats being removed. In order to identify lncRNAs, those RNAs preceded by either “NR_” or “XR_” in the Refseq database, those that were annotated with “lincRNA”, “processed transcripts”, “non-coding”, or “misc_RNA” in the Ensembl database, and those that were found within the NONCODEv5 database were pooled together. After removal of repeats, those non-coding RNA sequences that were > 199 bp long were then selected to yield a lncRNA alignment database. Next, BLAST (https://blast.ncbi.nlm.nih.gov/Blast.cgi) was carried out with the bidirectional best hit (BBH) criteria to map the 141,1080 probe set IDs for this dataset in order to yield an mRNA alignment database with an e-value threshold of 10e-5. Those probe sets aligned with these mRNAs were annotated as such, with the remaining probes being mapped to the lncRNA alignment database. Only those matched lncRNAs with an e-value threshold of 10e-5 were included in this analysis.

Identification of mRNAs and lncRNAs associated with BE responses

Those mRNAs and lncRNAs that were differentially expressed between responders and non-responders in both the male and female patient subgroups were next identified. These DE lncRNAs and DE mRNAs were identified via DE Analysis filtering, using filtering criteria of fold change >2 and p < 0.05 and fold change >1.2 and p < 0.05, respectively.

Functional enrichment analysis

DE mRNAs were separated into four groups according to the sex and direction in which they were differentially regulated (male-upregulated; male-downregulated; female-upregulated; female-downregulated). The DE mRNAs in these groups then underwent both Gene Ontology (GO: http://www.geneontology.org) and KEGG (http://www.genome.ad.jp/kegg/) pathway analyses using DAVID (http://david.abcc.ncifcrf.gov/). For the GO analyses, the biological processes associated with the input mRNAs in an annotated database were used to identify regulatory gene networks based upon the hierarchical relationships between these genes, while KEGG analyses offered insight into the pathways in which these DE mRNAs participate. Only those GO and KEGG terms meeting the following criteria were included in these analyses: p<0.001, FDR<0.05 [22].

We reasoned that those functional pathways that were selectively enriched in both male and female NSCLC patients may represent the pathways most closely associated with BE responses. We therefore compared such shared pathways between these two sex groups in order to assess whether the same genetic profiles were responsible or these BE response-associated phenotypes in both males and females. The same analysis was used to compare both up- and down-regulated KEGG and GO terms between male and female patients.

lncRNA-mRNA correlation network

Cytoscape (http://www.cytoscape.org) was used to generate two lncRNA-mRNA correlation networks (one each for males and females) based upon the predicted targets of DE lncRNAs and the available DE mRNA data. Normalized signal intensity values for DE lncRNAs and DE mRNAs, as generated above, was used for network construction. In addition, as most lncRNAs included in this analysis have not been functionally annotated, we conducted a co-expression analysis in order to identify those mRNAs that vary their expression in a pattern closely related to the expression pattern for a given lncRNA. Such analyses offer potentially insight into the poorly understood roles for these lncRNAs as regulators of protein-coding gene expression. The resultant co-expression network figure makes note of both the directionality of lncRNA/mRNA differential expression, as well as associated GO and KEGG terms and the correlative relationships between different transcripts. In addition, degree centrality was determined based upon the number of links between nodes, serving as an indicator of the relative biological importance of a given transcript.

PPI network analysis

The STRING database was used to screen for known or predicted interactions between proteins of interest in the present study, with a specific focus on interacting protein pairs associated with DE mRNAs in the lncRNA-mRNA correlation network. Only those interactions with a confidence cutoff > 0.4 were included in the resultant network, which was visualized using Cytoscape.

Survival analysis

We additionally examined the relationship between specific DE mRNAs identified in the present study and patient survival outcomes. To that end, the Cancer Genome Atlas (TCGA) was used to identify samples from LUSC (lung squamous cell carcinoma) and LUAC (lung adenocarcinoma) patients for which clinical follow-up data was available. Patients in these cohorts were then divided into two groups based on whether they expressed high or low levels of a given gene of interest, and the overall survival (OS) between these two groups was then compared via the Kaplan-Meier method, with log-rank tests used to gauge significance.

Statistical analyses

The Bioconductor limma package (v3.26.1) and R studio (v3.2.2) were used to identify DE lncRNAs and DE mRNAs between BE-responsive and -unresponsive patient samples in the present study. FDR values were used to control for multiple testing. P<0.05 was the significance threshold.

Results

Patient characteristics and clinical outcomes

We began by accessing the GSE37138 dataset containing exon array data from 42 NSCLC patients that had undergone targeted BE treatment. These patients were from 35–78 years old (median: 61), with all being of good performance status (WHO 0–1). Of these patients, 24 were male (57.1%) and 18 were female (42.8%), with almost all patients (92.9%) having stage IV disease. EGFR and KRAS mutations were detected in 4 (9.5%) and 6 (14.3%) of patients, respectively.

Identification of differentially expressed lncRNAs and mRNAs

We next compared the expression of lncRNAs and mRNAs in both male and female patients in this dataset in order to identify those transcripts that were differentially expressed as a function of patient BE response. In females, we identified 172 DE lncRNAs (113 up- and 59 down-regulated) and 1766 DE mRNAs (866 up- and 900 down-regulated), while in males we identified 78 DE lncRNAs (44 up- and 34 down-regulated) and 485 DE mRNAs (179 up- and 306 down-regulated). These data were then organized into both volcano plots (Fig 1) and were visualized following hierarchical clustering into a heatmap (Fig 2).

Fig 1

Differential expression of mRNAs and lncRNAs in NSCLC patients as a function of BE responses.

(A) Differences in lncRNA expression in female BE-responsive patients; (B) Differences in lncRNA expression in male BE-responsive patients; (C) Differences in mRNA expression in female BE-responsive patients; (D) Differences in mRNA expression in male BE-responsive patients. The color scale corresponds to relative transcript expression, with green and red corresponding to low and high expression, respectively. BE: bevacizumab/erlotinib.

Fig 2

Hierarchical clustering analysis of differentially expressed lncRNAs and mRNAs.

(A) Differences in lncRNA expression in female BE-responsive patients; (B) Differences in lncRNA expression in male BE-responsive patients; (C) Differences in mRNA expression in female BE-responsive patients; (D) Differences in mRNA expression in male BE-responsive patients. The color scale corresponds to relative transcript expression, with green and red corresponding to low and high expression, respectively. BE: bevacizumab/erlotinib.

Functional enrichment analyses of DE mRNAs

We next conducted GO- and KEGG-based functional enrichment analysis of the DE mRNAs identified above in order to understand the relationship between these genes and BE therapeutic responses. These DE mRNAs were analyzed separately for both males and females.

In the female dataset up-regulated DE mRNAs were enriched for 747 GO terms, including signal transduction (Fig 3A), while the most enriched term corresponding to down-regulated mRNAs was ‘apoptotic process’ (Fig 3B). In males, up-regulated mRNAs were most closely associated with nucleosome assembly, positive regulation of transcription from RNA, transcription from RNA polymerase II promoter, chromatin organization, and positive regulation of transcription, and DNA-templates (Fig 3C), while the most enriched term corresponding to down-regulated mRNAs was neutrophil degranulation (Fig 3D).

Fig 3

The 25 most-enriched GO terms associated with BE response-related differentially expressed mRNAs in lung cancer patients.

(A) Enriched GO terms associated with upregulated mRNAs in females; (B) Enriched GO terms associated with downregulated mRNAs in females; (C) Enriched GO terms associated with upregulated mRNAs in males; (D) Enriched GO terms associated with downregulated mRNAs in males; BE: bevacizumab/erlotinib.

A pathway enrichment analysis revealed that 150 annotated pathways were associated with upregulated DE mRNAs in the female dataset, with the ‘Ras signaling pathway’ being the most enriched (Fig 4A). In addition, 102 annotated pathways were associated with downregulated DE mRNAs in the female dataset, with ‘metabolic pathways’ being the most enriched (Fig 4B). In contrast, upregulated DE mRNAs from the male dataset were enriched in 16 pathways, of which ‘Systemic lupus erythematosus’ was the most strongly enriched (Fig 4C), while downregulated DE mRNAs from the male dataset were enriched in 67 pathways, of which ‘metabolic pathways’ were also the most enriched (Fig 4D).

Fig 4

The 25 most-enriched pathways associated with BE response-related differentially expressed mRNAs in lung cancer patients.

(A) Enriched pathways associated with upregulatedmRNAs in females; (B) Enriched pathways associated with downregulated mRNAs in females; (C) Enriched pathways associated with upregulated mRNAs in males; (D) Enriched pathways associated with downregulated mRNAs in males; BE: bevacizumab/erlotinib.

In order to identify those enriched GO terms and pathways most closely associated with BE responses, we identified overlapping terms between the male and female datasets, revealing 11 and 26 shared GO functions that were up- and down-regulated, respectively, and 4 and 14 pathways that were up- and down-regulated, respectively in both males and females. These overlapping pathways, including ‘apoptotic process’, ‘metabolic pathways’, ‘PI3K-Akt signaling’, and ‘pathways in cancer’ are thus likely associated with BE responsiveness (S1 and S2 Tables, and Fig 5). Among them, the positive/negative regulation of transcription from RNA polymerase II promoterGO functions and metabolic pathways enriched in both up-regulated and down-regulated groups. Despite these similarities with respect to GO and KEGG terms between male and female datasets, however, the genes which gave rise to these overlapping similarities differed in males and females. In the female group, 433 DE mRNAs corresponded to the up-regulated GO function type, 542 DE mRNAs were involved in the down-regulated GO terms and 191 DE mRNAs in both up- and down-regulated GO functions, whereas in males the numbers were 124, 262 and 69, respectively. Similarly, in females 83 and 227 DE mRNAs were associated with the up- and down-regulated KEGG pathways, respectively, whereas in males these numbers were 47 and 109, respectively. In addition, 132 mRNAs in females and 53 in males were included in both up- and down-regulated pathways. These mRNAs were considered to be likely to play roles in regulating patient responses to BE therapy, and so were pooled together for subsequent analyses, yielding 256 DE mRNAs in males and 690 in females (929 total).

Fig 5

Venn diagram analysis comparison of functional enrichment analysis results.

Significantly enriched (A) GO terms and (B) pathways overlapping between the female-upregulated, male-upregulated, female-downregulated, and male-downregulated DE mRNA datasets.

Co-expression of lncRNAs-mRNAs in BE response lung cancer patients

We next constructed a lncRNA-mRNA correlation network in order to better understand how lncRNAs influence BE responses in NSCLC patients, incorporating both the DE lncRNAs and BE response-associated mRNAS identified above into the resultant network (Fig 6A and 6B). In total, the female network was made up of 282 nodes (232 mRNAs and 50 lncRNAs) and 458 connections, while the male network was made up of 202 nodes (177 mRNAs and 25 lncRNAs) and 375 connections. The majority of the co-expressed lncRNAs and mRNAs were positively correlated with one another. In many cases, individual lncRNAs were co-expressed with multiple mRNAs, suggesting the potential for a complex regulatory relationship between these molecules. In order to gauge the relative importance of the lncRNAs and mRNAs in these networks, the degree of centrality was calculated (S3 and S4 Tables). The nodes that possess the highest degree of centrality in females were: TEK, CDH5, ERG, ASPA, LDB2, TIE1, ACVRL1 and EDNRB, while in males the nodes with the highest degree of centrality were: SNX27, BLZF1, ATP6V1A, CHML, STAM, ACOX1, UBQLN1 and CRK. These results thus confirmed that there were sex-dependent differences in lncRNA-mRNA regulatory relationships in NSCLC patient tumor samples.

Fig 6

Co-expression analysis of core mRNAs and correlated lncRNAs.

(A) A lncRNA-mRNA expression correlation network in females. (B) A lncRNA-mRNA expression correlation network in males. (C) A lncRNA-mRNA-GO-pathway correlation network in females. (D) A lncRNA-mRNA-GO-pathway correlation network in males. lncRNAs are represented by triangular nodes, while mRNAs are represented by circular nodes. GO terms are indicated by V-shaped arrow, while pathways are represented with rectangles. Red, blue, and gray nodes indicate upregulation, downregulation, or no significant regulation, respectively. Solid lines indicate positive direct correlations, while dotted lines correspond to negative correlations.

We further conducted the functional annotation of these co-expressed mRNAs in an effort to better understand lncRNA functionality (Fig 6C and 6D). The majority of lncRNA-associated mRNAs in this network were associated with GO functions relating to metabolic pathways and the positive regulation of transcription from RNA polymerase II. The key hub genes identified through these analyses differed between groups, with female hub genes including TEK, ERG, NONHSAT095695.2, LDB2, ENST00000552253.1, ENST00000588495.5, CDH5, and MEF2C, and male hub genes including ENST00000460756.5, G6PC3, ATF6, BMP2, SIX1, RAB7A, ENST00000575695.5, and SOX2 (S5 and S6 Tables). NONHSAT095695.2, ENST00000552253.1, and ENST00000588495.5 primarily regulated LDB2, ENST00000460756.5 primarily regulated ATF6, and ENST00000575695.5 primarily regulated UBQLN1. Together this network analysis thus highlights the potential for sex-specific differences in the drug response mechanisms in NSCLC patients.

PPI network construction

In an additional effort to understand these sex-dependent differences in BE responses, we generated a PPI network for the male and female patient cohorts, incorporating DE mRNAs that had been incorporated into the above lncRNA-mRNA co-expression network. As shown in Fig 7, TEK, LDB2, UBQLN1, and ATF6 all appeared to play key hub roles in their corresponding sex groups.

Fig 7

Protein-protein interaction network containing for differentially expressed mRNAs in the lncRNA-mRNA correlation network.

Proteins are represented by nodes, with edges corresponding to protein-protein interactions. Upregulated and downregulated DE mRNAs are represented in red and blue, respectively. DE mRNAs: differentially expressed mRNAs.

Lastly, we assessed whether the 929 DE mRNAs identified in our functional enrichment analyses were associated with different survival outcomes in LUAD and LUSC patients. To that end, a Kaplan-Meier-based approach was used to assess the prognostic value of these mRNAs in TCGA datasets. In total, we found that 266 of these mRNAs were significantly associated with LUAD patient survival, while 73 were associated with LUSC patient survival. In the LUSC patient cohort, 22 and 91 mRNAs were significantly associated with survival outcomes in females and males, respectively, while in LUAD 129 and 158 mRNAs were significantly associated with OS in females and males, respectively. Of the hub genes identified in females, we found that elevated expression of TEK, which was expressed at higher levels in biopsies from female BE-responsive patients, was correlated with better LUAD patient survival outcomes in females but not in males (Fig 8A and 8B). In contrast, decreased expression of UBQLN1, which was downregulated in biopsies from male BE-responsive patients, was associated with increased survival in male or total LUAD patients, but was unrelated to survival in female LUAD and LUSC patients (Fig 8C and 8D). Together these results highlight the potential relevance of sex-dependent mRNA expression profiles and corresponding differences in drug responses.

Fig 8

The relationship between TEK and UBQLN1 expression and LUAD patient overall survival.

Kaplan-Meier survival curves were used to assess (A) survival as a function of TEK expression in all LUAD patients; (B) survival as a function of TEK expression in female LUAD patients; (C) survival as a function of UBQLN1 expression in all LUAD patients; (D) survival as a function of UBQLN1 expression in male LUAD patients. LUAD: Lung Adenocarcinoma.

Discussion

In the present study, we established profiles of differential mRNA and lncRNA expression that were associated with EGFR wild-type NSCLC patient responses to BE combination therapy. We found that profiles of DE mRNAs and lncRNA differed between male and female patients, although there was clear overlap with respect to the molecular pathways associated with BE responses in both patient groups. We further conducted lncRNA-mRNA co-expression analyses that highlighted sex-dependent heterogeneity with respect to the differential patterns of gene expression underlying patient responses. Survival analyses further confirmed that many of the hub genes identified in these analyses were associated with sex-dependent differences in lung cancer patient survival outcomes.

About 10% of Caucasian NSCLC patients and 40–50% of East Asian NSCLC patients exhibit EGFR mutations [23]. EGFR tyrosine kinase inhibitors (TKIs) have been shown to offer therapeutic benefit even in patients that do no exhibit such activating EGFR mutations [24-26], with a subset of EGFR wild-type patients exhibiting clinical responses to these treatments. The relevance of many biomarkers has been investigated in these EGFR wild-type patients, such as the expression of EGFR, EGFR ligand, an amplification of EGFR, all of which are associated with TKI responses [27]. However, there are currently no clear biomarkers that are suitable for the identification of NSCLC patients that are likely to respond well to EGFR TKI treatment, particularly in EGFR wild-type patients.

In this study assess the relationships between lncRNA and mRNA expression in exon arrays from NSCLC patients in order to explore genetic signatures associated with responses to a combination of an EGFR-TKI and an anti-angiogenic agent. However, more potential differentially expressed long non-coding RNAs and coding genes can be discovered if RNA-seq data is used instead of exon array, we only used exon array data in this work due to the data availability. In total we identified 37 GO functions and 18 pathways that were putatively linked to BE responses, with prominent enriched functions including apoptotic processes, metabolic pathways, cancer-related pathways, and transcription-associated processes. These results suggested that these pathways are closely associated with lung cancer patient BE responses. Previous work has clearly demonstrated roles for these pathways in modulating cancer cell signal transduction, thereby altering the proliferation, growth, survival, and metabolism of tumor cells [28-30]. In a lncRNA-mRNA-GO-Pathway-network analysis, we further demonstrated the interrelatedness of these signaling pathways, highlighting a molecular network governing BE responses that were consistent with past experimental research [31].

Interestingly, we found that distinct lncRNAs and mRNAs were associated with BE responses in our male and female patient cohorts, with TEK, ERG, NONHSAT095695.2, LDB2, ENST00000552253.1, ENST00000588495.5, CDH5 and MEF2C being tightly linked to BE responses in females, and ENST00000460756.5, G6PC3, ATF6, BMP2, SIX1, RAB7A, ENST00000575695.5, and SOX2 being associated with these responses in males. In females, these core genes were upregulated in BE-responsive patients, whereas in males all of these core genes other than SIX1 and SOX2 were downregulated in BE-responsive patients. These results are of interest as they highlight sex-dependent differences in patient tumor responses, consistent with previous work demonstrating sexual dimorphism with respect to cancer patient physiological, pathophysiological, and pharmacological outcomes [32]. Work from the Moffitt Cancer Center indicated that men were associated with worse lung cancer outcomes than women (HR, 1.24 [95% CI, 1.09–1.41]) [33], and sex-related differences in patient immunotherapy have also been observed [34]. Indeed, immune checkpoint inhibitors have been shown to bolster the OS of advanced NSCLC patients in a sex-dependent manner [35]. Women have also been found to exhibit higher rates of EGFR mutations, thus making them more likely to respond to TKI treatment [7]. Work by Seo Yun Kim study 66 LUAC patients suggests that female EGFR wild-type patients are more likely to respond to EGFR TKI (odds ratio [OR], 3.10; 95% CI, 1.05–9.19; p = 0.041) than are males [36]. The specific mechanisms underlying these sex-dependent differences in pharmacological and physiological outcomes in men and women are, however, not well understood.

In this study we generated a lncRNA-mRNA co-expression network which allowed us to highlight male- and female-specific transcriptional signatures associated with BE responses. All of these hub genes have been linked with tumor growth and development, although they have not all been associated with drug responses. The TEK hub gene identified in females has been linked with pathologic complete responses to Bevacizumab in individuals with breast cancer [37]. Work by Song et al. further suggests that ERG is a biomarker predictive of docetaxel responsiveness in those with metastatic castration-resistant prostate cancer [38]. NONHSAT095695.2, ENST00000552253.1 and ENST00000588495.5 regulated LDB2, which inhibits liver cancer cell proliferation and migration via suppressing the expression of HEY1 [39].

Work by Xiaoman indicates that CDH5 can suppress the metastasis of breast cancer owing to the ability of its 3’-UTR serve as a ceRNA for STARD13 [40]. In addition, MEF2C has previously been identified as a possible leukemia-related oncogene [41].

In males, we found that ENST00000460756.5 primarily regulated ATF6, which is linked with signal transduction. In addition, ENST00000575695.5 primarily regulated UBQLN1, which is linked with hypoxia responses in cells. Previous studies have also explored the roles of G6PC3, ATF6, BMP2, SIX1, RAB7A, and SOX2 in tumors [42-47]

Together our results most clearly highlighted the importance of TEK, LDB2, ATF6 and UBQLN1 in NSCLC patient BE responses, underscoring the value of examining system-level gene interactions in order to better understand the determinants of lung cancer patient therapeutic responsiveness.

Responses to BE treatment may be associated with more favorable patient outcomes. We therefore assessed the relationship between BE response-related hub gene expression and lung cancer patient survival. In public microarray datasets, we found that the elevated expression of TEK, which was upregulated in female BE-responsive patients, was associated with increased survival in female lung cancer patients. Similarly, we found that decreased expression of UBQLN1, which was reduced in male BE-responsive patients, was associated with increased survival in male lung cancer patients. LDB2 and ATF6 were not included in TCGA lung cancer datasets. However, LDB2 has been shown to inhibit liver cancer cell proliferation and migration, functioning as a tumor suppressor [39], while ATF6 has been linked to the induction of intestinal dysbiosis and immune activity associated with the development of colorectal cancer [48]. Reduced ATF6 expression in the BE responders was thus likely associated with better survival.

Conclusion

In summary, the present study is the first to conduct global mRNA and lncRNA profiling in order to explore sex-related differences in the molecular mechanisms governing NSCLC patient responses to BE therapy. We identified a large number of differentially expressed lncRNAs and mRNAs when comparing datasets from patients that were and were not responsive to BE treatment. Through predictive and functional enrichment analyses, we were able to further uncover potential links between these lncRNAs and mRNAs, offering potential insight into the functional roles of these lncRNAs in the regulation of BE responses in NSCLC. While the present results highlight the unique cancer-associated pathways associated with BE responses in male and female NSCLC patients, more multi-omic studies will be integrated to our present work to discover more possible lncRNA-mRNA correlation in future, and further studies will be needed to validate these findings and to understand the underlying molecular basis for these findings.

(TIF)

Click here for additional data file.

GOs intersections.

(XLSX)

Click here for additional data file.

Pathways intersections.

(XLSX)

Click here for additional data file.

Nodes in the female lncRNA-mRNA expression correlation network.

(XLSX)

Click here for additional data file.

Nodes in the male lncRNA-mRNA expression correlation network.

(XLSX)

Click here for additional data file.

Nodes in the female lncRNA-mRNA-GO-pathway network.

(XLSX)

Click here for additional data file.

Nodes in the male lncRNA-mRNA-GO-pathway network.

(XLSX)

Click here for additional data file.

(XLSX)

Click here for additional data file.

(XLSX)

Click here for additional data file.

10 Aug 2020

PONE-D-20-16411

Integrated analysis of long non-coding RNAs and mRNA profiles reveals potential sex-dependent biomarkers of bevacizumab/erlotinib response in advanced lung cancer

PLOS ONE

Dear Dr. Tu,

Thank you for submitting your manuscript to PLOS ONE. After careful consideration, we feel that it has merit but does not fully meet PLOS ONE’s publication criteria as it currently stands. Therefore, we invite you to submit a revised version of the manuscript that addresses the points raised during the review process.

Please submit your revised manuscript by Sep 24 2020 11:59PM. If you will need more time than this to complete your revisions, please reply to this message or contact the journal office at plosone@plos.org. When you're ready to submit your revision, log on to https://www.editorialmanager.com/pone/ and select the 'Submissions Needing Revision' folder to locate your manuscript file.

Please include the following items when submitting your revised manuscript:

A rebuttal letter that responds to each point raised by the academic editor and reviewer(s). You should upload this letter as a separate file labeled 'Response to Reviewers'.

A marked-up copy of your manuscript that highlights changes made to the original version. You should upload this as a separate file labeled 'Revised Manuscript with Track Changes'.

An unmarked version of your revised paper without tracked changes. You should upload this as a separate file labeled 'Manuscript'.

If you would like to make changes to your financial disclosure, please include your updated statement in your cover letter. Guidelines for resubmitting your figure files are available below the reviewer comments at the end of this letter.

If applicable, we recommend that you deposit your laboratory protocols in protocols.io to enhance the reproducibility of your results. Protocols.io assigns your protocol its own identifier (DOI) so that it can be cited independently in the future. For instructions see: http://journals.plos.org/plosone/s/submission-guidelines#loc-laboratory-protocols

We look forward to receiving your revised manuscript.

Kind regards,

Rama Krishna Kancha

Academic Editor

PLOS ONE

Journal Requirements:

When submitting your revision, we need you to address these additional requirements.

1. Please ensure that your manuscript meets PLOS ONE's style requirements, including those for file naming. The PLOS ONE style templates can be found at

https://journals.plos.org/plosone/s/file?id=wjVg/PLOSOne_formatting_sample_main_body.pdf and

https://journals.plos.org/plosone/s/file?id=ba62/PLOSOne_formatting_sample_title_authors_affiliations.pdf

2. Please ensure that you refer to Figure xxxxx in your text as, if accepted, production will need this reference to link the reader to the figure.

Additional Editor Comments (if provided):

The manuscript entitled "Integrated analysis of long non-coding RNAs and mRNA profiles reveals potential sex-dependent biomarkers of bevacizumab/erlotinib response in advanced lung cancer" by Tu et al. demonstrates sex-specific differential expression of biomarkers in relation to drug response in lung cancer. The reviewers are of the opinion that the work is of significant interest. However, the reviewers raised important questions that need to be answered to bring more clarity and highlight the importance of the study.

[Note: HTML markup is below. Please do not edit.]

Reviewers' comments:

Reviewer's Responses to Questions

Comments to the Author

1. Is the manuscript technically sound, and do the data support the conclusions?

The manuscript must describe a technically sound piece of scientific research with data that supports the conclusions. Experiments must have been conducted rigorously, with appropriate controls, replication, and sample sizes. The conclusions must be drawn appropriately based on the data presented.

Reviewer #1: Partly

Reviewer #2: Partly

**********

2. Has the statistical analysis been performed appropriately and rigorously?

Reviewer #1: Yes

Reviewer #2: Yes

**********

3. Have the authors made all data underlying the findings in their manuscript fully available?

The PLOS Data policy requires authors to make all data underlying the findings described in their manuscript fully available without restriction, with rare exception (please refer to the Data Availability Statement in the manuscript PDF file). The data should be provided as part of the manuscript or its supporting information, or deposited to a public repository. For example, in addition to summary statistics, the data points behind means, medians and variance measures should be available. If there are restrictions on publicly sharing data—e.g. participant privacy or use of data from a third party—those must be specified.

Reviewer #1: Yes

Reviewer #2: Yes

**********

4. Is the manuscript presented in an intelligible fashion and written in standard English?

PLOS ONE does not copyedit accepted manuscripts, so the language in submitted articles must be clear, correct, and unambiguous. Any typographical or grammatical errors should be corrected at revision, so please note any specific errors here.

Reviewer #1: Yes

Reviewer #2: Yes

**********

5. Review Comments to the Author

Please use the space provided to explain your answers to the questions above. You may also include additional comments for the author, including concerns about dual publication, research ethics, or publication ethics. (Please upload your review as an attachment if it exceeds 20,000 characters)

Reviewer #1: The manuscript sounds 1. non coding genes transcriptional control is elusive and the work provides expression differences by sex dependent markers and 2. lung cancer dependent alterations with noncoding and mRNA profiles by a clinical drug perturbation is elaborated. Therefore, the manuscript is interest to the literature and I recommend the ms.

Reviewer #2: Gist/Summary: The authors come up with identification of several differentially expressed (DE) lncRNAs and mRNAs from publicly available gene expression omnibus (GEO) datasets associated with NSCLC with bevacizumab and erlotinib (BE) treatment. Then they make enrichment/pathway analyses based on gender and identify unique DEs

Strengths: Novel in what they claim that the authors mention on the the first of its kind of work on BE treatment ( which invariably is not so!)

Figures

Weaknesses: The work is based on GEO datasets and not original

Hazy texts in between

Major comments:

Wei et al. 2016 and Saito eta l. 2019, 2020 indeed mention this and identify DE genes ( if NOT lncRNAs) in BE treatment. The authors must crosscheck their works and cite them.

Lines 95-99: The authors must glorify lncRNAs and on why they could be used as biomarkers

Some of the texts go hand-i-hand with clinicaltrials.gov. The authors must check it

Lines 130 in Methods must be renamed. The methodology is out of vigour. As the work is public interpretation of GEO., I wonder why the authors mention about the ethics statement etc., which gives a falsehood to the story.

There must be a methodological flowchart

Have the authors carried out the bidirectional best blast hits whence inferring lncRNAs from BLAST?

The authors must mention whether or not the exon array enrichment they infer uniquely mapped to genes and the lncRNAs were indeed inferred from those genes.

Line 209 on interactions. The combined score of 0.4 doesn't make sense

The NONHSAT095695.2 appears to be largely expressed in liver/adrenal glands NOT lungs. How do the authors justify this?

The whole purpose of making a story for DEs for male and females is not warranted and justified in haste. For example, why males have SLE related DEs is not discussed. The same applies to the coalesce of pathways between males and females.

The multivenn diagram with union of intersection of venns showing 1 and 2 respectively for DE lncRNAs and De mRNAs were not discussed in great detail. For example, where are they localized to and their expression patterns

The conclusions and discussions must herald have point son what if multi-omic studies are integrated to it or say, RNA-Seq expression data is used instead of exon array

Minor comments

The abstract must have an expansion of EGFR

Lines 385 and line 296: Please rewrite or use comma

All figures must be high resolution

**********

6. PLOS authors have the option to publish the peer review history of their article (what does this mean?). If published, this will include your full peer review and any attached files.

If you choose “no”, your identity will remain anonymous but your review may still be made public.

Do you want your identity to be public for this peer review? For information about this choice, including consent withdrawal, please see our Privacy Policy.

Reviewer #1: Yes: Yusuf TUTAR

Reviewer #2: Yes: Prashanth Suravajhala

[NOTE: If reviewer comments were submitted as an attachment file, they will be attached to this email and accessible via the submission site. Please log into your account, locate the manuscript record, and check for the action link "View Attachments". If this link does not appear, there are no attachment files.]

While revising your submission, please upload your figure files to the Preflight Analysis and Conversion Engine (PACE) digital diagnostic tool, https://pacev2.apexcovantage.com/. PACE helps ensure that figures meet PLOS requirements. To use PACE, you must first register as a user. Registration is free. Then, login and navigate to the UPLOAD tab, where you will find detailed instructions on how to use the tool. If you encounter any issues or have any questions when using PACE, please email PLOS at figures@plos.org. Please note that Supporting Information files do not need this step.

21 Sep 2020

Dear Editor,

We have now revised our manuscript PONE-D-20-16411 taking into account the useful comments of the academic editor and the two reviewers.

We thank you and the reviewers for your careful reading and the suggestions to improve our initial manuscript.

We took some time to respond to each point raised by the academic editor and reviewers, the detailed responses to the reviewers’ comments are explained below. And we also uploaded our figure files to the Preflight Analysis and Conversion Engine (PACE) digital diagnostic tool (https://pacev2.apexcovantage.com/) to ensure that all of our figures except Supporting Information files meet PLOS requirements. And we checked all themanuscript to make sure that it meets PLOS ONE’s styl requirements, including those for file naming according to the PLOS ONE style templates.

We are confident that these additional modifications strenthened and ameliorated this manuscript.

In addition, the author information was modified according to the contribution of all the authors after the revision.

We hope that the reviewers and yourself will find our manuscript suitable for publication in PLOS ONE.

Best regards,

Chao TU

The detailed responses to the reviewers’ comments are explained below:

Wei et al. 2016 and Saito et al. 2019, 2020 indeed mention this and identify DE genes ( if NOT lncRNAs) in BE treatment. The authors must crosscheck their works and cite them.

We thank the reviewer for providing the references regarding BE treatment. We carefully searched these references, and found that Wei et al. 2018 (PMID:30032815) and Saito et al.2019 (PMID:30975627) did mention about the treatment of erlotinib plus bevacizumab(BE) in non-small cell lung cancer (NSCLC). Both papers investigated the efficacy and safety of erlotinib plus bevacizumab treatment in patients with non-small-cell lung cancer. Instead, the reference of Saito et al.2020 studied the effectiveness and safety of nivolumab treatment, not BE treatment, in NSCLC patients in real-world setting in Japan. However, none of these papers did the study on the profiling of differentially expressed genes or non-coding RNAs after BE treatment in NSCLC patients. Nevertheless, we discussed these papers in the manuscript and cited them as references according to the reviewer’s suggestion.

Lines 95-99: The authors must glorify lncRNAs and on why they could be used as biomarkers.

We agree with the reviewer that we should introduce more about lncRNAs and why they could be used as biomarkers.

LncRNA is a class of non-coding RNAs of more than 200 nucleotides with biological functions. They are widely involved in biological development and pathological process through chromatin remodelling, cis-regulation at enhancers and post transcriptional regulation. LncRNAs have been found to be specifically expressed in the oncogenesis and metastasis of many tumor tissues (June 2014, Nature Genetics, kiyer). Therefore, these lnRNAs are also considered to be associated with carcinogenesis or tumor inhibition, and are potential biomarkers that can predict the tissue carcinogenesis.

And we added the above introduction concerning why lncRNAs could be used as biomarkers in the manuscript.

Some of the texts go hand-i-hand with clinicaltrials.gov. The authors must check it.

We are sorry that we are not so sure about what the reviewer mentions in this point. We have written the text independent of clinicaltrials.gov. And we did the duplicate checking for our manuscript, and no similar content to clinicaltrials.gov was found.

Lines 130 in Methods must be renamed. The methodology is out of vigour. As the work is public interpretation of GEO., I wonder why the authors mention about the ethics statement etc., which gives a falsehood to the story.

We are sorry for the redundancy of Line 130 in Methods, and it has been renamed to “Data sets”. We agree with the reviewer that it’s not necessary to mention about the ethics statement in Methods as our work is based on the public GEO dataset, therefore, we removed those concerning ethics statement and written informed consent in this part (Please refer to the end of “Data sets” in Materials and Methods).

There must be a methodological flowchart

We have added the flowchart in “Supplementary Figure 1.tif” for the sample selection and data analysis used in the manuscript.

Have the authors carried out the bidirectional best blast hits whence inferring lncRNAs from BLAST?

Yes, we have carried out blast with the bidirectional best hit information method, which is more precise and possesses higher matching degree compared with the single-directional best hit information method.

The authors must mention whether or not the exon array enrichment they infer uniquely mapped to genes and the lncRNAs were indeed inferred from those genes.

We have removed those probes that were not uniquely mapped to genes, therefore, the exon array enrichment we used was uniquely mapped to genes and the lncRNAs were indeed inferred from those genes. We have added this information in Materials and Methods (Please refer to the part of “mRNA and lncRNA profiling” in Materials and Methods).

Line 209 on interactions. The combined score of 0.4 doesn't make sense

We thank the reviewer for this comment. However, we are not so sure about what the reviewer means here. When performing protein-protein interaction (PPI) network analysis, we used the interaction score of 0.4 according to the recommendation of the official developping team for STRING database in this paper from Nucleic Acids Research, 2017:

The STRING database in 2017: quality-controlled protein–protein association networks, made broadly accessible.

And this confidence cutoff of 0.4 is used widely by many other groups, below are some example references that used the interaction score of 0.4 when PPI network:

①Deciphering miRNA transcription factor feed-forward loops to identify drug repurposing candidates for cystic fibrosis. Genome Medicine 2014, 6:94

②Identifification of key pathways and genes in colorectal cancer using bioinformatics analysis. Med Oncol (2016) 33:111. DOI 10.1007/s12032-016-0829-6

③Screening and identification of key biomarkers in hepatocellular carcinoma: Evidence from bioinformatic analysis.

The NONHSAT095695.2 appears to be largely expressed in liver/adrenal glands NOT lungs. How do the authors justify this?

We thank the reviewer for this comment. In our present study, we only detected the differential expression of NONHSAT095695.2 through bioinformatic analysis by gene annotation, this needs to be validated by experiments. In our future work, we will detect the expression level of NONHSAT095695.2 in lungs of both male and female patients before and after BE treatment by RT-QPCR.

The whole purpose of making a story for DEs for male and females is not warranted and justified in haste. For example, why males have SLE related DEs is not discussed. The same applies to the coalesce of pathways between males and females.

We thank the reviewer for this point. Gender differences play an important role in the drug response in cancer, especially in lung cancer, for example, in the following papers they all focus their work on sex-dependent drug response in lung cancer:

①Conforti F, Pala L, Bagnardi V, et al. Sex-Based Heterogeneity in Response to Lung Cancer Immunotherapy: A Systematic Review and Meta-Analysis. J Natl Cancer Inst. 2019;111(8):772-781. doi:10.1093/jnci/djz094.

②Wang S, Zhang J, He Z, Wu K, Liu XS. The predictive power of tumor mutational burden in lung cancer immunotherapy response is influenced by patients' sex. Int J Cancer. 2019;145(10):2840-2849. doi:10.1002/ijc.32327.

③Tsiouda T, Sardeli C, Porpodis K, et al. Sex Differences and Adverse Effects between Chemotherapy and Immunotherapy for Non-Small Cell Lung Cancer. J Cancer. 2020;11(11):3407-3415. Published 2020 Mar 5. doi:10.7150/jca.40196.

We added the discussion about the coalesce of pathways between males and females in the manuscript. However, we didn't mention SLE-related DEs because no research in our lab is involved in the area of SLE. Sorry for this and thank you for your comprehension.

The multivenn diagram with union of intersection of venns showing 1 and 2 respectively for DE lncRNAs and De mRNAs were not discussed in great detail. For example, where are they localized to and their expression patterns

Thank you very much for pointing out this issue. We have checked the venn diagram in Fig 5A and Fig 5B, respectively, and found that the union intersection of the significantly enriched two GO terms overlapping the female-upregulated, male-upregulated, female-downregulated and male-downregulated DE mRNA datasets in Fig 5A are : “positive regulation of transcription from RNA polymerase II promoter ”and “negative regulation of transcription from RNA polymerase II promoter”, while the union intersection of the significantly enriched one shared pathway among those four groups in Fig 5B is : “Metabolic pathways”. The detailed mRNAs of the localization in the union intersection of shared GO terms and pathway mentioned above and their expression patterns are provided in Supplementary Table 7 and Supplementary Table 8, respectively. And we have added the discussion about the union of intersection of venns showing 1 and 2 respectively for shared GO functions and pathways in the manuscript according to the reviewer’s suggestion.

The conclusions and discussions must herald have point son what if multi-omic studies are integrated to it or say, RNA-Seq expression data is used instead of exon array

We thank the reviewer for this comment. It is sure that more potential differentially expressed long non-coding RNAs and coding genes can be discovered if RNA-seq data is used instead of exon array data. In future, we will integrate more multi-omic studies to our work in order to discover more possible lncRNA-mRNA correlation. We have added these points in the Discussion and Conclusion part. Please refer to the manuscript for the modifications.

Minor comments

The abstract must have an expansion of EGFR

We have added the full name of EGFR : “Epidermal growth factor receptor” in the abstract. Please refer to the abstract in the manuscript for the modification.

Lines 385 and line 296: Please rewrite or use comma

We have rewritten line 296 and line 385, please refer to the manuscript for the modifications.

All figures must be high resolution.

We have adjusted the figures with low resolution according to the PLOS ONE guidelines: https://journals.plos.org/plosone/s/figures. Now the resolution of all the submitted figures is at least 600dpi. We hope that the newly submitted figures will meet the requirements of PLOS ONE.

Submitted filename: Response to Reviewers.docx

Click here for additional data file.

30 Sep 2020

I ntegrated analysisof long non-coding RNAs and mRNA profiles reveal s  potential sex-dependent biomarkers of bevacizumab/erlotinib response in advanced lung cancer

PONE-D-20-16411R1

Dear Dr. Tu,

We’re pleased to inform you that your manuscript has been judged scientifically suitable for publication and will be formally accepted for publication once it meets all outstanding technical requirements.

Within one week, you’ll receive an e-mail detailing the required amendments. When these have been addressed, you’ll receive a formal acceptance letter and your manuscript will be scheduled for publication.

An invoice for payment will follow shortly after the formal acceptance. To ensure an efficient process, please log into Editorial Manager at http://www.editorialmanager.com/pone/, click the 'Update My Information' link at the top of the page, and double check that your user information is up-to-date. If you have any billing related questions, please contact our Author Billing department directly at authorbilling@plos.org.

If your institution or institutions have a press office, please notify them about your upcoming paper to help maximize its impact. If they’ll be preparing press materials, please inform our press team as soon as possible -- no later than 48 hours after receiving the formal acceptance. Your manuscript will remain under strict press embargo until 2 pm Eastern Time on the date of publication. For more information, please contact onepress@plos.org.

Kind regards,

Rama Krishna Kancha

Academic Editor

PLOS ONE

Additional Editor Comments (optional):

Reviewers' comments:

2 Oct 2020

PONE-D-20-16411R1

Integrated analysis of long non-coding RNAs and mRNA profiles reveals potential sex-dependent biomarkers of bevacizumab/erlotinib response in advanced lung cancer

Dear Dr. Tu:

I'm pleased to inform you that your manuscript has been deemed suitable for publication in PLOS ONE. Congratulations! Your manuscript is now with our production department.

If your institution or institutions have a press office, please let them know about your upcoming paper now to help maximize its impact. If they'll be preparing press materials, please inform our press team within the next 48 hours. Your manuscript will remain under strict press embargo until 2 pm Eastern Time on the date of publication. For more information please contact onepress@plos.org.

If we can help with anything else, please email us at plosone@plos.org.

Thank you for submitting your work to PLOS ONE and supporting open access.

Kind regards,

PLOS ONE Editorial Office Staff

on behalf of

Dr. Rama Krishna Kancha

Academic Editor

PLOS ONE