Identification of the targets of hematoporphyrin derivative in lung adenocarcinoma using integrated network analysis

Hematoporphyrin derivative (HPD) has a sensibilization effect in lung adenocarcinoma. This study was conducted to identify the target genes of HPD in lung adenocarcinoma. RNA sequencing was performed using the lung adenocarcinoma cell line A549 after no treatment or treatment with X-ray or X-ray + HPD. The differentially expressed genes (DEGs) were screened using Mfuzz package by noise-robust soft clustering analysis. Enrichment analysis was carried out using “BioCloud” online tool. Protein–protein interaction (PPI) network and module analyses were performed using Cytoscape software. Using WebGestalt tool and integrated transcription factor platform (ITFP), microRNA target and transcription factor (TF) target pairs were separately predicted. An integrated regulatory network was visualized with Cytoscape software. A total of 815 DEGs in the gene set G1 (continuously dysregulated genes along with changes in processing conditions [untreated—treated with X-ray—X-ray + treated with HPD]) and 464 DEGs in the gene set G2 (significantly dysregulated between X-ray + HPD-treated group and untreated/X-ray-treated group) were screened. The significant module identified from the PPI network for gene set G1 showed that ribosomal protein L3 (RPL3) gene could interact with heat shock protein 90 kDa alpha, class A member 1 (HSP90AA1). TFs AAA domain containing 2 (ATAD2) and protein inhibitor of activated STAT 1 (PIAS1) were separately predicted for the genes in gene set G1 and G2, respectively. In the integrated network for gene set G2, ubiquitin-specific peptidase 25 (USP25) was targeted by miR-200b, miR-200c, and miR-429. RPL3, HSP90AA1, ATAD2, and PIAS1 as well as USP25, which is targeted by miR-200b, miR-200c, and miR-429, may be the potential targets of HPD in lung adenocarcinoma.

types of lung cancers, it is likely to show distant metastasis at an early stage [6]. Therefore, evaluation of the mechanisms underlying lung adenocarcinoma progression is desirable to improve treatment outcomes.
Porphyrins may selective accumulate in malignant tumors [7], and hematoporphyrin derivative (HPD) has been used as a photosensitizer in the photodynamic therapy (PDT) of lung cancer [8]. It was demonstrated that 99mTc-hematoporphyrin-linked albumin nanoparticles (99mTc-HP-ANPs) may be used for the PDT and radiodiagnosis of lung cancer [9]. The combination of PDT and photosensitizers HPD and 5-aminolevulinic acid (ALA) may increase the curative rate for skin cancers, distinctly cut down the photosensitive period, and reduce the dose of pro-toxic HPD [10]. Hematoporphyrin conjugated with the antibodies directed to vascular endothelial growth factor shows high antitumor activities in patients with Lewis lung carcinoma [11]. The expression level of interleukin-6 affects the cellular sensitivity to PDT, and the combination of PDT and interleukin-6 may serve as a novel strategy for the therapy of Lewis lung carcinoma [12]. Wang et al. used X-ray as the energy source for PDT activation and suggested that the X-ray-induced photodynamic therapy (X-PDT) may be used as a novel therapeutic method against human cancers [13]. Austerlitz et al. found that the response of the Fricke solution doped with hematoporphyrin and irradiated with lowenergy X-rays was enhanced in PDT [14]. However, the mechanism underlying HPD effects in lung adenocarcinoma are not investigated.
As a commonly studied lung adenocarcinoma cell line, A549 was used as a model in the present study. RNA sequencing was applied to the untreated A549 cells as well as those treated with X-ray or the combination of X-ray and HPD. The differentially expressed genes (DEGs) were screened and enrichment analysis, proteinprotein interaction (PPI) network, and module analyses as well as integrated network analysis were carried out to identify the important genes affected by HPD.

