HORMAD2 methylation‐mediated epigenetic regulation of gene expression in thyroid cancer

Abstract This study is aimed to investigate the methylation level of candidate genes and its impact on thyroid carcinoma (THCA) development. Infinium Human Methylation 450 BeadChip Arrays by Illumina (Illumina HM450K) was the most popular CpG microarray platform widely used in biological and medical research. The methylation level of differentially expressed genes and their corresponding CpG sites were analysed by R programme. The expression of HORMAD2 was evaluated by qRT‐PCR and Western blot, while the methylation level was examined via methylation‐specific PCR. Cell viability, metastasis, cell cycle and apoptosis were detected by MTT assay, transwell and wound healing assay and flow cytometry, respectively, after treatment with 5‐aza‐2′‐deoxycytidine (5‐Aza). Tumour formation assay was used to analyse thyroid tumour growth in nude mice in vivo. The methylation levels of all 116 differentially expressed genes were analysed. HORMAD2 was significantly hypermethylated and its mRNA expression was inhibited in THCA cells. After treatment with 5‐Aza, HORMAD2 expression was up‐regulated in THCA cells and its overexpression can suppress thyroid cancer cell viability, mobility and invasiveness remarkably. Up‐regulation of HORMAD2 in THCA cells could prolong G0/G1 phase and shorten S phase to impede cell mitosis as well as promote thyroid cancer cells apoptosis. Furthermore, tumour formation assay showed that increased HORMAD2 level impeded tumour growth in vivo. Hypermethylation of HORMAD2 could induce THCA progression, while hypomethylation of HORMAD2 retard cell growth and mobility and facilitate apoptosis through increasing its mRNA expression.

incidence of thyroid carcinoma is still on the rise for uncertain reasons. Some researchers have supposed that it might be associated with epigenetic events. [3][4][5] As THCA recurrence is a common event in thyroid carcinoma patients, especially in the early stage. 6,7 Therefore, it is imperative to figure out the pathogenesis of THCA for developing effective diagnostic and therapeutic strategies.
DNA methylation is an epigenetic process linked to the regulation of several biological events, including transcriptional regulation of gene expression, X chromosome inactivation, genomic "imprinting," silencing endogenous retroviruses and so forth. 8 The alterations in DNA methylation status in cancer cells are characterized by promoter CpG island hypermethylation and diffuse genomic hypomethylation, which has emerged as an important mechanism involved in cancer progression. Hypermethylation of promoter regions of crucial genes can repress relevant gene expression and influence biological functions, thereby facilitating tumorigenesis. 9 For instance, Wehbe et al 10 12 Thus, aberrant promoter methylation may be an effective way to detect candidate targets for THCA treatment.
In recent years, the deeper insights into THCA and the rapid development of molecular detection technology have allowed analyses of THCA at the molecular level. Genome-wide expression analysis has been successfully adopted to identify molecular signatures, improving the diagnosis and prognosis of several types of tumours. 13 Considerable researches have explored the mechanism of certain genes in multiple types of cancers via genome-wide microarray analysis of gene expression profiling. For instance, Hou et al 14  profiles. HORMA-domain containing 2 (HORMAD2) was a conserved meiotic chromosomal protein-coding gene in many organisms, including yeast, plants, nematodes and mammals. 16 Nonetheless, few studies focused on the molecular mechanism of HORMAD2 from the perspective of whole genome DNA methylation and epigenetic regulation of gene expression.
The purpose of present research is to probe into the methylation level of candidate genes and the effects on THCA progression. We first identified differentially expressed genes (DEGs) and relevant CpG sites via bioinformatics analysis via R programme. Then, DNA methylation level in tumour and normal tissues was detected by methylation-specific PCR (MSP). Subsequently, we observed significantly hypermethylated HORMAD2, of which the expression was examined by qRT-PCR and Western blot. Furthermore, 5 0 -aza-2 0 -deoxycytidine was utilized to investigate the impact on the methylation level of HORMAD2 and biological functions in THCA.

| DNA methylation profiling
A portion of these DNA samples were isolated from thyroid cancer tissues and adjacent normal tissues in 507 patients (136 male and 371 female). The DNA methylation statuses of 108 DNA samples in total that were hybridized in the Infinium Human Methylation 450 BeadChip were evaluated, and then we integrated the data sets and removed probes that did not exist in the data sets. Furthermore, approximately 200 000 CpGs were abandoned according to SNPs, repeats and multiple mapping sites. The final set incorporated above 180 000 unique probes.

| CHARM (comprehensive high-throughput arrays for relative methylation)
Comprehensive high-throughput arrays for relative methylation was applied to McrBC analysis, the array design and computational algorithms are fractionation method-independent and make this a simple, general, relatively inexpensive tool suitable for genome-wide analysis, and in which individual samples can be assayed reliably at very high density, allowing locus-level genome-wide epigenetic discrimination of individuals, not just groups of samples. The process of establishing the array was as follows: (i) all the CpGs in the genome were identified and any region of 300 bp with no CpGs was removed. (ii) Probes with multiple matches including fuzzy matches as defined by NimbleGen were discarded to the genome. (iii) Any region of 300 bp between sequential probes was divided into two new regions, while any region with fewer than 15 probes was removed. (iv) Regions were tiled using 50-mers 35 bp apart.

| Data processing and bioinformatics analysis
The b values, or the percentage of CpGs at a given site that were methylated, were calculated for every sample at each CpG site.
Observations with detection P > .05 were set to missing, and any CpG site with missing data was omitted from the analysis. The R programme, included in the SVA package, was used to adjust for the specific chip. A parametric empirical Bayes framework, that is a particularly robust method for dealing with small samples sizes, was utilized to adjust data for batch effects. For each site and each individual, b values were calculated by the following formula: b = signal Meth /(signal meth + signal Unmeth ).

