Transcriptome Sequencing Reveals Gene Expression Differences in the Response of Malignant Melanomas Cells to Hypoxia

Background Malignant melanoma (MM) is the most deadly type of skin cancer, with 5-year survival rate of less than 16%. HIF-1α overexpression is associated with poor prognosis in many cancers including MM. Hence, we characterized differentially expressed genes (DEGs) in the response of MM cells to normal and hypoxia. Methods We rst successfully constructed cell hypoxia model and then performed RNA-seq to explore the changes of gene transcription in MM cells during hypoxia. The highest expression of the six genes was detected using qRT-PCR and western blot assays. We explored the binding sites between BIRC7 promoter and HIF-1α by dual-luciferase assay. Cellular function assays were used to observe the role of BIRC7 in the effect of hypoxia on tumor progression. Results We found that compared with the transcriptome data of the control group, a total of 2601 DEGs were identied in the hypoxic group. There were 1517 genes with signicantly higher expression and 1084 genes with lower expression in the hypoxic group. Among them, OSCAR, BIRC7, HBA1, SFN, GOLT1A, and BEX2 were signicantly up-regulated in the hypoxic group. BIRC7 expression was most signicantly up-regulated and a downstream factor of HIF-1α. We highlighted that knockdown of BIRC7 reversed the positive effects of HIF-1α on A875 and M14 cells. Our ndings demonstrated that BIRC7 was a downstream factor of HIF-1α and reversed the effect of hypoxia on promoting tumor progression of MM cells. Sequencing. We found that signicant differences in DEG between normal and Hypoxic-treated MM cells.


