Gene Expression Analysis of Human Papillomavirus-Associated Colorectal Carcinoma

Purpose Human papillomavirus (HPV) antigens had been found in colorectal cancer (CRC) tissue, but little evidence demonstrates the association of HPV with oncogene mutations in CRC. We aim to elucidate the mutated genes that link HPV infection and CRC carcinogenesis. Methods Cancerous and adjacent noncancerous tissues were obtained from CRC patients. HPV antigen was measured by using the immunohistochemical (IHC) technique. The differentially expressed genes (DEGs) in HPV-positive and HPV-negative tumor tissues were measured by using TaqMan Array Plates. The target genes were validated with the qPCR method. Results 15 (31.9%) cases of CRC patients were observed to be HPV positive, in which HPV antigen was expressed in most tumor tissues rather than in adjacent noncancerous tissues. With TaqMan Array Plates analyses, we found that 39 differentially expressed genes (DEGs) were upregulated, while 17 DEGs were downregulated in HPV-positive CRC tissues compared with HPV-negative tissues. Four DEGs (MMP-7, MYC, WNT-5A, and AXIN2) were upregulated in tumor vs. normal tissues, or adenoma vs. normal tissue in TCGA, which was overlapped with our data. In the confirmation test, MMP-7, MYC, WNT-5A, and AXIN2 were upregulated in cancerous tissue compared with adjacent noncancerous tissue. MYC, WNT-5A, and AXIN2 were shown to be upregulated in HPV-positive CRC tissues when compared to HPV-negative tissues. Conclusion HPV-encoding genome may integrate into the tumor genomes that involved in multiple signaling pathways. Further genomic and proteomic investigation is necessary for obtaining a more comprehensive knowledge of signaling pathways associated with the CRC carcinogenesis.


Introduction
Colorectal cancer (CRC) is the primary cause of cancer mortality. Although progress has been made, the long-term survival of patients with metastatic disease is still poor [1]. Both genetic and environmental factors contribute to the pathogenesis of CRC; about 75% of CRC cases are not hereditary and occur spontaneously [2]. The gut microbiota is closely correlated with the progression of CRC [3,4]. The dysbiosis is associated with the genesis and evolution of CRC, but the relationship remains unclear [5][6][7][8]. Accumulating evidence shows that colonizing microbes can drive cancer development and progression by direct or indirect effects on host tissues, potentially through inflammatory pathways or carcinogenic microbial metabolites; a causative link was found between chronic inflammation and CRC [9][10][11]. As microorganisms establish a persistent infection in host cells, either in the form of silent or active infections, barrier deterioration may be a major contributor to colorectal tumorigenesis by microbial products that trigger tumorelicited inflammation [12,13].
HPV carries out its life cycle in either the mucosal or the cutaneous stratified squamous epithelia. An oncogenic HPV infection is associated with several malignancies, including cervical cancer [14,15] and lung cancer [16]. To date, 200 distinct HPV genotypes have been identified, of which at least 18 belong to the "high-risk" group that is chiefly responsible for the development of cancer [17,18]. Recently, studies have demonstrated the possible involvement of HPV in CRC [19], owing to the detection of HPV antigen or DNA in CRC tissues [20][21][22]. However, the findings are not consistent. HPV DNA was positive in 60 (83.3%) of the 72 total cancerous colorectal samples, and no HPV DNA was observed in any of the noncancerous tissues in one study [22], while in another study, HPV DNA also found in 53% of healthy mucosal tissues [23]. The frequency of HPV DNA in tumor tissues was higher than that in nontumor colorectal tissues in a Chinese population [24]. Evidence of HPV existing in CRC suggests its role in CRC carcinogenesis. Persistent with a high-risk type of HPV infection is responsible for the premalignant lesions and progression of CRC, which correlated with the late clinical stage [25][26][27]. Theoretically, the effect of a virus on cancer may be mediated through the integration of the mutagenic viral genome into the host genome, expression of oncogenic viral proteins, or inhibition of tumor-suppressive genes [28][29][30]. HPV-encoding genome may integrate into the tumor genome, which is considered an essential step in the malignant progression [31]. Currently, there is little evidence demonstrating the association of HPV infection with oncogenic mutations in CRC [32], and the type of oncogenes involved with HPV remains to be elucidated.
Gene analysis with microarray technology has shown great potential in discriminating sophisticated gene profiling, simultaneously mapping thousands of genes in a single sample, and giving a measurement of articulated gene expression patterns [33]. To clarify the molecular mechanism underlying HPV-associated gene mutations in CRC, we used Taq-Man Array Plates to detect differentially expressed genes (DEGs) in HPV viral-positive and HPV viral-negative CRC tissues. To overcome the limitation of the small sample size, we also took advantage of the bioinformatics technique; the significant genes associated with CRC screened among the pool of all DEGs in The Cancer Genome Atlas (TCGA). Finally, we validated target genes in CRC tissues by conventional quantitative assays (RT-PCR) and aimed to elucidate whether there is a specific gene mutation pattern involving HPV-related CRC. Identification of genome aberrations in CRCs may help to clarify the mechanism of pathogenesis involved with virus-related CRC.  The CRC patients were selected based on  histologically confirmed CRC at TNM stage I-IV disease with  or without distant metastases. An individual with a history of polyps, adenomas, or other diseases related to cancer was excluded. Tumor types were evaluated according to histology, immunohistochemistry (IHC), and chromatin assays. Clinical and histopathologic staging at diagnosis was determined in all patients by combining histopathologic findings with surgical records and operative imaging. All selected tissues were based on histopathological diagnoses; adjacent noncancerous tissue was collected from an area that was 15 cm distal from the tumor. Fresh frozen tissues were obtained from CRC patients before they underwent curative surgical resection. Tumor and the nontumor regions were distinguished using H&E matched slides and microdissected to pinpoint the tumor as well as nontumor areas. Samples were kept at -80°C until measurement. The demographic, clinical, and histopathological data of patients were recorded. Clinical data including race, gender, associated past medical history, medication use, family history of colonic cancer, and laboratory tests were collected for all patients.

