Weighted Gene Coexpression Network Analysis Identifies Cysteine-Rich Intestinal Protein 1 (CRIP1) as a Prognostic Gene Associated with Relapse in Patients with Acute Myeloid Leukemia.

BACKGROUND Acute myeloid leukemia (AML) is associated with a high relapse rate and poor prognosis. This study aimed to use weighted gene coexpression network analysis (WGCNA) of gene coexpression networks to identify candidate prognostic biomarker genes in patients with AML and to investigate the expression of these genes in the human U937 cell line in vitro. MATERIAL AND METHODS RNA-seq data were retrieved from the Cancer Genome Atlas (TCGA) and included bone marrow samples and survival data of patients with AML (N=151), patients who did not relapse after treatment (N=119), and patients with relapse (N=40). Differentially expressed genes were identified, WGCNA was used to detect functional modules, and survival analysis was performed. The Cell Counting Kit-8 (CCK-8) assay investigated the proliferation of U937 cells transfected with short hairpin RNAs (shRNAs), shCRIP1, shHIST1H1C, and shHIST1H1E. RNA-seq analysis identified gene expression following CRIP1 knockdown. RESULTS Eighty-two genes were associated with both relapse and prognosis in patients with AML. There were two prognosis-related gene modules in the coexpression network. In the coexpression network, the histone cluster 1 H1 family member gene, HIST1H1C had the maximum relapse fold change, HIST1H1E had the lowest survival p-value, and the cysteine-rich intestinal protein 1 (CRIP1) gene had the most edge numbers and was significantly associated with poor prognosis (P=0.0165786). RNA-seq data showed that there was a significant difference in gene expression after CRIP1 knockdown in U937 cells. CONCLUSIONS WGCNA of gene coexpression networks identified CRIP1 as a potential prognostic biomarker gene in patients with AML.

Background

Acute myeloid leukemia (AML) is characterized by the abnormal proliferation and differentiation of the immature myeloid cells [1,2]. The pathogenesis of AML is a complex process involving several molecular events and signaling pathways [2]. Currently, chemotherapy and allogeneic bone marrow transplantation are the main treatments for AML [3]. Although the clinical management of AML has made considerable progress in the past few decades, many patients experience poor prognosis and relapse on follow-up [4]. Although early diagnosis of AML can improve patient prognosis, relapse following a period of remission remains a major clinical challenge [5,6]. Therefore, accurate assessment of prognosis of patients with AML is important for clinical management, and indicators of prognosis may guide clinical decisions for the treatment of patients with AML [7]. The prediction of accurate clinical outcome in patients with AML depends on accurate predictive and prognostic biomarkers for patients with an increased risk of relapse and poor prognosis.

Aberrant gene expression in human malignancy is associated with oncogenesis, cell proliferation, and resistance to chemotherapy [8]. Clinically, aberrantly expressed genes associated with malignancy have a potential role as diagnostic or prognostic biomarkers. Abnormal gene expression frequently occurs in AML, and have roles in the pathogenesis of AML [4]. Identification of molecular events associated patient prognosis may also provide a deeper understanding of the pathogenesis of AML, and several abnormally expressed genes have been identified in patients. For example, in 2018, Zhang et al. showed that overexpression of the H19 gene was a prognostic and predictive biomarker associated with reduced overall survival (OS) and a lower rate of complete remission (CR) in patients with AML [9].

Therefore, this study aimed to use weighted gene coexpression network analysis (WGCNA) of gene coexpression networks to identify candidate prognostic biomarker genes in patients with AML and to investigate the expression of these genes in the human U937 cell line in vitro.

Material and Methods

Data acquisition

RNA-seq data relating to patients with acute myeloid leukemia (AML) were obtained from The Cancer Genome Atlas (TCGA) database () and included bone marrow samples and survival data of patients with AML (N=151), patients who did not relapse after treatment (N=119), and patients with relapse (N=40). The RNA-seq data included only patients who had complete clinical information. This study was fully compliant with the publication guidelines provided by TCGA. Data were acquired from TCGA, which is a public database, and ethics committee approval for data analysis was not required.

Identification of differentially expressed genes associated with patient relapse in AML

