1. IntroductionDespite improvements in early detection and treatment methods during the last two decades, colorectal cancer (CRC) remains the 3rd most common cause of cancer-related death worldwide with approximately 150,000 new CRC cases diagnosed in the United States annually [1]. Among them, approximately 20% of patients will have distant metastasis, and around 30% of patients with stage II and III disease will develop metastasis [2]. The five-year survival rate of CRC patients with distant metastasis is less than 15% [2]. The incidence of CRC in men and women under the age of 50 has steadily increased in the past two decades [1,2]. For patients with metastatic CRC (mCRC), chemotherapy remains the mainstay of treatment, but eventually, all patients develop resistance to therapy and experience treatment failure due to the intra-tumoral heterogeneity of CRC [3].Precision medicine guided by tumor genomic profiling has transformed the cancer diagnosis, prognosis, and treatment paradigm over the past two decades. However, recent estimates suggest that fewer than 10% of cancer patients benefit from this approach [4,5]. The primary issues facing genomic profiling include firstly, actionable genomic alterations are not detected in a vast majority of cases [6], and even when they are, secondly, a significant proportion of patients fail to experience an antitumor response to the indicated targeted therapy [7].One such target is the amplification of ERBB2 (HER2), which occurs in approximately 3% of patients with metastatic CRC (mCRC) and 5% of patients with wild-type NRAS and KRAS tumors [8,9]. Several ERBB2-targeted therapies are either in different phases of clinical trials or approved for use in patients with ERBB2-positive breast and gastric and gastroesophageal tumors [8]. However, despite recommendations of the National Comprehensive Cancer Network guidelines and clinical evidence from phase II trials that anti-ERBB2 therapies improve disease outcomes in ERBB2-positive mCRC patients, no ERBB2-directed approved therapies for patients with CRC are currently approved for clinical use [8,9,10].The role of HER2 in carcinogenesis is most well-characterized in breast cancer [11,12]. HER2+ breast cancer is a historically aggressive subtype of breast cancer with a five-year survival rate of 30% [11,13,14]. The discovery that amplification or overexpression of ERBB2 was associated with extremely poor survival in breast cancer led to efforts that resulted in the development of a monoclonal antibody (mAb) to HER2, trastuzumab [14,15]. However, whether ERBB2 overexpression-mediated carcinogenesis follows similar mechanisms in breast and colon tissue is unknown.Thus far, the majority of the studies related to ERBB2 in cancers have focused on identifying the landscape of genomic amplification in ERBB2 and defining therapeutic regimens to target these amplifications [8]. Our earlier study has shown that ERBB2 mRNA and protein overexpression correlates with more aggressive colorectal cancer in the North African population [16]. However, the effect of the overexpression of ERBB2 on the global transcriptomic profiles within CRC patients is not known. Therefore, the objective of the current study is to characterize whole transcriptomic changes associated with ERBB2 overexpression in CRC cell lines and patient samples with a view to gain a deeper understanding of the role of ERBB2 in CRC pathogenesis. 2. Materials and Methods 2.1. Patients and Tissue Specimens

Ethical approval for the study was provided by the Research and Ethics Committee (REC) of University Hospital Sharjah (UHS-HERC-055-25022019). All methods were performed in accordance with the relevant guidelines based on the Declaration of Helsinki and the Belmont Report. We obtained written informed consent from all study participants. This is a retrospective study of 14 patients with primary CRC. Patients with secondary cancers were excluded, whilst all primary CRC patients were included, regardless of age, gender, or tumor stage. The initial diagnosis was performed prior to and independently of our study to determine the Tumor, lymph Nodes, and Metastasis (TNM) score. Tissues were sectioned from formalin-fixed, paraffin-embedded (FFPE) biopsies for molecular and immunohistochemical analysis.

2.2. ImmunohistochemistryTo begin with, 3 μm sections from the FFPE of 14 CRC patients’ biopsies were immunohistochemically stained using the rabbit monoclonal antibody for HER2 (1:4000 dilution; ab214275, Abcam, Waltham, MA, USA) according to the manufacturer’s instructions. An experienced pathologist (R.H.) scored the stained slides, following the consensus recommendations for HER2 scoring for CRC [17,18]. Briefly, scoring was performed on a 4-point scale—0, 1+, 2+, 3+ focusing on intensity and extent according to the Allred scoring system [19]. In this study, 0 and 1+ intensity were taken to be negative for HER2 expression and the study focused on the assessment of membranous HER2 expression. 2.3. Cell Culture

The CRC cell lines HT29 and HCT116, and two normal colon cell lines: CCD33 and CCD841, were obtained from Bio Medical Scientific Services (BIOMSS, Al Ain, United Arab Emirates), and cultured in Dulbecco’s Modified Eagle’s Medium (DMEM) supplemented with 10% fetal bovine serum (Sigma Aldrich, St. Louis, MO, USA), 1% Penicillin/Streptomycin (Sigma) and 20 mM L-Glutamine (Sigma) at 37 °C in 5% CO2 incubator.