Cell cultivation
The lung adenocarcinoma cell line A549 was acquired from the Cell Bank of Chinese Academy of Sciences. The cells were cultured in Dulbecco's modified Eagle's medium (DMEM) (GIBCO) supplemented with 1% penicillin/streptomycin double-antibody (GIBCO) and 10% fetal bovine serum (FBS, GIBCO) in a humidified 5% CO 2 incubator (Thermo) at 37 °C. After the medium was discarded, the cells were washed once with cold phosphate-buffered saline (PBS) and treated with pancreatin (GIBCO) at 37 °C for 2 min. The suspension was treated with complete medium to neutralize pancreatin, and the mixture was centrifuged (1000 rpm, 5 min). The supernatant was discarded and the cells were resuspended in FBS-supplemented DMEM and cultured in a humidified 5% CO 2 incubator (Thermo) at 37 °C.

Flow cytometry assay
To detect the apoptosis of A549 cells, flow cytometry assay was conducted according to the previously described method [15]. Cells were counted and cultured into six-well plates (ABI, 2 × 10 5 cells/well) in a humidified 5% CO 2 incubator (Thermo) at 37 °C overnight. The cells were subsequently treated with different concentrations of HPD and different doses of X-ray. After the medium was discarded, the cells were digested with pancreatin (GIBCO), followed by treatment with fresh medium to deactivate pancreatin and centrifugation (1000 rpm, 6 min). The supernatant was discarded and the cells were washed once with PBS and resuspended in 1× binding buffer (BD Biosciences; 400 μL 1× binding buffer for the control group and 100 μL 1× binding buffer for other groups). A total of 100 μL of the above solution was transferred into flow tubes and treated with 5 μL of fluorescein isothiocyanate (FITC)-Annexin V (BD Biosciences) and 5 μL of propidium iodide (PI, BD Biosciences, 50 μg/mL) (the  control group was divided into unstained, Annexin V-stained, PI-stained, and Annexin V + PI-stained groups). After being incubated (in the absence of light) for 15 min at room temperature, the cells were treated with 400 μL of 1× binding buffer (BD Biosciences) and analyzed with a flow cytometer (BD Biosciences).

RNA extraction and RNA-seq library construction
The cells were counted, seeded into 100 mm culture dishes (1 × 10 6 cells/well), and cultured in a humidified 5% CO 2 incubator (Thermo) at 37 °C overnight. After treatment with different concentrations of HPD or different doses of X-ray (the control group was left untreated; the X-ray group was treated with 10 Gy X-ray and cultured for 24 h; X-ray + HPD group was treated with 10 Gy X-ray + 10 μg/mL of HPD and cultured for 24 h; each group had three replicates), the cells were washed twice with cold PBS. Total RNA was extracted using Trizol reagent (TaKaRa) following the manufacturer's instruction and quantified with a spectrophotometer (Nanodrop). RNA-seq library was constructed with NEBNext ® Ultra ™ RNA Library Prep Kit for Illumina ® (New England Biolabs) and sequencing was performed on Illumina Hiseq 4000 (PE150) (Illumina). The sequencing data were deposited into the Sequence Read Archive (SRA) database under the accession number of SRP091521.