The RNA-seq data of 119 bone marrow samples from patients with primary AML and 40 bone marrow samples from patients with relapse of AML using DESeq2 package were screened for differentially expressed genes associated with relapse of AML [10]. Differentially expressed genes were identified with a threshold of P-value of <0.01 after adjustment, and |log2FC| >1.

Identification of genes associated with patient prognosis in AML from functional enrichment analysis

Univariate Cox proportional hazards regression analysis was performed using survival R package to explore the genes associated with patient prognosis in AML. Genes with a P-value <0.05 were identified as prognostic genes. The Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways enriched by these prognostic genes using were investigated using the clusterProfiler package [11]. A P-value <0.05 was considered to be statistically significant.

Genes associated with both recurrence and prognosis in patients with AML

The online Venn diagram tool () was used to obtain genes associated with both recurrence and prognosis obtained from the intersection of the two groups of genes.

Coexpression network analysis and overall survival (OS) analysis of patients with AML

The coexpression network of genes associated with both recurrence and prognosis was constructed using the Weighted Gene Coexpression Network Analysis (WGCNA) package, as previously described [12]. The genes that were extracted from each gene coexpression module were used to perform the survival analysis. Overall survival analysis was conducted using the survival package. A P-value <0.05 was considered statistically significant.

Cell culture and transfection

U937 cells were cultured in RPMI 1640 medium with 10% fetal bovine serum (FBS) and grown at 37°C in an incubator containing 5% CO2. Scrambled shRNA (shControl), CRIP1 shRNA, HIST1H1C shRNA, HIST1H1E shRNA, and pcDNA3.1-Control (Vector), pcDNA3.1-CRIP1, pcDNA3.1-HIST1H1C and pcDNA3.1-HIST1H1E were used (Invitrogen Life Technologies, Carlsbad, CA, USA).

Quantitative reverse transcription-polymerase chain reaction (RT-qPCR)

Total RNA was extracted from U937 cells with TRIzol reagent (Invitrogen Life Technologies, Carlsbad, CA, USA). To examine the expression of HIST1H1C, HIST1H1E, and CRIP1, total RNA was transcribed into cDNA by PrimeScript RT Master Mix (Takara, Dalian, China), according to the manufacturer’s instructions. The RevertAid First Strand cDNA synthesis Kit (Thermofisher Scientific, Waltham, MA, USA) was used for reverse transcription into cDNA using SYBR Green Supermix (Bio-Rad Laboratories, Hercules, CA, USA). The experimental results were automatically calculated using BIO-RAD CFX Manager version 3.1 RT-qPCR analysis software (Bio-Rad Laboratories, Hercules, CA, USA). All experiments were independently performed in triplicate. The primer sequences of the target genes are shown in Table 1. Glyceraldehyde-3-phosphate dehydrogenase (GAPDH) served as an internal control. The relative mRNA expression levels of the target genes were calculated by the 2−ΔΔCt method.

Table 1

The primer sequences used for polymerase chain reaction (PCR).

Target gene	Primer sequence
CRIP1	5′-AAGTGCGACAAGGAGGTGTAT-3′5′-AAGTGCGACAAGGAGGTGTAT-3′
HIST1H1C	5′-CCGCCTCTAAAGAGCGTAGC-3′5′-AGACCAAGTTTGATACGGCTG-3′
HIST1H1E	5′-CTGGGGAAGCCAAGCCTAAG-3′5′-CTGCTTTCGCCTTTTTCGGG-3′
GAPDH	5′-GGAGCGAGATCCCTCCAAAAT-3′5′-GGCTGTTGTCATACTTCTCATGG-3′

Western blot

Total protein was extracted from U937 cells using RIPA lysis buffer (Sigma-Aldrich, St. Louis MO, USA). The extracted protein was separated by sodium dodecyl sulfate-polyacrylamide gel electrophoresis (SDS-PAGE) and transferred to a polyvinylidene fluoride (PVDF) membrane (Merck Millipore, Burlington, MA, USA), followed by blocking with a 5% solution of dried skimmed milk powder for 1 h. The membrane was incubated with the primary antibodies overnight at 4°C, and incubation for 1 h at room temperature with the secondary antibody. Proteins were visualized using ECL chemiluminescence and analyzed using the ChemiDoc XRS+ electrophoresis imaging system (Bio-Rad Laboratories, Berkeley, CA, USA). Primary antibodies included anti-CRIP1 (Abcam, Cambridge, MA, USA), anti-HIST1H1C (Abcam, Cambridge, MA, USA), and anti-HIST1H1E (Abcam, Cambridge, MA, USA), and anti-GADPH (Abcam, Cambridge, MA, USA). GADPH served as an internal control.