| Overall survival analysis
MethSurv performed univariable and multivariable survival analysis based on patient methylation levels for any CpG site (probe) using Cox proportional hazards models. In the univariable analysis, survival LIN ET AL.

| 4641
analysis was performed with probe methylation levels as explanatory variable and survival time as the response variable. In multivariable analysis, in addition to methylation status, clinical covariates such as age, sex and stage can be used. To assess differences in survival, methylation levels of patients were dichotomized into higher (the methylation b value higher than the cut-off point) and lower groups (the methylation b value lower than the cut-off point). Hazard ratio (HR) with 95% CI is derived from Cox fitting. The differences in survival between lower and higher methylated groups of patients were visualized by Kaplan-Meier survival curve.

| RNA isolation and qRT-PCR
Total RNA was isolated using Trizol reagent (Life Technologies), and then evaluated by gel electrophoresis and spectrophotometric analysis. First strand cDNA was generated through 5 mg of total RNA using the Superscript III-reverse transcriptase kit (Invitrogen) based on the instructions, followed by dilution of the reaction mixture to 100 mL with water. 2.5 mL of diluted cDNA mixture was employed for PCR reaction. The cycling condition was 95°C

| Methylation-specific PCR (MSP)
Bisulfite-treated DNA (50 ng) was mingled with AmpliTaq Gold polymerase (Applied Biosystems), mgCl2, and deoxynucleotide triphosphates for amplified reaction of MSP. The PCR amplification was carried out for 40 cycles at the annealing temperature of 60°C for 30 seconds. The final products were subjected to electrophoresis in a 2% agarose gel, and then recorded using a Molecular Imager (Bio-Rad, Hercules, CA, USA). The methylation-specific primers were listed in Table 2.

| 5-aza-2 0 -deoxycytidine treatment
Thyroid cancer cell lines were treated with 5-aza-2 0 -deoxycytidine (Sigma) to demethylate HORMAD2 gene. For demethylation, the cells were cultured in the growth medium, to which 5-aza-2 0 -deoxycytidine was added at a concentration of 2 mmol/L for 3 days before extraction of DNA using the Puregene DNA isolation kit (Gentra Systems Inc., Munich, Germany).

| Western blot
A total of 20 lg proteins were added in 10% sodium dodecyl sulfate-polyacrylamide gel, purified by centrifugalization, separated by electrophoresis and then shifted onto nitrocellulose membranes. The membranes were sealed with 5% skim milk for 1 hour and incubated with antibodies: rabbit anti-HORMAD2 (ab89961, Abcam, Cambridge, UK, dilution: 1:5000) overnight at 4°C and then with HRP-conjugated mouse anti-rabbit secondary antibody (1:2500) for 1 hour at 37°C. After washing with Tris-Buffered Saline Tween-20 (TBST) and Tris-Buffered Saline (TBS), the protein bands were visualized via enhanced chemiluminescence (ECL) kit (Bio-Rad). b-actin was used as an internal control.

| Flow cytometry analysis
The cell cycle and apoptosis were detected by flow cytometry. For apoptosis, both adherent and non-adherent cells were collected T A B L E 1 qRT-PCR primer sequence

| Statistical analysis
The experimental data were analysed by Graphpad 6.0 (Version 6, CA, USA). Results were displayed in form of mean AE standard deviation. Comparison between the two sets of data was performed using Student's t test, while the differences between the three or more data were assessed by One-way ANOVA. All experiments were repeated thrice. P < .05 signified a statistical significance.

| Genome-wide methylation data for thyroid cancer tissues
Multi-dimensional scaling most 1000 variable positions showed that tumour and normal tissues had significantly different methylation level ( Figure 1A). The density plot used in methylation analysis illustrated conspicuous difference between type I probes (blue line); type II probes (orange line) ( Figure 1B). Each sample analysed displayed significant difference of methylation level between normal and tumour groups ( Figure 1C). Unsupervised hierarchical clustering analysis was based on DNA methylation for 395 735 probes in control and tumour tissues ( Figure 1D). The heatmap of hierarchical clustering analysis represented CpG sites methylation levels from completely methylated (dark blue) to unmethylated (green) ( Figure 1E). These results indicated that the sample data we used exerted differential methylation level from two endpoints.

| The methylation level of top 1000 differentially methylated CpG sites in imprinted genes
Probes were implemented within the Bioconductor package

| HORMAD2 was hypermethylated in thyroid cancer tissues
The differentially methylated region DMR_19 enriched with HOR-MAD2 was significantly hypermethylated in thyroid cancer tissues ( Figure 3A) and the heatmap displayed top 20 different methylation genes enriched in THCA tissues, in which HORMAD2 was hypermethylated ( Figure 3B). ALL 116 genes enriched by DMR-related CpGs had different methylation level, which indicated that HORMAD2 was significantly hypermethylated in THCA tissues ( Figure 3C). Survival prognosis curve revealed that high methylation level of HORMAD2 was correlated with poor prognosis in THCA patients ( Figure 3D).
Based on the methylation levels, RFS of 54 THCA patients was a continuous variable, which could be considered as prognostic factors. RFS referred to the time between initial diagnosis and recurrence or death, with follow-up censored at last contact if no event had occurred. The methylation level of top 9 CpGs located in HORMAD2 were shown in Boxplot for cg01141459, cg04046669, cg13245431, cg14509403, cg15209808, cg16686158, cg17632937, cg2184594 and cg21890667, respectively ( Figure 4A-I). The above results illustrated the methylation level in THCA group was remarkably higher compared with control group.

| The expression of HORMAD2 was downregulated by hypermethylation in thyroid cancer
The expression of HORMAD2 was detected by qRT-PCR and Wes-