Introduction
Malignant melanoma (MM) is the most deadly type of skin cancer, the incidence of MM accounts for 10% of all skin cancer, and its mortality accounts for more than 80% of skin cancer (Faraj et al., 2020). The 5-year survival rate of patients with metastatic MM is less than 16% (Anton et al., 2018). Studies have shown that the skin of humans and mice are mildly hypoxic, and the oxygen content in the skin is su cient to stabilize HIF-1α and express hypoxic markers, such as carbonic anhydrase IX (CAIX) and glucose transporter 1 (Glut-1) (Johan et al., 2011). Recent research also con rms this claim through quantitative evaluation of skin oxygenation. The study nds that pO 2 in the dermis is 10%, epidermis is 10-0.5%, and the oxygen tension in the hair follicle is between 2.5% and 0.1% (Spence and Beck, 2005). Therefore, the melanocytes present in the dermal-epidermal junction of humans and some transgenic mice or in the hair follicles of mice are under anoxic environment (Gaustad et al., 2017). Hypoxia stabilizes the expression of HIF-1α, which is characteristic of many tumors (Ernestina et al., 2013, Guan et al., 2015, Kim et al., 2011, Nakayama et al., 2002, Sasabe et al., 2005. The emergence of next-generation sequencing (NGS) technology and the development of bioinformatics provide a useful platform for exploring Differently Expressed Genes (DEG) (Zhigalova et al., 2015). In addition, NGS technology can simultaneously assess the mRNA transcription patterns of all genes in various species. RNA sequence analysis (RNA-Seq) using NGS technology has been widely used for high-throughput studies of gene expression (Hu et al., 2016, Nueda et al., 2018. RNA-Seq utilizes high-throughput sequencing methods, and then aligns short sequence reads to the reference genome to analyze complementary DNA (cDNA). This new technology makes it possible to 2. Materials And Methods 2.1 Cell culture and treatment. M14 and A875 cells were obtained from the Cell Bank of the Chinese Academy of Sciences (Shanghai, China). A chamber lled with 95% N 2 and 5% CO 2 simulates the hypoxic environment in melanoma cells. Cells were culture in anoxic chamber for 0 h, 6 h, 12 h, 24 h, and 48 h, respectively. The cells were divided into the following groups: negative control (NC) group, hypoxia group (M14 cells with hypoxia treatment for 12 h), hypoxia + HIF-1α inhibitor group (M14 cells transfected with HIF-1α inhibitor in hypoxia condition for 12 h), siBIRC7 group (M14 cells transfected with siRNA-BIRC7), HIF-1α-OE group (M14 cells transfected with HIF-1α overexpression vector), siBIRC7 + HIF-1α-OE group (M14 cells co-transfected with siRNA-BIRC7 and HIF-1α overexpression vector).
2.2 Western blot assay. Cells were centrifuged (12,000 r/min) at 4°C for 20 min. SDS-PAGE was conducted and proteins were transferred to the PVDF membranes. Then the membranes were blocked and proteins were detected by anti-HIF-1α, anti-Tubulin, anti-lamin B1, anti-GAPDH, anti-BIRC7 (Abcam, USA) before probing with secondary antibodies for 1 h at room temperature. ECL reagent was used to visualize protein bands. Each experiment included triplicate measurements.
2.3 Immuno uorescence staining (IF) assay. The cells were xed with 4% PFA for 15 min. Then 0.05% Triton X-100 was added to permeabilize cells. Cells were blocked with 5% BSA, incubated with primary antibodies overnight at 4°C and FITC-conjugated secondary antibodies at room temperature 1 h, and stained with DAPI for 6 min. Finally, cells were mounted with Aqua Mount. Each experiment included triplicate measurements.
2.4 RNA preparation and RNA sequencing. Total RNA was isolated using Trizol (CWBIO, Beijing, China). RNA quality was assessed with Agilent 2100 Bioanalyzer (Agilent, USA). cDNA library was constructed with Oligo (dT) magnetic beads (Invitrogen, Carlsbad, CA, USA) and random hexamers. The Agilent 2100 Bioanalyzer was used to detect the range of library inserts and the ABI StepOnePlus Real-Time PCR System was used to quantify the library concentration. After passing the quality inspection, sequencing was carried on with the Illumina HiSeq sequencer.
2.5 Data analysis. The original image data le obtained by Illumina HiSeq was converted into the original data by CASAVA base calling analysis, and the result was stored in the FASTQ le format.
2.6 Reads Trimming and Reference Genome Alignment. After the raw reads discarding low-quality readings through perlscript to obtain clean reads, we used HISAT (Kim et al., 2015) to compare the clean reads to the reference genome. For the repeated experiments of biology-free biology, we used the PossionDis algorithm for differential gene detection, and screened genes with |log2(FoldChange)|> 1&qvalue < 0.001 as differentially expressed genes. 2.7 GO and KEGG pathway analysis. Gene Ontology (GO) consists of three main branches, namely, Biological Process, Molecular Function and Cellular Component. Genes or proteins can nd the corresponding GO number through ID correspondence or sequence alignment. The GO number can be used to correspond to Term, that is, functional category or cell location. If the species under study has a relevant GO annotation database, the database wass directly identify exons and introns of genes, locate their boundaries, and the 59 and 39 ends of genes, thereby enabling people to fully understand the complexity of eukaryotic transcriptomes. In addition, RNA-Seq can identify transcription initiation sites (TSS) and new splice variants, and accurately quantify the expression of exons and splice isoforms (Pan et al. ).
In this study, we explored gene expression differences in the response of MM cells to normal and hypoxia by Transcriptome Sequencing. We found that signi cant differences in DEG between normal and Hypoxic-treated MM cells.
analyzed by GO, if not, Blast2GO was used to obtain the GO entry corresponding to each gene. The degree of KEGG enrichment was measured using Rich factor, Qvalue obtained by Fisher's exact test, and the number of genes enriched in this pathway.
Each experiment included triplicate measurements.
2.10 CCK8 assay. The proliferation of M14 cells was determined using CCK8 (Solarbio Science & Technology, Beijing, China). After transfection with si-BIRC7, HIF-1α, and si-BIRC7 + HIF-1α, cells were planted into a 96-well plate at a density of 2,000 cells per well. M14 cells was incubation for 0, 24, 48 and 72 h, then CCK8 solution (10 µl/well) was added in cells for 2 h before the detection. The OD value at 450 nm was recorded, then growth curve was drawn.
2.11 Cell invasion and migration assays. Transwell assy was used to evaluate the invasion and migration of M14 cells using a 24-well transwell insertion system (BD, USA). For invasion assay, matrigel was added to the upper chamber. For migration assay, there was no Matrigel in the transwell chamber. After transfection, cells were planted into the upper chamber, and 500 µl complete medium was added to the lower chamber. After incubation, cells were stained and counted.
2.12 Flow cytometry. Cells from different treatment groups were digested with trypsin. Annexin V binding buffer (Thermo Fisher Scienti c, Shanghai, China) was added to resuspended cells. Then cells were added with Annexin V/FITC mix and PI dye. After incubation in dark for 15 min, cell apoptosis was analyzed. Each experiment included triplicate measurements.
2.13 Statistical analysis. All data were performed using SPSS 18.0 software. The differences of the data in different group was performed using the Student's t test or a one-way ANOVA test. P < 0.05 was considered statistically signi cant.

Results
3.1 Construction of cell hypoxia model. To clear the effect of hypoxia on M14 cells, we create an oxygen-de cient M14 cells. The results in Fig. 1a showed that the expression of HIF1α was induced by hypoxia in 6 h, 12 h, 24 h and 48 h, while the expression of HIF1α in 12 h was highest. IF assay showed the same result that the expression of HIF1α in hypoxia group was signi cantly increased and mainly distributed in the nucleus compared with NC group (Fig. 1b). Subsequently, the results shown in Fig. 1c, suggested that in the cells treated with hypoxia, HIF1α was highly expressed in nucleus compared with cytoplasm. The results indicated that cell hypoxia model was successfully constructed.
3.2 Analysis of RNA-Seq data. To explore transcriptional effects caused by hypoxia, we performed RNA-seq of M14 cells in two conditions: untreated and exposed to hypoxic conditions. More than 95% of reads could be mapped to the genome, and Q30 > 93%. More than nine-tenths of reads could be mapped to the reference genome. In addition, more than four-fths of mapped reads had unique genomic locations (Table 1). The rst 6 up-regulated genes in the hypoxia group were listed in Table 2 compared with NC group. OSCAR, BIRC7, HBA1, SFN, GOLT1A, and BEX2 were signi cantly up-regulated in the hypoxic group (p = 0.000 < 0.01). By comparing the transcriptome data of the two groups, a total of 2601 DEGs were identi ed. Compared with the control group, there were 1517 genes with signi cantly higher expression and 1084 genes with lower expression in the hypoxic group (Fig. 2a). We used the FPKM value of the different genes under different experimental conditions to conduct hierarchical clustering analysis (Fig. 2b). The regions with different colors represent different clustering information. The gene expression patterns in the same group were similar and may have similar function or participate in the same biological process. These results suggested that hypoxia caused extensive transcriptome response in MM cells. GO annotation indicated that DEGs were mainly involved in single biological processes, cellular processes and biological regulation of biological processes, cell and cell part of cellular components, binding and catalytic activity of molecular functions (Fig. 3a). KEGG pathway analysis indicated that DEGs were mainly involved p53 signaling pathway and the PI3K/Akt signaling pathway (Fig. 3b).
3.3 BIRC7 expression is most signi cantly up-regulated. To verify the gene expression data by RNA-seq, we performed qRT-PCR to detect mRNA levels of six genes. As shown in Fig. 4a, six genes in hypoxia group were up-regulated compared with NC group, and BIRC7 expression was most signi cantly up-regulated. Then western blot results showed that BIRC7 was overexpressed in hypoxia group and low expression in hypoxia + HIF-1α inhibitor group (Fig. 4b-4c). These results suggested that BIRC7 may be a downstream effector of HIF-1α. Based on the TCGA database, we analyzed the box plot of BIRC7 on GEPIA database (Loftus et al., 2017). As shown in Fig. 4d, the expression of BIRC7 in SKCM was signi cantly higher than that in normal samples.
3.4 BIRC7 is a downstream effector of HIF-1α. To clear whether BIRC7 was a downstream effector of HIF-1α, we predicted the binding sites on BIRC7 promoter with HIF-1α using StarBase software and dual-luciferase assay. As shown in Fig. 5a, there was a binding site between BIRC7 promoter and HIF-1α. The relative luciferase activity enhanced signi cantly after the transfection of BIRC7 promoter (wt) with hypoxia compared with the transfection of BIRC7 promoter (mut) with hypoxia (Fig. 5b). Then as shown in Fig. 5c, the expression of BIRC7 was decreased in A875 and M14 cells transfected with siBIRC7. CCK8 and colony formation assays results showed that siBIRC7 reversed the effect of hypoxia on promoting the growth of A875 and M14 cells (P < 0.01, Fig. 5d and 5e).
We further analyzed the migration of A875 and M14 cells. The results showed in Fig. 6a that the migrated cell was increased in A875 and M14 cells treated with hypoxia, and decreased after the knockdown of BIRC7 compared with NC group (P < 0.01). siBIRC7 reversed the effect of hypoxia on promoting the migration of A875 and M14 cells. Then as shown in Fig. 6b and 6c, the percentage of apoptotic cells was increased in siBIRC7 group and decreased in the cells treated with hypoxia compared with NC group (P < 0.01). siBIRC7 reversed the effect of hypoxia on inhibiting the apoptosis of A875 and M14 cells.

Discussion
In this study, gene expression differences in the response of MM cells to normal and hypoxia were analyzed. Our data showed for the rst time that transcriptome gene expression pro ling of M14 cell in normal and hypoxia through RNAseq technology. Total Mapped Reads account for more than 95% of Total Reads, indicating high utilization of transcriptome data, so the selected reference genome assembly meets the requirements for information analysis.
Our study also identi ed that a total of 2601 genes were shown to be signi cantly differentially expressed between the two groups. GO terms and KEGG pathway analysis indicated that DEGs were mainly involved in biological regulation and cellular process of biological processes, cell, cell part and organelle of cell components, and the binding and catalytic activity of molecular function. Our ndings were basically consistent with transcriptome analysis of other studies. According to our data, six up-regulated genes were screened from DEG, and qRT-RCR and western blot assays were used to determine BIRC7 as an up-regulated gene. We con rmed that the binding sequence of BIRC7 and HIF-1 was ACGTGCAGT through dual luciferase assay, suggesting that the transcription factor HIF-1 regulated the expression of BIRC7.  hypoxia, HIF1α was highly expressed in nucleus compared with cytoplasm. Immunostaining for HIF1α (green). The nucleus was stained by DAPI (blue). **represents P<0.01 compared with NC.     BIRC7 reverses the effect of hypoxia on tumor migration and apoptosis. (A) The number of migrated cell was increased in A875 and M14 cells treated with hypoxia, and decreased after the knockdown of BIRC7 compared with NC group. (B and C)The percentage of apoptotic cells was increased in siBIRC7 group and decreased in the cells treated with hypoxia compared with NC group. *represents P<0.05, **represents P<0.01 compared with NC.