RNA-seq analysis

RNA-seq data analysis was performed to investigate the differential expression of genes after CRIP1 knockdown in the U937 cell line compared with the control group. Total RNA was extracted from U937 cells using TRIzol reagent (Life Technologies, Carlsbad, CA). Electrophoresis was used to assess RNA integrity and DNA contamination on denaturing agarose gels. The RNA-seq library was constructed. SOAPnuke version 1.5.2 was used for data filtering including removal of inclusions, as well as unknown base N content >5% and low quality reads [13]. After obtaining clean reads, Hierarchical Indexing for Spliced Alignment of Transcripts (HISAT) was used to compare clean reads with the reference genome sequence [14]. Bowtie2 was used to compare clean reads with the reference gene sequence for alignment [15,16]. Whole transcriptome sequencing was performed, including the randomness, coverage, and distribution sequencing saturation of the transcript. The length of the reference transcriptome sequence was counted using a Perl script. We predicted genes with the ability to encode transcription factors (TF). The transcription factor family of the gene was classified and using AnimalTFDB2.0. Bowtie2 was used to compare clean reads with gene sequences. Then, RSEM, a software package for RNA-seq reads computational genes and expression levels of transcript isoforms was used to calculate gene expression levels for each sample.

Pearson’s correlation coefficient of all gene expression levels between each sample was calculated and reflected in the form of a heat map using the cor () function in the R software to reflect the correlation of gene expression between samples. The correlation coefficient reflected the similarity of the overall gene expression between samples, as the higher the correlation coefficient, the more similar the gene expression level. A box plot, a density map, and a stacked histogram were drawn based on the expression of the fragments per kilobase of exon model per million reads mapped (FPKM) information of each sample. Differentially expressed genes were screened based on a false discovery rate (FDR) ≤0.001 and a fold change (FC) >2. A Venn diagram was used to show the gene between the two groups. The FPKM values of each comparative group of differential genes were clustered and analyzed by the pheatmap package in the R software. GO and KEGG enrichment analysis were performed using the phyper function in the R software.

Cell counting kit-8 (CCK-8) assay

U937 cells with or without knockdown of HIST1H1C, HIST1H1E, and CRIP1 were seeded in a 96-well plate. The cell viability was measured using the CCK-8 assay (Dojindo, Tokyo, Japan) after cell culture for 0, 1, 2, 3, and 4 days.

Validation of CRIP1, HIST1H1C, and HIST1H1E using the Gene Expression Profiling Interactive Analysis (GEPIA) database

Gene Expression Profiling Interactive Analysis was used (GEPIA; ) using data from TCGA and Genotype-Tissue Expression (GTEx) datasets to analyze CRIP1, HIST1H1C, and HIST1H1E expression differences between AML and normal samples. Overall survival (OS) analysis was performed. A P-value <0.05 was considered statistically significant.

Results

Differentially expressed genes associated with relapse in patients with acute myeloid leukemia (AML)

Based on RNA-seq data from The Cancer Genome Atlas (TCGA), we compared the differentially expressed genes between the primary bone marrow-derived AML samples (119 cases) and the relapsed bone marrow-derived samples (40 cases). Finally, a total of 24,991 genes were obtained, and 1,014 differentially expressed genes were screened under the conditions of a false discovery rate (FDR) <0.01 and |log2FC| >1 using the R package, DESeq2.

Identification of genes associated with prognosis in AML and potential biological pathways

Univariate Cox proportional hazards regression analysis showed that 2,043 genes were associated with prognosis of AML (P<0.05). Using the R package, clusterProfiler, the Kyoto Encyclopedia of Genes and Genomes (KEGG) database pathway enrichment analysis was performed, and the results showed that these prognostic genes were mainly enriched in signaling pathways for alcoholism, systemic lupus erythematosus, and viral carcinogenesis (Table 2).

Table 2