Immunohistochemical Analyses.
Immunohistochemical staining was performed on paraffin sections 4 μm thick. The sections were cut, dewaxed in xylene, and rehydrated in decreasing concentrations of ethanol. HPV was measured in the tumor tissues and tumor-adjacent tissues using an IHC technique according to the manufacturer's instructions. Considering the HPV's pathogenesis, we had chosen the anti-HPV high-risk subtype antibodies, which is a monoclonal antibody, reacting with HPV subtypes 6, 11, 16, 18, 31, 33, 42, 51, 52, 56, and 58. The details of antibody are HPV mouse monoclonal antibody IgG1, ab75574, Lot. GR207368-5, Abcam, U.K. In negative control sections, the primary antibody was omitted. Virus-infected monocytes or lymphocytes from the same patient and cervical squamous cell carcinoma were used as the positive control for staining. The sections were graded based on the estimated percentages of virus protein-positive cells. IHC staining was evaluated and graded by two pathologists in isolation from the investigator who performed the real-time PCR.

CRC-Associated Oncogene Measurement.
Genes expressed in HPV-positive and HPV-negative tumor tissues were measured by using TaqMan Array Plates. Fresh specimens were obtained from the initial diagnosis of CRC patients who underwent curative surgical resection; adjacent noncancerous tissue was taken far mostly from the tumor. Total RNA was extracted from the tissue using a PureLink RNA Mini Kit (Applied Biosystems, USA) according to the manufacturer's instructions. The RNA concentration was measured using the NanoDrop 1000 spectrophotometer (Thermo Fisher Scientific), and RNA quality was evaluated by using the BioAnalyzer 2100 microcapillary electrophoresis system (Agilent Technologies, Inc., Santa Clara, CA, USA). The RNA concentration was adjusted 75 ng/L, so that the 2 BioMed Research International RNA concentration of each sample was consistent, and then, reverse transcription was performed. RNA was converted into cDNA by using a High-Capacity cDNA Reverse Transcription Kit (Applied Biosystems, USA). The cDNA microarray analysis was performed using the TaqMan® Array Human Colorectal Cancer Metastasis 96-well plate that contains 92 assays for CRC metastasis-associated genes and four assays to candidate endogenous control genes. Each plate includes predefined probes and endogenous controls, ready for accurate assessment of an entire gene signature in one simple experiment (Thermo Fisher Scientific). Each patient's sample was measured in the tumor tissues and tumoradjacent tissues simultaneously, and the gene expression level was calculated by gene level in tumor tissue minus gene level in tumor-adjacent tissue. Each cDNA was analyzed using real-time PCR on the "mixed" TaqMan Array Cards, and Ct values for target genes were obtained. All Ct values were normalized with an inner control (GAPDH) gene, and ΔCt values were calculated. Relative gene expression was computed using the relative expression level for each gene. The 2ΔCt values were calculated by SDS Software v2.3 and exported. Relative expression values were mean-centered, and heat maps were generated using the "Cluster" and "Tree-View" software tools.