Data preprocessing and DEG screening
The Prinseq-lite (http://edwar ds.sdsu.edu/cgi-bin/prins eq/prins eq.cgi) tool [16] and FASTX_ToolKit (http:// hanno nlab.cshl.edu/fastx _toolk it/) [17] were applied for the quality control of the raw data. Barcode and adaptor sequences in the reads were removed. The reads with N content larger than 5% were filtered out. Bases with continuous quality under 10 at 5′ or 3′ end were discarded. The reads with low quality (having over 20% bases with quality lower than 20) and those shorter than 30 nt were removed. The clean reads obtained from the three groups of samples were mapped to GRCH38 human genome using TopHat software (version 2.0.8) [18]. The fragments per kilobase million (FPKM) and read count matrix of the genes were acquired using StringTie tool [19], and gene annotation information was obtained from GENCODE database (version 24, http://genom e.imim.es/genco de/) [20]. To identify DEGs, noise-robust soft clustering analysis was performed using the fuzzy c-means clustering algorithm in Mfuzz package (http://www.bioco nduct or.org/packa ges/relea se/bioc/html/Mfuzz .html) [21]. Specific clusters were selected based on gene expression trends. Both minSTD and acore parameters were set as 0.5.

Functional and pathway enrichment analysis
Gene ontology (GO, http://www.geneo ntolo gy.org) database can use structured vocabularies for noting genes or gene products from three aspects (MF, molecular function; BP, biological process; and CC, cellular component) [22]. The Kyoto Encyclopedia of Genes and Genomes (KEGG, http://www.genom e.ad.jp/kegg) database links genomic information with functional information through the investigation of gene functions [23]. "BioCloud" online tool (http://www.biocl oudse rvice .com) is developed for settling computing problems of high-throughput biological data. Using "BioCloud" online tool, DEGs were subjected to GO functional and KEGG pathway enrichment analyses with the threshold of p-value < 0.05.

PPI network and module analyses
Search Tool for the Retrieval of Interacting Genes (STRING, http://strin g-db.org/) is a database that collects the PPIs involving more than 1100 organisms [24]. Based on STRING database [24], the PPIs among the proteins corresponding to DEGs were analyzed with a combined score > 0.4 as the cut-off criterion. PPI network was subsequently visualized using Cytoscape software (http://www.cytos cape.org) [25], and the hub nodes [26] in the PPI network were screened by calculating their connectivity degrees. Module analysis for PPI network was conducted using the Molecular Complex Detection (MCODE) plugin [27] of Cytoscape software. In addition, enrichment analysis was performed for the nodes of significant modules using "BioCloud" online tool.

Integrated network analysis
WEB-based gene set analysis toolkit (WebGestalt, http:// www.webge stalt .org) [28] was used to predict the genes involved in the PPI network at p < 0.001 and the number Fig. 4 The protein-protein interaction network constructed for the genes in gene set G1. Red and green circles represent the upregulated and downregulated genes, respectively of target genes ≥ 4 as the thresholds. Using Cytoscape software [25], the microRNA (miRNA) target regulatory network was constructed. According to the integrated transcription factor platform (ITFP, http://itfp.biosi no.org/itfp) [29], the transcription factors (TFs) targeting DEGs and the differentially expressed TFs were predicted. The TF target regulatory network was visualized with Cytoscape software [25]. The PPI network, miRNA target regulatory network, and TF target regulatory network were integrated, and an integrated network was constructed with Cytoscape software [25].

Statistical analysis
One-way analysis of variance and two-tailed t-test were applied for statistical analysis using GraphPad prism software (GraphPad Software, San Diego, CA). Data were shown as the mean ± standard error of the mean (SEM). A value of p < 0.05 was considered statistically significant.

Effects of HPD and X-ray on the proliferation and apoptosis of A549 cells
A549 cells were treated with different concentrations of HPD and different doses of X-ray and their proliferation rate was calculated. HPD in combination with X-ray significantly suppressed the proliferation rate of A549 cells (p < 0.05, Fig. 1a). In particular, 10 μg/mL of HPD and 10 Gy X-ray was the lowest concentration-dose combination that showed significant effects on A549 cell proliferation (p < 0.001). A549 cells were treated with 10 μg/ mL of HPD and 10 Gy X-ray, and their apoptosis rate was analyzed with flow cytometry. The results of flow cytometry analysis showed that the combination of HPD (10 μg/mL) and X-ray (10 Gy) significantly increased the apoptosis of A549 cells (p < 0.05, Fig. 1b). Therefore, the combination of 10 μg/mL of HPD and 10 Gy X-ray was used in the subsequent experiments.

Data preprocessing and DEG analysis
Data sequencing was carried out with quality control (Table 1), and the sequences were mapped to GRCH38 human genome ( Table 2). The gene expression matrix was processed with Mfuzz package to reveal a total of 14 clusters (Fig. 2). According to the experimental design, only two types of clusters were selected for analysis. One type of clusters showed continuous upregulation (cluster 2 and 3) or downregulation (cluster 7 and 14) of gene expression along with the change in processing conditions (untreated-treated with X-ray-treated with X-ray + HPD) (containing a total of 815 genes that were included in gene set G1). Another type of cluster included the significantly upregulated (cluster 13) or downregulated (cluster 11) genes under the processing condition of X-ray + HPD in comparison with the processing condition of untreated and X-ray treatment (containing a total of 464 genes that were included in gene set G2).

Functional and pathway enrichment analysis
The genes in gene set G1 and G2 were separately subjected to functional and pathway enrichment analyzes. The top five enriched terms in BP, CC, MF, and KEGG categories are shown in Fig. 3. The genes in gene set G1 were mainly enriched in double-strand break repair (BP), membrane (CC), ATP binding (MF), and proteoglycans in cancer (KEGG) (Fig. 3a). The genes in gene set G2 were mainly enriched in positive regulation of transcription from RNA polymerase II promoter (BP), nucleus (CC), zinc ion binding (MF), and estrogen signaling pathway (KEGG) (Fig. 3b).

PPI network and module analyses
The PPI network constructed for the genes in gene set G1 had 210 nodes and 333 interactions (Fig. 4). On the other hand, the PPI network constructed for the genes in gene set G2 had 135 nodes and 164 interactions (Fig. 5). The top 10 nodes with high degrees in PPI networks are listed in Table 3 and included heat shock protein 90 kDa alpha, class A member 1 (HSP90AA1, degree = 31) and ribosomal protein L3 (RPL3, degree = 10). Based on MCODE plugin, one significant module was identified from the PPI network constructed for the genes in gene set G1 that included 13 nodes and 37 interactions (such as HSP90AA1-RPL3) (Fig. 6). The terms enriched for the genes in the module are listed in Table 4, and mainly included translational initiation (BP, p-value = translational initiation), cytosol (CC, p-value = 1.28E−09), poly(A) RNA binding (MF, p-value = 1.70E−06), and ribosome (KEGG, p-value = 4.57E−06). No significant module was identified from the PPI network constructed for the genes in gene set G2.

Integrated network analysis
The miRNAs of the genes implicated in the PPI networks constructed for the genes in gene set G1 (Table 5) and G2 (Table 6) were predicted. The TFs targeting the genes in gene set G1 (ATPase family, AAA domain containing 2 [ATAD2]) and G2 (protein inhibitor of activated STAT 1 [PIAS1]) were also analyzed ( Table 7). An integrated network was constructed for the genes in gene set G1 that had 259 nodes (including 25 TFs and 31 miRNAs) and 687 pairs (Fig. 7). The Table 3 The top 10 nodes in the protein-protein interaction networks constructed for the genes in gene set G1 and the genes in gene set G2    Table 5 The miRNAs targeted the genes implicated in the protein-protein interaction networks constructed for the genes in gene set G1  Table 8.

Discussion
In this study, the CCK-8 assay result revealed the significant suppression in the proliferation activity of A549 cells in response to the treatment with the combination of HPD and X-ray. The combination of 10 μg/mL of HPD and 10 Gy X-ray was selected as the lowest concentration/dose that achieved a significant increase in the apoptosis of A549 cells, which might be a limitation in terms of analysis of the data. In Europe, hematoporphyrin is the most commonly used photosensitizer for the treatment of advanced lung cancer [30,31]. Previous studies have demonstrated that HPD-PDT may inhibit proliferation and induce apoptosis of A549 cells, thereby inducing effective killing of adenocarcinoma cells [32,33]. After the optimization of the combination treatment of HPD and X-ray, a series of bioinformatic analyses were performed with the RNA-seq data. Through noise-robust soft clustering analysis, 815 genes that showed continuous upregulated or downregulated expression along with the change in processing conditions (untreated-treated with X-ray-treated Table 6 The miRNAs targeted the genes implicated in the protein-protein interaction networks constructed for the genes in gene set G2 with X-ray + HPD) were included in the gene set G1. A total of 464 genes that were significantly upregulated or downregulated under the processing condition of X-ray + HPD were included in the gene set G2. The significant module identified from the PPI network constructed for the genes in gene set G1 revealed the interaction between RPL3 and HSP90AA1. RPL3 functions in the response of cells to oxaliplatin-and 5-fluorouracil-induced nucleolar stress and may be used to improve the therapeutic effects of these drugs against cancers [34]. The chemotherapy curative effect of actinomycin D is determined by RPL3 status in cancers shorting of p53; thus, high level of RPL3 may be useful for the treatment of lung and colon cancers [35].
The frequencies of mutant genotypes of HSP90AA1, HSP90AB1, and HSP90B1 are reported to be significantly higher in the patients with non-small cell lung cancer (NSCLC) in the Turkish population [36]. Downregulation of HSP90 expression correlated with increased overall survival of patients with NSCLC, and HSP90 inhibitor exerts an antiproliferative effect on NSCLC cell lines [37,38]. These observations suggest that RPL3 interacting with HSP90AA1 may be associated with the sensibilization effect of HPD in lung adenocarcinoma. ATAD2 and PIAS1 were separately predicted as the TFs targeting the genes from the gene sets G1 and G2. Caron et al. demonstrated that ATAD2 overexpression Fig. 7 The integrated network constructed for the genes in gene set G1. Red and green separately indicate the upregulated and downregulated genes, respectively. Circles and triangles represent genes and transcription factors, respectively. Blue diamonds indicate miRNAs may promote the malignant transformation of lung and breast cancers by affecting the basic properties of chromatin [39]. Wang et al. found that ATAD2/AAA + nuclear coregulatory cancer associated (ANCCA ) may serve as a promising biomarker for the treatment and prognosis of squamous cell lung carcinoma [40]. PIAS1 contributes to cytoplasm-nuclear distribution of focal adhesion kinase by interacting with it, and focal adhesion kinase activity in the nucleus facilitates survival and progression of NSCLC via promotion of DNA repair regulation and cellextracellular matrix interaction [41,42]. PIAS1 mediates oncogenic signaling by promoting promyelocytic leukemia (PML) degradation, and PIAS1 and PML expression is negatively correlated in NSCLC cell lines [43]. Therefore, ATAD2 and PIAS1 may be involved in the action mechanism of HPD in lung adenocarcinoma.
In the integrated network for the genes in the gene set G2, USP25 was targeted by miR-200b, miR-200c, and miR-429. miR-200c may serve as a tumor suppressor in NSCLC through the inhibition of USP25 expression and may be applied for therapeutic purposes [44]. The overexpression of miR-200c and miR-141 is associated with the short overall survival of patients with lung adenocarcinoma via angiogenesis and mesenchymal-epithelial transition [45]. The low expression of miR-200b is reported to induce E2F transcription factor 3 overexpression and increase the chemoresistance of patients with lung adenocarcinoma to docetaxel [46]. Zhu et al. suggested that the serum levels of miR-29c and miR-429 may be used as non-invasive biomarkers for patients with early stage NSCLC [47]. Lang et al. found that miR-429 contributes to cell proliferation and metastasis and regulates several tumor suppressor genes in patients with NSCLC, serving as a possible therapeutic target [48]. These observations suggest that USP25 targeted by miR-200b, miR-200c, and miR-429 may also function in the action process of HPD in lung adenocarcinoma.

Conclusion
A total of 815 DEGs in gene set G1 were identified along with a change in processing conditions (untreated-treated with X-ray-treated with X-ray + HPD). In addition, 464 DEGs in gene set G2 were screened under the processing condition of X-ray + HPD. RPL3, HSP90AA1, ATAD2, as well as PIAS1 and USP25, which is targeted by miR-200b, miR-200c, and miR-429 may show correlations with the sensibilization effect of HPD in lung adenocarcinoma. Further validation with experimental research is warranted to confirm the roles of these genes in the sensibilization effect of HPD in lung adenocarcinoma.