The main Kyoto Encyclopedia of Genes and Genomes (KEGG) database enrichment pathways of 2,043 genes associated with prognosis in patients with acure myeloid leukemia (AML).

ID	Description	P-value	Count
hsa05034	Alcoholism	1.20E-08	42
hsa05322	Systemic lupus erythematosus	3.23E-07	32
hsa05203	Viral carcinogenesis	3.54E-05	37

Identification of genes associated with both relapse and prognosis in patients with AML

By overlapping differentially expressed genes associated with relapse and genes associated with prognosis, 82 genes were identified as being associated with both relapse and prognosis in AML, as shown in Figure 1. These 82 genes are listed in Table 3. Most genes were highly expressed in patients with relapse of AML and were closely associated with prognosis.

Figure 1

Venn diagram showing 82 genes associated with both recurrence and prognosis in patients with acute myeloid leukemia (AML).

Table 3

The 82 genes associated with both relapse and prognosis in patients with acute myeloid leukemia (AML).

Total

Genes

82

KCNG1 LOC101929243 NEIL1 HIST1H1E HIST1H2AD ASCL5 NANOS3 HIST1H2BJ METTL12 NUAK1 CD151 ZFPM2 HIST1H3E IGLL5 CHCHD5 PPDPF SOX15 CAHM LOC100996720 HIST2H2AC KCNH3 HIST2H3A HIST1H2AM LINC01635 HIST1H2BI C2CD4D CLEC17A CRIP1 SH3BGRL3 HIST1H3D HS3ST2 CSTB RPP25L STX1B HIST1H2BD TGFB1I1 HIST1H1C HIST3H2BB HIST1H2BO BCYRN1 HIST2H2BE RENBP FLJ20021 GPR162 FAM46B DBNDD2 HAGHL HIST1H2AE METRN THAP8 HIST1H2AL HIST2H2AA4 ZGLP1 NUDT14 HDHD5-AS1 FAM166B GABRA2 GREM1 NTRK3 GPR88 HSD3B7 PACRG LINC00309 GADD45G CD5L SLC35E4 TSPO C4orf47 LOC283440 TLE6 PPP2R2B HSPA1B MAPK15 HIST2H2BF BOLA1 LINC00886 ACOX2 PTGIS KRT23 HIST2H3C JOSD2 HIST1H2BF

Construction of the coexpression network

There were 44 nodes in the coexpression network. These genes were associated with both relapse and prognosis in AML and were mainly divided into two modules, including the histone family module and the transcription factor module (Figure 2). In the coexpression network, the histone cluster 1 H1 family member gene, HIST1H1C had the maximum relapse fold change (FC) and the minimum relapse FDR. HIST1H1E had the lowest survival P-value, and CRIP1 had the most edge number in the network. Survival analysis results of three genes are shown in Figure 3. High expression of the three genes was associated with poor prognosis in patients with AML.

Figure 2

The coexpression network based on 82 genes associated with both recurrence and prognosis of acute myeloid leukemia (AML). There are 44 nodes in the coexpression network and two modules, including the histone family module and the transcription factor module.

Figure 3

Results of the overall survival (OS) analysis associated with the expression of CRIP1, HIST1H1E, and HIST1H1C. (A) CRIP1 (P=0.0165786). (B) HIST1H1E (P=0.0015068). (C) HIST1H1C (P=0.0111662).

CRIP1, HIST1H1C, and HIST1H1E and proliferation of U937 cells

To determine whether CRIP1, HIST1H1C, and HIST1H1E could affect cell proliferation, the Cell Counting Kit-8 (CCK-8) assay was performed in U937 cells. Firstly, U937 cells were transfected with shCRIP1, shHIST1H1C, and shHIST1H1E, respectively (Figure 4A–4D). Quantitative reverse transcription-polymerase chain reaction (RT-qPCR) and Western blot were performed to determine the effects of transfection. As shown in Figure 4E–4G, the expression of CRIP1, HIST1H1C, and HIST1H1E was significantly inhibited after transfection at the mRNA expression level. Western blot showed that the protein expression levels of CRIP1, HIST1H1C, and HIST1H1E were also significantly reduced in U937 cells transfected with shCRIP1, shHIST1H1C, and shHIST1H1E (Figure 4H–4J). The CCK-8 assay results showed that U937 cell proliferation was inhibited after transfected by shCRIP1, while cell proliferation was promoted after transfected by shHIST1H1C or shHIST1H1E (Figure 5). The findings showed that expression of CRIP1, HIST1H1C, and HIST1H1E could affect cell proliferation of U937 cells in vitro.