2.4. TransfectionThe wild-type ERBB2 expression construct was a gift from Mien-Chie Hung (Addgene plasmid #16257; https://www.addgene.org/16257, accessed on 15 December 2020) [20]. The cell lines were transfected with 5µg of pcDNA3-ERBB2 plasmid construct using Lipofectamine 3000 reagent (ThermoFisher Scientific, Cambridge, MA, USA) according to the manufacturer’s instructions. Cells transfected with the empty pcDNA3 vector served as the experimental control. The ERBB2 expression level was checked 24 h post-transfection at the mRNA and protein levels using qRT-PCR and Western blotting, respectively. 2.5. RNA Isolation

RNA extraction was carried out from three sequential (3 µm) sections from the same FFPE block. A needle macrodissection was carried out to enrich the tumor’s content. This was carried out by marking the tumor areas on the slides and carefully removing the unmarked non-tumor areas using a sterile needle, following which the marked areas were collected for molecular analysis. RNA was extracted using the RNA RecoverAll kit (ThermoFisher Scientific, Cambridge, MA, USA) according to the manufacturer’s instructions. Genomic DNA removal was ensured by treating the RNA with Turbo DNase (ThermoFisher Scientific). For RNA extraction from cell lines, cells were pelleted at 12 × 103 g for 5 min and rinsed thrice with ice-cold 1X sterile PBS. RNA was extracted from the cell pellet as described above using the RNA RecoverAll kit, followed by genomic DNA removal using Turbo DNase. All RNA samples were stored at −800 °C until further use.

2.6. Quantitative Reverse Transcriptase-PCR (qRT-PCR)qRT-PCR was performed using the Superscript First-strand Synthesis system (ThermoFisher Scientific). Real-time qPCR was performed in triplicates, using SYBR green (Solis BioDyne, Tartu, Estonia), on Quant Studio 3 (Applied Biosystem, Waltham, MA, USA). ERBB2 and the reference genes (18S ribosomal RNA) were pre-amplified using the following primer sets [18]. ERBB2_sense: 5′-ACATGCTCCGCCACCTCTACCA-3′; ERBB2 Antisense: 5′-GGATCTGCCTCACTTGGTTGTG-3′; 18SrRNA_sense: 5′-TGACTCAACACGGGAAACC-3′; 18SrRNA_antisense: 5′-TCGCTCCACCAACTAAGAAC-3′. A total of 40 PCR cycles were performed consisting of 15 s denaturation at 95 °C and a combined annealing and extension cycle of 10 min at 60 °C. The threshold cycle value (Ct) was normalized against the Ct value of internal control 18 s RNA. 2.7. Cell Lysis and Western BlotCell lysis and Western blot were performed as described previously [21]. Anti-HER2 (1:1000; Abcam) was used to probe the blots. All blots were subsequently stripped and re-probed for β-Actin (1:5000; Abcam, Waltham, MA, USA) to confirm equal loading. 2.8. Next-Generation RNA Sequencing

RNA sequencing was carried out on the indicated samples using a targeted AmpliSeq Transcriptome panel on Ion S5 XL System (ThermoFisher Scientific, Cambridge, MA, USA). In brief, ~30 ng of Turbo DNase treated RNA was used for cDNA synthesis using a SuperScript VILO cDNA Synthesis kit (ThermoFisher Scientific) followed by amplification using Ion AmpliSeq gene expression core panel primers. The enzymatic shearing was performed using FuPa reagent to obtain amplicons of ~200 bp and the sheared amplicons were ligated with the adapter and the unique barcodes. The prepared library was purified using Agencourt AMPure XP Beads (Beckman Coulter, Indianapolis, IN, USA) and the purified library was quantified using an Ion Library TaqMan™ Quantitation Kit (Applied Biosystems, Waltham, MA, USA). The libraries were further diluted to 100 pM and pooled equally with four individual samples per pool. The pooled libraries were amplified using emulsion PCR on Ion OneTouch™ 2 instrument (OT2) and the enrichment was performed on Ion OneTouch™ ES following the manufacturer’s instructions. Thus, prepared template libraries were then sequenced with Ion S5 XL Semiconductor sequencer (ThermoFisher Scientific, Cambridge, MA, USA) using the Ion 540™ Chip.

2.9. Bioinformatics AnalysesRNA-seq data were analyzed using Ion Torrent Software Suite version 5.4 and the alignment was carried out using the Torrent Mapping Alignment Program (TMAP). TMAP is optimized for aligning the raw sequencing reads against the reference sequence derived from the hg19 (GRCh37) assembly, and the specificity and sensitivity were maintained by implementing a two-stage mapping approach by employing BWA-short, BWA-long, SSAHA [22], Super-maximal Exact Matching [23] and the Smith–Waterman algorithm [24] for optimal mapping. Raw read counts of the targeted genes were performed using Samtools (Samtools view–c–F 4–L bed_file bam_file) and the number of expressed transcripts was confirmed after Fragments Per Kilobase Million (FPKM) normalization. For technical variations, code-set content normalization was performed with the geometric median for all genes. Principal component analysis (PCA) was performed using the indicated samples with R statistical software. Differentially expressed gene (DEG) analysis was performed using R/Bioconductor package DESeq2 with raw read counts from the RNA sequencing data [25,26]. Genes with less than ten normalized read counts were excluded from further analysis. A fold change of 2 was set as the cutoff for differentially expressed gene identification. p 2.10. Analyses of Publicly Available Transcriptomic Data Sets for Breast CancerIn order to compare the biological pathways and differentially expressed genes between the ERBB2 over-represented in breast cancer (BC) and CRC, transcriptomic data sets of BC were searched and retrieved from the GEO (Gene Expression Omnibus) database (https://www.ncbi.nlm.nih.gov/geo/, accessed on 3 October 2022). The datasets were searched based on “Breast cancer” and “ERBB2” keywords. Then, datasets including BC patients with variable ERBB2 expression based on the same platform Affymetrix Human Genome U133 Plus 2.0 Array were considered. Sixteen well-matched datasets were available, out of which fourteen were excluded for further analysis. Exclusion criteria were datasets performed in vitro cancer cell lines or in vivo study models using non-human species, repeated samples in super-series, and datasets exhibiting poor ERBB2 expression values uncharacteristic of HER2+/HER2−. Two datasets that met the criteria were selected, including GSE29431 and GSE48391 (Supplementary Figure S1). The raw data and the probe annotation files were downloaded for further analysis. 2.11. Breast Cancer Microarray Data AnalysisA total of 65 BC patients were selected in our analysis including 48 samples with low ERBB2 and 17 samples with high ERBB2 expression. The Affymetrix microarray represents more than 54,000 probes where each gene is represented with different probes. The raw data were processed and normalized using in-house R script as previously described [27]. For normalization and adaptive filtering, Affymetrix Microarray Suite 5 (MAS5) and Gene Chip Robust Multiarray Averaging (GCRMA) packages in Bioconductor/R software were applied. Probes with a MAS5 value > 50 and coefficient of variation (CV) 10–100% in GCRMA among all samples of each dataset were identified to get only common variant probes. The filtered probes were then annotated and collapsed into the gene names list based on the maximum expression of probes for each gene. The unchanged probes, positive control probes, and unassigned probes were excluded from the downstream analysis. The mapped gene expression lists were subjected to Gene Set Enrichment Analysis (GSEA) to identify the activated and enriched biological pathways between high and low ERBB2-BC patients. 2.12. GSEAGSEA was carried out separately for all resulting gene sets from the above different transcriptomic analyses including CRC patients, CRC cell lines, normal colon cell lines, and BC patients. First, the absolute GSEA was performed to identify the significantly enriched pathways among sets related to the C2: curated gene sets; C5: ontology gene sets including molecular function (MF) and biological process (BP); C6: oncogenic signature gene sets; and C7: immunologic signature gene sets. The results of the GSEA were ranked and selected according to the p27,28]. Next, the selected significant pathways were further analyzed to identify the differentially enriched genes and the leading edge genes in each pathway. In order to further reduce the set of resulting genes, a systematic cross-reference of each gene enriched within statistically significant pathways was carried out. Finally, the genes with the highest frequency across the multiple significant pathways enriched between the HER2 positive and HER2 negative samples were identified. 2.13. Statistical Analysis

Functional data are presented as mean ± SD, except where otherwise stated. When two groups were compared, the student’s t-test was used unless otherwise indicated. p < 0.05 was considered as statistically significant.

4. DiscussionOur earlier work has shown that HER2 overexpression is correlated with more aggressive disease in CRC patients [16], indicating that stratification of patients according to HER2 status might be beneficial in the early detection and subsequent therapeutic management of patients with metastatic CRC. The results from the current study reveal that HER2 overexpression is associated with distinct global transcriptomic profiling, the characterization of which might lead to the identification of putative diagnostic and prognostic biomarkers.

In our study, there were fewer up and downregulated genes in HT29 compared to HCT116, indicating that the transcriptomics perturbation post-HER2-overexpression was less in HT29 in comparison to HCT116. The expression levels of APC, GSK3B, and CTNNB1 were not significant, which indicates that perhaps there may be a link between ERBB2 expression and TP53 mutational status (rather than APC), since HT29 harbors TP53 mutation.

About half of the top 100 differentially expressed genes between HER2+ and HER2− CRC patients exhibited a more prominent difference in expression and were overexpressed in the HER2+ samples. A vast majority of them were non-coding RNA (small nucleolar RNAs, snoRNAs) that are normally involved in the biogenesis of other RNAs. That they were among the top differentially expressed genes indicates that expression of the snoRNAs may be regulated by HER2 and that they might be involved in the CRC pathogenesis in HER2+ patients. Indeed, it has been reported that numerous snoRNAs, including tumor-promoting and tumor-suppressing snoRNAs, are not only dysregulated in tumors but also show associations with clinical prognosis [30]. In addition, aberrant expression of snoRNAs has been reported in cell transformation, tumorigenesis, and metastasis, indicating that snoRNAs may be considered as biomarkers and/or therapeutic targets of cancer [31]. Even CRC associations between snoRNAs and CRC development have been reported [32,33,34]. For example, SNORA21 promotes CRC cell proliferation by regulating cancer-associated pathways such as Hippo and Wnt signaling pathways, and overexpression of SNORA21 has been reported to be associated with distant metastasis in CRC [35].Colorectal cancer demonstrates hyperactivation of the Wnt pathway, which is involved in tumorigenesis, stemness, and metastatic progression [36,37]. One of the enriched pathways in our analysis was the regulation of stem cell differentiation. HER2-overexpression in gastric cancer cells results in increased stemness and invasiveness [38]. Furthermore, this increased HER2-mediated stemness is regulated by Wnt/β-catenin signaling [39]. Our analysis also revealed significant enrichment of pathways related to positive and negative regulation of Wnt signaling and the enrichment of the Wnt signaling pathway genes WNT3A and WNT5A. These results would thus indicate that HER2 overexpression in CRC cells might result in poorly differentiated tumors that are more invasive.One of the leading-edge genes that were differentially expressed between HER2+ and HER2− patients and cell lines was a homeobox gene of the NKL subclass, cardiac transcription factor NKX2-5. NKX2-5 is one of the earliest known transcription factors required for cardiac cell specification and proliferation [40,41,42].NKX2-5 is expressed in several types of tumors [43,44,45,46], but its precise role in tumorigenesis is unknown. Another family member of the NKL, NKX2-1, has been reported to mediate p53-induced tumor suppression [47,48,49,50], Indeed, in the context of CRC, NKX2-5 functions as a conditional tumor suppressor gene via activating the p53-mediated p21WAF1/CIP1 expression [51]. It has been predicted via bioinformatic analysis and confirmed by chromatin immunoprecipitation analysis that in hepatocellular carcinoma cells the promoter of ERBB2 binds to the transcription factor NKX2-5, resulting in a negative regulatory effect [52]. Interestingly, promoter hypermethylation of NKX2-6 has been identified as a candidate biomarker associated with differential methylation in HER2+ breast cancer and breast carcinogenesis [53]. Whether a similar mechanism is operant for the HER2-mediated downregulation of NKX2-5 remains to be determined.During cardiogenesis, NKX2-5 potentiates Wnt signaling by regulating the expression of the R-spondin3 [54]. In the current study, we also saw the enrichment of pathways related to positive and negative regulation of Wnt signaling. Analysis of the RNAseq data showed significantly more NKX2-5 expression in HER2− patients compared to HER2+ patients (Supplementary Table S3). Whether HER2 overexpression drives negative regulation of NKX2-5-mediated Wnt signaling, ultimately resulting in well-differentiated less invasive tumors in CRC, remains to be determined.

It is important to note that the transcriptomic analysis in CRC cell lines was conducted post-transient transfection in CRC cell lines. It is highly plausible that short-term transfection (to an unphysiological expression level) in the CRC cell lines will induce transcriptomic changes that will be different from those observed in cell lines that have evolved under HER2-overexpression selective pressure. Hence, comparing our data to gene expression changes observed in HER2-amplified or HER2-mutant cell lines would indeed be informative in this context. However, this shortcoming is potentially offset to a large extent by three facts—(a) our analysis also involved HER2− and HER2+ patient samples, which are a better model compared to any of the homogeneous cell line models; (b) our usage of two different CRC cell lines and two normal colon cell lines; and (c) our observation that similar pathways were being enriched when the patient samples and the HER2-overexpressing cell lines were compared.

Despite the current study having inherent weaknesses in low sample size and lack of validation in a wider population, the results do highlight the importance of transcriptional profiling of HER2+ and HER2− CRC patients in identifying potential biomarkers that play a role in CRC pathogenesis. Additional studies are warranted in different population cohorts as the incidence of HER2+ CRC patients varies widely based on geographical location. It would be intricate and intriguing to investigate whether and how HER2 regulates NKX2-5 in wild-type and p53 mutant CRC and their subsequent effects on Wnt signaling and CRC invasives.