Bioinformatic Analyses.
To avoid the limitation of a small sample size that may not provide a robust prediction, we utilized the TCGA database. The gene expression profile of GSE20916-RAW.tar was downloaded from Gene Expression Omnibus (GEO, http://www.ncbi.nlm.nih.gov/geo/). The GSE20916 dataset contained 145 samples, including 105 CRC tumor tissues and 40 normal colon epitheliums. The hierarchical clustering analysis was used to categorize the data into two groups that had similar expression patterns in CRC and normal colon epitheliums. We used a classical t -test to identify DEGs with a change ≥ twofold and defined P < 0:05 to be statistically significant. The Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses were performed to analyze the DEGs at the functional level. The DEGs were screened using Expression Console Affymetrix.
2.6. Validation of Target Gene Selection. Target genes were selected after analyzing the data from the TaqMan Array Plates and data mining with bioinformatics. Then, we validated the target genes using tumor tissues of 47 CRC patients. TaqMan quantitative PCR was used to determine the expression of the DEGs in the cancerous and adjacent noncancerous tissues (AXIN2 (ID. Hs01063168_m1), MMP-7 (ID. Hs01042796_m1), MYC (ID. Hs99999003_m1), WNT5A (ID. Hs00998537_m1), and GAPDH (ID. Hs99999905_ m1), Thermo Fisher Scientific, USA). Specifically, amplified RNA from each sample was reverse-transcribed using Super-Script II reverse transcriptase and random primers (Invitrogen). cDNA from each reverse-transcribed sample was PCR amplified using TaqMan quantitative PCR. The quantification of the gene copy number, together with the human albumin gene, was performed with the 7500 Real-Time PCR System (Applied Biosystems, Courtaboeuf, France).
The expression values for each gene were normalized with GAPDH.

Statistical Analyses.
Statistical analyses were performed using SPSS (Version 10.0). Pearson's chi-squared analysis and Fisher's exact test were employed to compare the differences in categorical variables between patient groups. The relative gene expression levels were compared between various subgroups using the Kruskal-Wallis test (42 groups) and Mann-Whitney U test (two groups). DEGs were determined by Welch's t-test, and precise targets were selected from the DEGs using the prediction analysis for microarrays (PAM) software package. Significance was established when the statistical tests returned P values < 0.05. Software R (v3.1.1) was used to determine the significance in the Taq-Man Array Plates, comparing genes expressed in CRC-HPV+ vs. CRC HPV tumor tissues (heat map).

The Clinicopathological Characteristics of CRC Patients.
To understand whether CRC patients associated with HPV would have unique features, we investigated their clinicopathological characteristics in this study; the parameters included age, gender, tumor location, histological grading, lymphatic metastasis, and clinical staging. The 47 patients included 4 cases of mucinous adenocarcinoma and 43 cases of nonmucinous adenocarcinoma. Males comprised a more significant proportion of CRC patients compared to females. The majority of tumor regions were located in the sigmoid colon and rectum (Table 1).

HPV Antigens in CRC Tissue and Adjacent Noncancerous
Tissue. IHC staining showed that 15 (31.9%) cases were HPV positive, while 68.1% were HPV negative. HPV antigen was expressed in most cancerous tissues but not in adjacent noncancerous tissues in some CRC patients (Figure 1). At TNM stage I, all adjacent noncancerous tissues were viral negative, but at stage IV, half of the adjacent noncancerous tissues were viral positive. At stages II-III, HPV expression was significantly higher in cancerous tissue compared to adjacent noncancerous tissue. Our results indicated that cancerous tissues had higher levels of HPV expression than adjacent noncancerous tissue.

Bioinformatics Analyzed the DEGs between CRC Tumor
Tissue and Nontumor Tissue. Additionally, we downloaded 145 raw CEL files from the NCBI Gene Expression Omnibus (GEO) database and normalized them with GeneChip RMA (GCRMA). Gene expression data were normalized for each tissue type by computing the Robust Multichip Average (RMA) directly from the Affymetrix. CEL were files obtained for tumors, adenoma, and healthy samples. A total of 3500 genes were identified, and the oncogenes associated with CRC were in the range. From these data, we found 16 oncogenes with a significantly different expression between the tumor and healthy tissues. Compared with healthy tissue, nine genes were upregulated between onefold and twofold in the CRC group: MMP1, MMP7, MMP9, WNT5A, IL-6, PTGS2, MYC, BIRC5, and CCND1. Seven genes upregulated between 0.5-fold and onefold: AXIN2, VEGFA, TLR2, STAT1, LEF1, MEH2, and TLR4. The downregulated genes (0.5-fold to 1.0-fold) include IL-6R, TLR3, PIGER2, MAP2K6, ARRB1, PTGER4, CDKN1A, MAPK3, and TLR7. The DEG expression heat map is shown in Figure 3. When adenoma is compared with the healthy tissue, we found that six genes increased between onefold and twofold: MMP1, MMP7, WNT5A, MYC, AXIN2, and TLR4, while six genes increased between 0.5-fold and onefold: MSH2, TP53, NOS2, CCND2, BAX, and PTGER2. The downregulated genes include IL-6R, TLR3, ARRB1, and TLR7.
3.5. Target Gene Validation. To expand the sample size, the TCGA database was exploited to search the DEGs with significantly different expression in CRC and non-CRC tissues; combining with the DEGs from our microarray analysis, we found that MMP7, WNT5A, MYC, and AXIN2 genes upregulated in tumor vs. healthy tissues, adenoma vs. healthy tissue in TCGA; these genes were overlapped between the TaqMan Array Plates data and TCGA data. To validate the results, we detected them in the 47 samples of CRC using the qRT-PCR assay. For all the four genes, the results were concordant with the TaqMan Array Plates data. Initially, a comparison was made between cancerous tissue and adjacent noncancerous tissue. The results showed that MMP-7, MYC, WNT-5A, and AXIN2 were upregulated in cancerous tissue compared to adjacent noncancerous tissue, but only MMP-7 and WNT-5A reached statistical significance (P < 0:05) ( Figure 4). Then, we compared the same four genes between viral positive and negative groups. MYC, WNT-5A, and AXIN2 were shown to be upregulated in HPV-positive CRC tissues when compared to HPV-negative tissues, while MMP7 was shown to be downregulated in HPV-positive CRC tissues when compared to HPV negative tissues, but the difference did not reach the level of statistical significance (P > 0:05) ( Figure 5).
3.6. DEGs' Principal Component Analysis. We then applied the hierarchical clustering analysis and principal component analysis (PCA) to categorize all the genes into different groups. Unsupervised two-dimensional hierarchical cluster analysis of relative gene expression data (2DCt method) reveals distinct clusters of cancer samples and healthy control samples. With PCA analyses, the 92gene signature discriminated four samples into two significantly separated areas: viral-positive and viral-negative groups. Our results revealed that isolated areas considerably accurately classify the patterns of viral-positive and viral-negative CRC ( Figure 6). The two-way hierarchical clustering was also used in the data from TCGA, and the results revealed significantly separated areas between CRC cancer tissue and healthy tissues (Figure 7). Our data show that multiple genes participate in the process of viral-related CRC.
Comparative analysis of the DEGs reported by independent studies shows a relatively limited degree of overlap due to variations in population and technical methods [36,37]. Therefore, it is still trying to delineate the underlying genetic events for individual CRC patients precisely. TaqMan Array Plates technology provides a highly effective way to analyze the expression changes of mRNA in CRC [38]. In this study, the DEGs detected by TaqMan Array Plates were related to colorectal carcinogenesis, cell cycle progression, invasion, antiapoptosis, cell adhesion and proliferation, and carbohydrate metabolism. Specifically, we found that 39 DEGs were upregulated, while 17 DEGs were downregulated when comparing HPV-positive CRC tissues with HPV-negative tissues. The DEGs in HPV-positive tissues were involved in signaling pathways including activators of Wnt/β-catenin, phosphoinositide 3-kinase (PI3K), mitogen-activated protein kinase (MAPK), mammalian target of rapamycin (mTOR) signaling, and inhibitors of TGFβ and TP53 signaling, which were also the most frequent changes in CRC progression. The DEGs are also involved in cell invasion, adhesion, and tissue remodeling which may play an essential role in the carcinogenesis and progression of metastases in CRC [34,[39][40][41]. With PCA analyses, we found that genes expressed in viral-  BioMed Research International positive and viral-negative samples have located in distinctly different areas, which might accurately classify viral-positiveand viral-negative-related CRC DEGs. Our data is partially overlapping with TCGA data. TCGA data showed that DEG expression among CRC, adenoma, and healthy tissue was different. DEGs in CRC tumor tissue have also significantly separated areas with genes expressed in healthy tissue. The Wnt/β-catenin signaling pathway controls cell proliferation and differentiation in the intestinal crypt microenvironment [42]. In the intestinal epithelium, terminally differentiated intestinal epithelial cells (IECs) were constitutively replaced by newly divided IECs from intestinal stem cells (ISCs) located in the crypts. This biological process was tightly controlled by Wnt signaling [43]. Deregulations of Wnt/β-catenin signaling were associated with CRC [44,45]. The APC protein, encoded by the APC gene, plays a vital role in maintaining homeostasis of the gut epithelium and thus, when mutated, can impart a functional consequence capable of driving cancerous cell proliferation [46]. In this study, APC and WNT5A were upregulated in the HPVpositive group as compared to the HPV-negative group, indicating that HPV infection may activate the Wnt/βcatenin signaling pathway. Evidence from the previous study also supports that microbes play a role in dysregulating the Wnt/β-catenin signaling pathway as a means of tumor initiation [47]. AXIN proteins act as tumor suppressors that function through inhibition of Wnt signaling [48]. AXIN2 is a primary transcriptional target of β-catenin-dependent Wnt signaling, which acts as a scaffold factor in the βcatenin destruction complex. The two AXIN proteins have similarity in dysregulation of the Wnt pathway that believed to be an essential first step in the genesis of CRCs, while aberrant Wnt pathway activation is in roughly 90% of sporadic CRCs [34,49].
In our study, we found that WNT5A was upregulated while AXIN2 was downregulated in the HPV-positive group 9 BioMed Research International as compared to the HPV-negative group. However, both AXIN2 and WNT5A were upregulated in the HPV-positive group when we verified the results with qPCR in additional samples. The data from TCGA showed that AXIN2 was upregulated in CRC tumor tissues. Our results revealed that differences in gene expression might be extended in individual patients with CRC; therefore, a large cohort is required to eliminate a distinctive difference in CRC oncogenes. The differences between microarray and q-PCR should be carefully examined to reveal the role of AXIN and Wnt signaling in HPV-associated CRC pathogenesis. Genes mutations in the Wnt/β-catenin signaling pathway include the c-MYC proto-oncogene (MYC) [50]. c-Myc is a target of APC signaling, which aberrantly activated in most colorectal cancers [51]. With the qPCR method, we found that MYC upregulated in all CRC cancerous tissues including HPV-positive cancerous tissues. MYC encodes a transcription factor that stimulates cellular proliferation and growth by controlling the expression of genes whose products regulate metabolism, ribosome biogenesis, and cell cycle progression [52].
Interferon (IFN) has a wide range of biological functions, including antiviral, antiproliferative, and immunomodulatory properties. IFN-induced cellular antiviral response is the primary defense mechanism against viral infections. Consequently, viral infection may disrupt the secretion of interferon, IFNG, INFGR1, and TNF genes upregulated in HPV-positive CRC tissues as compared to HPV-negative CRC tissues in this study. The upregulated IFNG and INFGR1 may trigger the STAT (signal transducers and acti-vators of transcription) proteins that mediate some biological functions and form part of the signaling cascades [53]. Therefore, the intimate relationship between viral infection, IFN/TNF/IL-6R/STAT1expression, and the associated signaling pathway is valuable in elucidating the mechanism of virus-associating CRC.
Striking experimental data have shown that gut microbes play a crucial role in vascularization of the intestinal mucosa and wound-healing processes [54]. MMPs and VEGFA are the most critical prometastatic and proangiogenic factors that inhibit cancer cell apoptosis, decrease cell adhesion, and induce angiogenesis, resulting in the promotion of the development and progression of CRC [55]. KRAS is a small GTPase, playing a significant role in the EGFR-MEKMAP kinase pathway. KRAS mutations are detectable relatively early in colon carcinogenesis, which is relevant to tumor progression and seems to be associated with distant metastases [56]. The PI3K/Akt signaling pathway affects diverse cell activities such as proliferation, differentiation, migration, and apoptosis, which is a crucial antiapoptosis pathway that regulates cell survival and death [57]. PI3K/Akt signaling also mediates inflammatory reactions by interacting with NF-κB [58]. PTEN is a molecular switch that antagonizes the PI3K/Akt signaling pathway by dephosphorylating PIP3 to PIP2, which suppresses cell proliferation and promotes apoptosis. Downregulation of PTEN results in an enrichment of PIP3 and a high expression of Akt [59]. In this study, we also found that MMPs, VEGFA, JAK1, EGF, KRAS, AKT1, MAPK3, and MAP2K4 were upregulated in

10
BioMed Research International HPV-positive CRC. The functions of these genes in virusassociated CRC need further investigation. Currently, traditional clinical and pathological parameters are not always sufficient to discriminate high-risk from low-risk CRC, and validated molecular markers with prognostic value are not yet available [33]. A recent report demonstrated that analysis of the microbiome could serve as a screening biomarker for CRC [60]. High-risk HPVs are present in 53.84% of invasive colorectal cases, which is associated with an invasive and metastatic phenotype [61]. E6/E7 of HPV type 16 induces cellular transformation and migration in human healthy colorectal mesenchymal cells but not epithelial ones, accompanied by the upregulation of D-type cyclins and cyclin E as well as Id-1 in these cells [62]. Our results revealed that the HPV antigen was detected more strongly in tumor tissue than in adjacent nonneoplastic tissue. Moreover, all adjacent noncancerous tissues that were viral negative at TNM stage I were changed into half of the adjacent noncancerous tissues that were viral positive at stage IV. Our results indicated that HPV might be associated with the CRC stage. When viral proteins begin to accumulate, progeny viral genomes are replicated and used for secondary transcription. HPV integration sites, with a preference for sites of known genomic fragility, are distributed randomly over the whole genome in one study, and the majority of integrated HPV genomes appear to be actively transcribed [63]. In this pilot study, DEGs involved various signaling pathways. A larger cohort to identify genes from HPV-associated CRC is necessary to obtain comprehensive knowledge of the signaling pathways associated with CRC carcinogenesis, discover novel targets for cancer therapy, and develop new biological drugs [33].
The main limitation of this study is that no significance was observed between deregulated genes from HPVpositive vs. HPV-negative CRC cases of the cohort. We could not exclude the reason for the small sample size. However, there may be a more complicated mechanism involved. Further study with large sample size is needed to elucidate the HPV-associated CRC.

Conclusions
HPV-induced cervical cancer has been widely confirmed clinically, and the anti-HPV vaccine has made significant progress in the prevention and treatment of cervical cancer. At present, studies have found that HPV viral antigens also 11 BioMed Research International exist in colon cancer tissues. The relationship between HPV infection and CRC carcinogenesis is very worth studying.
In this study, we demonstrated that tumor tissues had higher levels of HPV antigen than adjacent nonneoplastic tissue. Moreover, we found that 39 differentially expressed genes (DEGs) were upregulated, while 17 DEGs were downregulated when compared to HPV-positive CRC tissues with HPV-negative tissues. Combined with the database of the Cancer Genome Atlas (TCGA), we screen four DEGs and tested with the PCR method, and there was some overlap between our data and TCGA data. Our data revealed that HPV could take part in CRC carcinogenesis through gene mutation. The further genomic and proteomic investigation is necessary for obtaining a more comprehensive knowledge of signaling pathways associated with CRC carcinogenesis.

Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.

Ethical Approval
The study was approved by the Ethics Committee of Shantou University Medical College. The human investigations performed after were approved by an institutional review board.

Consent
Informed written consent was obtained from each subject or each subject's guardian.