Figure 4

The effects of transfection of U937 cells transfected with short hairpin RNAs (shRNAs), shCRIP1, shHIST1H1C, and shHIST1H1E. (A–D) Representative images of U937 cells transfected with shNC, shCRIP1, shHIST1H1C, and shHIST1H1E, respectively. (E–G) Quantitative reverse transcription-polymerase chain reaction (RT-qPCR) results show the mRNA expression levels of CRIP1, HIST1H1C, and HIST1H1E in U937 cells transfected with shCRIP1, shHIST1H1C, and shHIST1H1E. (H–J) Western blot results show the protein expression levels of CRIP1, HIST1H1C, and HIST1H1E in U937 cells transfected with shCRIP1, shHIST1H1C, and shHIST1H1E. ** P<0.01; *** P<0.001.

Figure 5

Cell counting kit-8 (CCK-8) assay shows the proliferation of U937 cells transfected with short hairpin RNAs (shRNAs), shCRIP1, shHIST1H1C, and shHIST1H1E.

Gene expression following CRIP1 knockdown in U937 cells using RNA-seq

RNA-seq analysis was performed in U937 cells after CRIP1 knockdown to explore the potential function of CRIP1 further. The box plot shows the distribution of gene expression levels in two samples (Figure 6A), and the degree of dispersion of the distribution of data are shown. The density map depicts the trend of the gene levels in the sample as the expression level changes and reflects the concentration range of gene expression in the sample (Figure 6B). To show the number of genes in each fragments per kilobase of exon model per million reads mapped (FPKM) interval of each sample, statistical analysis was performed on the number of genes in three cases, FPKM ≤1, FPKM 1–10, and FPKM ≥10 (Figure 6C). The volcano plot and hierarchical cluster analysis results show the differentially expressed genes in Figure 7A and 7B. As shown in Figure 7C, these differentially expressed genes were mainly enriched in cellular processes and cell binding. The KEGG enrichment analysis results showed that these genes were mainly enriched in signal transduction, infectious diseases, and viral cancers (Figure 7D).

Figure 6

Gene expression in U937 cells transfected by short hairpin CRIP1 (shCRIP1) using RNA-seq. (A) The box plot. (B) The density map. (C) The stacked histogram.

Figure 7

Differentially expressed genes in U937 cells transfected with short hairpin CRIP1 (shCRIP1) using RNA-seq. (A) Volcano plot of differentially expressed genes. Red indicates upregulated genes, and blue indicates down-regulated genes. (B) Hierarchical cluster analysis of differentially expressed genes. (C) The top 20 Gene Ontology (GO) results enriched by differentially expressed genes. (D) The top 20 Kyoto Encyclopedia of Genes and Genomes (KEGG) database enrichment pathways enriched by differentially expressed genes.

Using the GEPIA database, CRIP1 was found to be upregulated in AML (Figure 8A). Also, HIST1H1C and HIST1H1E were both down-regulated in AML (Figure 8B, 8C). High CRIP1 expression was associated with poor survival time (Figure 8D), and high expression levels of HIST1H1C and HIST1H1E were associated with poor prognosis (Figure 8E, 8F). Therefore, CRIP1, HIST1H1C, and HIST1H1E were genes associated with reduced prognosis in patients with AML in this study.

Figure 8

Validation of CRIP1, HIST1H1E, and HIST1H1C in acute myeloid leukemia (AML) using the Gene Expression Profiling Interactive Analysis (GEPIA) database. (A–C) The expression of CRIP1, HIST1H1E, and HIST1H1C in AML compared with normal samples. * P<0.05. (D–F) Overall survival (OS) analysis of CRIP1, HIST1H1E, and HIST1H1C in AML.

Discussion

Acute myeloid leukemia (AML) is a biologically and clinically heterogeneous disease [17]. Although there have been recent advances in supportive care and prognostic risk stratification, the prognosis for patients with AML remains poor [18,19]. Therefore, there is a need to identify novel prognostic markers and therapeutic targets for AML. In the present study, 24,991 genes associated with relapse in patients with AML were obtained from The Cancer Genome Atlas (TCGA) database by comparing 119 bone marrow samples from patients at primary diagnosis with AML, and 40 bone marrow samples from patients with relapsed. Under the condition of a false discovery rate (FDR) <0.01 and |log2FC| >1, there were 1,014 differentially expressed genes identified that were screened. Also, 2,043 genes were identified as prognostic markers, which had P<0.05 in univariate Cox proportional hazards regression analysis. These 2,043 genes were associated with prognosis and were mainly enriched in several signaling pathways that included alcoholism, systemic lupus erythematosus (SLE), and viral carcinogenesis. It is of interest to note that disulfiram (DSF), which is used to treat alcohol dependency, also has an anti-cancer effect. A recently reported in vitro study using human AML cell lines, showed that DSF combined with copper increased the expression of reactive oxygen species (ROS), resulting in cell cycle arrest and AML cell apoptosis [20]. Also, the incidence of AML is increased in patients with autoimmune diseases, such as SLE, and is also increased after exposure to cytotoxic agents [21,22].

In this study, 82 genes were identified that were associated with both recurrence and prognosis in AML (Table 3). Using weighted gene coexpression network analysis (WGCNA) of gene coexpression networks, we identified two prognosis-related coexpression modules, which included the histone family module and the transcription factor module. Based on gene expression profiling, WGCNA has been widely used to identify functional coexpression gene modules and has provided insight into the biological signaling networks associated with phenotypic traits of interest [23,24]. In patients with relapse of AML, the high expression of the genes associated with poor prognosis may result in disturbances in transcriptional regulation, DNA repair, DNA replication, and chromosomal stability, as well as abnormal karyotype due to abnormal expression of the protein [25-28]. The karyotypes of most patients with AML are abnormal, and abnormal key genes may lead to the incorrect assembly of chromosomes in cells of patients with AML, errors in the conformation of DNA molecules, or transcription errors in genetic information [29]. In the coexpression network described in the present study, genes were highly expressed in bone marrow samples from patients who underwent relapse of AML, and HIST1H1C had the largest recurrence fold change (FC) and the smallest recurrence FDR, while HIST1H1E had the smallest survival P-value, and CRIP1 had the largest number of edges in the network. High expression of these three genes was associated with poor prognosis. Also, the findings from the Cell Counting Kit-8 (CCK-8) assay showed that these three genes could affect the proliferation of U937 cells.

Previously published studies have shown that CRIP1 may affect prognosis and cellular processes in human disease. In 2013, Ludyga et al. showed that down-regulation of CRIP1 was significantly associated with worse prognosis in patients with breast cancer [30]. In 2014, Xie et al. proposed a combined method using CRIP1 protein and sequence data for tissue samples from patients with breast cancer as a sensitive detection method that may be used clinically [31]. CRIP1 as a prognostic biomarker has also been investigated in prostate cancer, intestinal-type gastric carcinoma, and endometrial carcinoma [32-34]. However, the prognostic role of CRIP1 in AML has not been previously reported.

In the present study, RNA-seq was performed to investigate further the differences in gene expression following knockdown of CRIP1 in the U937 cell line, which showed that there was a significant difference in gene expression compared with the normal cell line. KEGG and GO enrichment analysis showed that these differentially expressed genes were enriched in several signaling pathways, including cell growth, signal transduction, replication and repair, and immune processes that may be associated with the development and progression of AML. Differentially expressed genes that are enriched in immune pathways may be future targets for immunotherapy in AML [35]. In this study, analysis of the Gene Expression Profiling Interactive Analysis (GEPIA) database showed that CRIP1 was upregulated in patients with AML and was associated with a shorter survival time compared with patients with AML and down-regulated CRIP1.

Conclusions

This study aimed to use weighted gene coexpression network analysis (WGCNA) of gene coexpression networks to identify candidate prognostic biomarker genes in patients with AML and to investigate the expression of these genes in the human U937 cell line in vitro. WGCNA of gene coexpression networks identified CRIP1 expression in bone marrow samples to be associated with both relapse and reduced prognosis in AML. These preliminary findings support the need for future clinical studies to evaluate the role of CRIP1 as a prognostic biomarker in patients with AML.