Microarray-based differential expression profiling of long noncoding RNAs and messenger RNAs in formalin-fixed paraffin-embedded human papillary thyroid carcinoma samples.

BACKGROUND: Long noncoding RNAs (lncRNAs) can regulate the expression of genes at almost every level. The altered expression of lncRNAs was observed in many kinds of cancers. Until recently, few studies have focused on the function of lncRNAs in the context of papillary thyroid carcinoma (PTC).
METHODS: In the current study, we collected seven PTC and nodular goiter tissue samples and explored mRNA and lncRNA expression patterns in these samples by microarray.
RESULTS: We observed aberrant expression of 94 lncRNAs and 99 mRNAs in the seven PTC samples as compared to the nodular goiter tissue [fold change (FC) ≥2.0; P<0.01]. To confirm these microarray results, quantitative polymerase chain reaction (q-PCR) was performed to assess the expression of three randomly selected differentially expressed mRNAs and lncRNAs, confirming our microarray findings significantly. We then performed gene ontology (GO) and kyoto encyclopedia of genes and genomes (KEGG) analyses to systematically characterize the twelve significantly differential genes. A co-expression analysis revealed that the lncRNAs n382996, n342483, and n409114 were closely related to the regulation of MT1G, MT1H, and MT1F.
CONCLUSIONS: In the present study a string of novel lncRNAs associated with PTC were identified. Further study of these lncRNAs should be performed to identify novel target molecules which may improve diagnosis and treatment of PTC. 2019 Translational Cancer Research. All rights reserved.

Introduction

In recent decades, the incidence of thyroid cancer has increased in most regions of the world (1,2). Among different forms of thyroid cancer, papillary thyroid carcinoma (PTC) is the most common form, accounting for 85–95% of thyroid tumors (3,4). Several population-based studies have revealed that while rates of thyroid cancer have increased, small tumors are now more likely to be detected and there have been no significant changes in mortality associated with this disease (5,6). A family history of thyroid cancer, exposure to radiation, and abnormal intake of iodine are all PTC risk factors (7). Recent studies have indicated that genomic instability, alterations in epigenetic regulation, and subsequent inappropriate gene expression have key roles in the regulation of PTC (8,9). However, the mechanisms underlying PTC pathogenesis have not been fully elucidated to date. In clinical contexts, nodular goiters are known to be capable of developing into poorly differentiated thyroid carcinoma (PDC), but the genomic mechanisms governing this transformation remain uncertain. A thorough investigation of the mechanisms underlying this transformation could thus offer new clinical insights into the pathogenesis or treatment of PDC.

The majority of genes are regulated through complex epigenetic mechanisms, with a diverse array of noncoding RNAs (ncRNAs) participating in these regulatory processes (10,11). In the context of PTC, the functions of specific microRNAs have been relatively well elucidated, but alterations in PTC-associated long noncoding RNAs (lncRNAs) have also been found to lead to the altered regulation of key genes, and these events are less well understood (12,13). Studies have shown that aberrant gene expression and the tumorigenesis can be caused by dysregulated expression of lncRNAs. In addition, in many kinds of cancers, changes in lncRNA expression can accelerate the invasion and metastatic progression of the disease (14). At present, the study on the role of lncRNAs in PTC progression is in its early stages, and further studies are needed to explore the lncRNAs associated with PTC and their mechanisms of action.

Intracellular mammalian metallothioneins (MTs) are low molecular weight proteins (~6 kDa) initially isolated from horse renal cortex tissue in 1957 by Margoshes and Valee (15,16). In humans, of the 17 different genes encoding MTs, 13, 2, 1, and 1 gene each encode MT-1, MT-2, MT-3, and MT-4, respectively (17-19). MT is primarily located in the cytoplasm of adult tissues, while in fetal/early neonatal tissues it is also found in the nuclei of cells (20). MT located in the nucleus has been observed during cell proliferation and differentiation (21). Expression of MT family genes and altered localization of MT in cells during development suggest that MT may be a proto-oncogene. The carcinogenic potential of MT has been evaluated in workshops and symposiums (22), and several reports have studied MT expression in human tumors (23).

In the current study, the expression of lncRNAs and mRNAs were detected in seven PTC and seven nodular goiter samples via a microarray analysis, and differentially expressed lncRNAs and mRNAs associated with PTC were identified. In addition, we demonstrated that genes associated with these lncRNAs and mRNAs were known or suspected risk factors for tumorigenesis. These differentially expressed lncRNAs and mRNAs may thus be novel molecular biomarkers useful for the diagnosis and/or treatment of PTC.

Methods

Patient samples and hematoxylin & eosin (H&E) staining

Seven patients with a follicular variant of papillary thyroid carcinoma (FVPTC) (2 males and 5 females, mean age: 39.43±13.06 y) and seven patients with nodular goiter, none of whom had undergone radiation, chemotherapy, or immunotherapy treatments, were recruited from Beijing Haidian Hospital Beijing Haidian Section of Peking University Third Hospital in China between January 2011 and January 2013. The diagnosis of PTC was made by at least two experienced oncologists. Sample collection was conducted according to the following criteria: the minimum tumor diameter was greater than 2 cm. The Formapure nucleic acid extraction kit (Agencourt Biosciences, Beverly, MA, USA) was then used to extract and purify the RNA from formalin-fixed, paraffin-embedded (FFPE) tissue sections prepared from these tumor samples. TRIzol (Invitrogen, Carlsbad, CA, USA) was used to extract RNA from frozen tissues and the RNeasy Protect kit (Qiagen, Valencia, CA, USA) was used to purify RNA from these samples. To remove contaminating DNA from RNA samples, DNase I treatment (Ambion, Austin, TX, USA) was additionally used. RNA concentrations were determined using a Nanodrop ND-1000 spectrophotometer (Nanodrop Technologies, Rockland, DE, USA), and electropherograms were used to evaluate the integrity of RNA, with the RNA 6000 PicoAssay for the Bioanalyzer 2100 (Agilent Technologies, Santa Clara, CA, USA) being chosen to determine the RNA integrity number, RIN42 (a correlative measure of mRNA intactness). Hematoxylin and eosin (H&E) staining was conducted as previously described. Briefly, FFPE tissue samples were sectioned at 3–5 µm thickness and stained with stained with hematoxylin (Sigma H3136) for 10 min and with eosin (Sigma E4382) for 1 min to better visualize diagnostic areas.

Microarray and computational analysis

Sample labeling and array hybridization were performed according to provided protocols. Briefly, after the rRNA was removed, the mRNA-ONLY™ Eukaryotic mRNA Isolation kit (Epicentre, Madison, WI, USA) was used to purify mRNA from total RNA. Subsequently, each sample was amplified and transcribed into fluorescent cRNA along the entire length of the transcripts without 3' bias utilizing a random priming method (Flash RNA Labeling kit; Arraystar, Rockville, MD, USA). The RNeasy Mini kit (Qiagen, Inc., Valencia, CA, USA) was then used to purify the labeled cRNAs. A NanoDrop ND-1000 was used to measure concentrations of the labeled cRNAs and their specific activity (pmol Cy3/µg cRNA). Labeled cRNAs were then hybridized onto the Affymetrix GeneChip® Human Transcriptome Array 2.0. lncRNAs were carefully constructed using the following public transcriptome databases: RefSeq (NCBI Reference Sequence Database; www.ncbi.nlm.nih.gov/refseq/); NCBI unigene (https://www.ncbi.nlm.nih.gov/unigene); Ensembl (http://ensemblgenomes.org/); UCSC (genome.ucsc.edu/index.html); LNCRNA-DB (http://www.lncrnadb.org/); as well as landmark publications (24). To identify each transcript accurately, a specific exon or splice junction probe was applied, and subsequently the hybridized arrays were washed and fixed. To analyze acquired array images, the Affymetrix DNA Microarray Scanner and Feature Extraction software was used to scan the fixed hybridized arrays. Subsequently, a variation of the reads/Kb/Million method and Z-score analysis were performed to normalize the raw data. Differences between microarray data for the two groups were compared via Student’s t-test. Significance of these microarray results was analyzed based on fold change (FC), and P values were corrected based on a false discovery rate (FDR)]. The threshold value used to define differentially expressed lncRNAs and mRNAs was a FC of ≥2.0 or ≤0.5, with P<0.05 as the threshold of statistically significance. SPSS v13.0 (SPSS, Inc., Chicago, IL, USA) was used for all statistical analyses.

Gene ontology (GO) and pathway analysis

To determine the functions of differentially expressed mRNAs and their associated pathways, we used the local FunNet software to analyze the GO term and pathway enrichment based on the available Gene Ontology platform (www.geneontology.org), which provides three structured networks of defined items: biological processes, cellular components, and molecular functions. In this analysis, the P value denotes the significance of enrichment of a given GO term in the differentially expressed mRNA list (cut-off: P<0.05). To highlight important pathways, pathways enriched in differentially expressed mRNAs were analyzed based on the latest kyoto encyclopedia of genes and genomes pathway annotations (www.genome.jp/kegg/). For the KEGG analysis, P values denote the significance of the pathway associated with the latest data available in the KEGG database (cut-off, P<0.05).

Hierarchical clustering analysis and co-expression analysis

Hierarchical clustering was performed to analyze the significantly differentially expressed lncRNA and mRNA datasets as previously described by Eisen et al. (25). Briefly, a log2 transformation was performed for the relative ratios of significantly aberrantly expressed lncRNAs and mRNAs, and then similarity was defined based on the Euclidean distance similarity metric, with the complete linkage clustering method chosen to assemble hierarchical clusters. The Search Tool for the Retrieval of Interacting Genes/Proteins (STRING) database of physical and functional interactions (26) was then used to analyze protein-protein interactions (PPIs), and the Cytoscape software (version 3.0.0) was used to represent the resultant PPI network (27).

Reverse transcription-quantitative polymerase chain reaction (RT-qPCR)

Differential expression of select mRNAs identified in the initial microarray was verified by RT-qPCR. In total, three mRNAs (down-regulated: MT1G, MT1F, and MT1H) were chosen for a q-PCR reaction in twenty PTC and nodular goiter samples. TRIzol was used to extract total RNA from forty FFPE specimens, and this RNA was reverse transcribed into cDNA using the PrimeScript RT reagent kit (Invitrogen), according to provided instructions. Samples were also treated with DNase (1 µL; 1U/µL) to remove contaminating DNA. qPCR reactions were then carried out in a total of 50 µl volume per well. qPCR settings were as follows: 95 °C for 2 min; 95 °C for 15 sec, 57 °C for 15 sec, and 72 °C for 1 min for a total 40 cycles. The Taqman probes used for the three mRNAs were as follows: ACTB: 5'-CCTTTGCCGATCCGCCGCCCGTCCA-3'; MT1G: 5'-GAAGTGCAGCTGCTGCGCCTGATGT-3' and 5'-GCTCCTGTGCCGCTGCAGGTGTCTC-3'; MT1F: 5'-TTATACCACCTTGACCCATTTGCTA-3'; and MT1H: 5'-TACAACTCCGACTCATTTGCTACAT-3'. Data were collected using the ABI analytical thermal cycler.

Results

Differential expression profiles of lncRNAs and mRNAs in PTC

shows the typical pathological characteristics of patients with PTC. PTC was identified in patient samples via H&E staining, and the ratio of cancer cells to total cells in individual tissue sections was at least 80%. PTC tissue and nodular goiter tissues in each patient were collected, fixed, and paraffin embedded. The human PTC and nodular goiter tissue samples were then analyzed by microarray to explore the differential expression of lncRNAs and mRNAs between these two sample types, with goiter samples serving as the experimental control tissue. The normalized expression values from the nodular goiter tissue samples were used to calculate the FCs and P values for each gene. Relative to the nodular goiter tissue samples, there were significantly differences in the expression of 94 lncRNAs and 89 mRNAs in PTC samples (FC ≥2.0; P<0.01). Among the 94 differently expressed lncRNAs, 15 were up-regulated and 79 were down-regulated in PTC samples. Among the 89 abnormally expressed mRNAs, 12 were up-regulated and 77 were down-regulated (, http://fp.amegroups.cn/cms/tcr.2019.02.12-1.pdf). Based on the differential expression levels of the identified lncRNAs and mRNAs, both were arranged into groups via a hierarchical clustering analysis (). The expression modes of lncRNA and mRNAs were further explored using a dendrogram (FC ≥2.0; P<0.01; ). n340297 (FC =3.7) was the most down-regulated lncRNA, while n332905 (FC =2.7) was the most up-regulated lncRNA. MT1G (FC =10.9) was the most down-regulated mRNA, while SLC34A2 (FC =4.3) was the most up-regulated mRNA. In total we have summarized 40 significantly differential expressed items, including 20 significantly up-regulated expressed items and 20 significantly down-regulated expressed items (). Notably, multiple members of the MT gene family were all down-regulated in PTC samples as compared with nodular goiter tissue samples ().

Figure 1

H&E staining of patient PTC and nodular goiter tissue samples. (A) The typical pathological characteristics of patients with nodular goiter; (B) the papillary structure of patients with PTC; (C) 1, nuclear groove; 2, ground glass nuclei; 3, intranuclear inclusions. H&E, hematoxylin & eosin; PTC, papillary thyroid carcinoma.

Figure 2

Differentially expressed gene identification. Red dots on the right indicate genes that are significantly up-regulated, while those on the left indicate genes that are significantly down-regulated. The X-axis represents the fold change of genes in PTC samples relative to nodular goiter samples. The Y-axis represents the statistical significance of changes in gene expression. PTC, papillary thyroid carcinoma.

Figure 3

Clustering analyses of significantly differential expressed mRNAs and lncRNAs. Hierarchical clustering of significantly differential expressed mRNAs and lncRNAs in PTC and nodular goiter patient samples. Each column represents one sample. Each row displays the changes of significantly differential expressed mRNAs and lncRNAs using color coding based on the relative ratio listed in Supplementary file 1. lncRNAs, long noncoding RNAs; PTC, papillary thyroid carcinoma.

Table 1

Forty significantly differential expressed items

Primary ID	Fold changes	Regulation	P value	Locus type
MT1G	10.94924	Down	6.45E-04	Coding
MT1H	7.715949	Down	7.24E-04	Coding
MT1F	7.657472	Down	0.001268	Coding
TPO	7.596827	Down	0.002575	Coding
SLC2644	4.001563	Down	0.003759	Coding
IYD	3.731338	Down	0.001858	Coding
n340297	3.698557	Down	0.007797	Non-coding
n407932	3.568332	Down	0.004669	Non-coding
n342369	3.523873	Down	5.32E-04	Non-coding
n336033	3.39254	Down	0.008311	Non-coding
n334406	3.364263	Down	0.007475	Non-coding
LRP2	3.35106	Down	0.00515	Coding
HSPA5	3.230271	Down	0.009589	Coding
n344917	3.22871	Down	0.003499	Non-coding
DIO1	3.213385	Down	0.008882	Coding
n339733	3.177889	Down	9.03E-04	Non-coding
n336012	3.147624	Down	0.003718	Non-coding
n340776	3.144596	Down	0.006244	Non-coding
ATP5H	3.117614	Down	0.002836	Coding
nN335616	3.087649	Down	0.004551	Non-coding
uc021yst.1	6.794604	Up	0.009216	Coding
SLC34A2	4.349818	Up	0.002504	Coding
ENST00000408542	3.050107	Up	0.004546	Non-coding/miRNA
uc021vey.1	2.737648	Up	0.004626	Coding
uc022avi.1	2.737648	Up	0.004626	Coding
n332905	2.712229	Up	6.49E-04	Non-coding
DCSTAMP	2.650783	Up	0.006883	Coding
uc021wsr.1	2.580299	Up	0.007814	Coding
ENST00000517241	2.471056	Up	0.007653	Non-coding/snRNA
n332412	2.446805	Up	2.29E-04	Non-coding
ENST00000408429	2.311055	Up	0.007908	Non-coding/miRNA
ENST00000401335	2.259779	Up	0.008486	Non-coding/miRNA
ENST00000517054	2.25374	Up	0.008552	Non-coding/snRNA
MIR412	2.230374	Up	0.003681	Non-coding/miRNA
ENST00000516485	2.22826	Up	0.006558	Non-coding/snRNA
CAMK2N1	2.21727	Up	9.50E-04	Coding
DTX4	2.213393	Up	0.004628	Coding
ENST00000408375	2.155374	Up	0.008106	Non-coding/miRNA
uc022bin.1	2.138922	Up	0.006837	Coding
ENST00000408745	2.110623	Up	0.009567	Non-coding/miRNA

Table 2

Members of the MT gene family that down-regulated in PTC samples

Gene symbol	Fold changes	Regulation	P value	Locus type
MT1G	10.94924	Down	6.45E-04	Coding
MT1H	7.715949	Down	7.24E-04	Coding
MT1F	7.657472	Down	0.001268	Coding
MT1L	2.810005	Down	0.001391	Non-coding/non-coding RNA
MT1JP	2.789976	Down	5.35E-04	Non-coding/non-coding RNA
MT1X	2.528777	Down	0.006529	Coding
MT1M	2.348179	Down	7.47E-04	Coding
MT1CP	2.328058	Down	9.09E-04	Non-coding/pseudogene
MT1E	2.255994	Down	0.002713	Coding
MT1A	2.147448	Down	6.96E-04	Coding
MT1B	1.49531	Down	0.001379	Coding

PTC, papillary thyroid carcinoma.

GO and pathway analysis

GO terms were further assigned to the significantly differentially expressed genes based on their sequence similarities with known proteins in the UniProt database. Genes were annotated both with GO terms as well as with the InterPro and Pfam domains contained therein. GO annotation and enrichment analysis of twelve differentially expressed sequences including MT1G, MT1H, MT1F, TPO, uc021yst, SLC34A2, SLC26A4, IYD, n340297, n407932, n342369 and n336033 was conducted using the local FunNet software (28), which corrected for gene length bias. GO terms with corrected P values <0.05 were considered to be significantly enriched among these differential expressed genes (http://fp.amegroups.cn/cms/tcr.2019.02.12-2.pdf, http://fp.amegroups.cn/cms/tcr.2019.02.12-3.pdf, http://fp.amegroups.cn/cms/tcr.2019.02.12-4.pdf). The bar charts of GO biological processes, GO cellular components, and GO molecular functions distributions are shown in . With respect to molecular functions, 96 GO terms were assigned among which the top three significantly enriched terms were protein binding (GO: 0005515, P value =6.14E-22), RNA binding (GO: 0003723, P value =4.35E-05), and structural constituent of ribosome (GO: 0003735, P value =3.76E-12). With respect to cellular components, 70 GO terms were assigned of which the most significantly enriched terms were cytoplasm (GO: 0005737, P value =5.65E-03), membrane (GO: 0016020, P value =1.86E-03), and cytosol (GO: 0005829, P value =8.28E-13). With respect to biological processes, 237 GO terms were assigned among which gene expression (GO: 0010467, P value =7.93E-16) and metabolic process of cellular protein (GO: 0044267, P value =1.86E-19) were the top two significantly over-represented terms ().

Figure 4

GO biological process-enriched categories of the twelve most differentially expressed mRNAs and lncRNAs. The bar corresponds to the percentage of differentially expressed genes in relation to all annotated genes in the respective category. GO, gene ontology; lncRNAs, long noncoding RNAs.

In vivo, biological functions arise from complex and interdependent regulation of multiple different signaling pathways. Pathways enrichment analysis can offer some clues with respect to the biochemical and signal transduction pathways in which differentially expressed genes may participate. The KEGG database can thus be used to understand the high-level functions of a given biological system, such as a cell, an organism, or an ecosystem using large-scale molecular datasets generated by genome sequencing and other high-through put experimental technologies (http://www.genome.jp/kegg/). The local FunNet software was used to test the statistical enrichment of differentially expressed genes in particular KEGG pathways. In this study, 11 differentially expressed genes were found to be involved in 18 pathways (http://fp.amegroups.cn/cms/tcr.2019.02.12-5.pdf). shows the results of the pathway enrichment analysis, clearly revealing that protein processing in endoplasmic reticulum (04141, P value =4.3E-14) was the top enriched item. Eight differentially expressed genes identified in our study participate in this pathway. In addition, oxidative phosphorylation, Huntington’s disease, and Alzheimer’s were also significant enriched in this study.

Figure 5

KEGG analyses of the twelve most differentially expressed mRNAs and lncRNAs. The bar corresponds to the percentage of differentially expressed genes in relation to all annotated genes in the respective category. KEGG, kyoto encyclopedia of genes and genomes; LncRNAs, long noncoding RNAs.

RT-qPCR and co-expression analysis

To verify the microarray data, three significantly down-regulated mRNAs (MT1G, MT1F, and MT1H) mRNAs were randomly selected for an RT-qPCR analysis in 20 PTC and 20 nodular goiter tissue samples. The results of this analysis were consistent with the microarray data (). Moreover, to identify the regulatory networks dictating relationships between the differentially expressed lncRNAs and mRNAs and their potential substrates a co-expression analysis was performed using the STRING database. To improve the accuracy of this analysis, we set the confidence level (score) to a high value (0.900). Then, using the Cytoscape software, the network of all lncRNAs and mRNAs with their potential substrates was extracted from the whole interaction network and reconstructed. Based on Pearson correlation values (Pearson correlation >0.90, P value <0.01), we further screen out 57 lncRNAs and mRNAs which participated in the regulatory network. Of note, n382996, n342483 and n409114 were closely associated with MT1G, MT1H and MT1F. Based on these results, we speculate that these lncRNA may participate in the regulation of the MT gene family during PTC development or progression ().

Figure 6

q-PCR analysis of MT1G, MT1F, and MT1H expression. q-PCR, quantitative polymerase chain reaction. ***, P<0.001.

Figure 7

Co-expression analyses of mRNAs and lncRNAs. Interactions of mRNAs and lncRNAs were extracted by searching the STRING database using a confidence cutoff value >0.900 and P value <0.01. The interaction network was then reconstructed using the Cytoscape software. The yellow squares, the wathet blue circles, and the green circles represent lncRNAs, mRNAs and MT mRNAs, respectively. Only the significantly differential expressed mRNAs and lncRNAs identified in this study were used for construction of this network. LncRNAs, long noncoding RNAs.

Discussion

PTC is a common endocrine cancer that occurs more often in women (29). The incidence of PTC has been increasing worldwide in recent decades. PTC tumorigenesis is an intricate process with many differentially expressed genes contributing to the development of disease. Many studies have extensively explored the molecular mechanism of PTC. However, the exact pathogenesis of this disease is not fully understood. In clinical settings, nodular goiter tissues can develop into PDC, but the genomic mechanisms underlying this transformation are poorly understood. Studies of gene expression in thyroid tumors offer insights into the molecular differences among different thyroid disease types. In a valuable study performed by The Cancer Genome Atlas Research Network, the genomic landscape of PTC has been comprehensively described, revealing that compared with other solid tumor types, PTC has a lower frequency of somatic mutations. CHEK2, ATM, TERT, miR-21, and miR-146b all contribute to the clinical classification and progression of this disease. Eszlinger et al. have described the differences in gene expression between the functional and nonfunctional benign thyroid tumors (30). In a separate study performed by Huang et al. (31) using high density DNA microarrays to analyze PTC samples, certain markers identified in this study were validated as PTC markers (32). Recent evidence strongly suggests that lncRNAs play a vital role in regulating gene expression. Although lncRNAs were initially considered to be a form of transcriptional noise, studies have more recently revealed their significant role in cellular development and various diseases (33,34), including a role for lnc-IL7R in inflammation (35), and for HOTAIR and Gas5 in breast tumors and metastases (36,37). However, in the pathogenesis and development of PTC, the role of lncRNAs is not well understood. Our present study revealed differentially expression profiles for both lncRNAs and mRNAs in PTC tumor samples, comparing the expression level of these transcripts with those in nodular goiter tissue. In total, 94 lncRNAs and 89 mRNAs were significantly differentially expressed in PTC samples. Among the 94 lncRNAs, 15 were up-regulated and 79 were down-regulated. Among the 89 mRNAs, 12 were up-regulated and 77 were down-regulated ( and http://tcr.amegroups.com/public/system/tcr/supp-tcr.2019.02.12-1.pdf). A dendrogram revealed the lncRNA and mRNAs expression patterns among samples (FC ≥2.0; P<0.01; ). n340297 (FC =3.7) was the most down-regulated lncRNA, while n332905 (FC =2.7) was the most up-regulated lncRNA. Additionally, MT1G (FC =10.9) was the most down-regulated mRNA and SLC34A2 (FC =4.3) was the most up-regulated mRNA.

The results of the present study revealed an altered lncRNA and mRNA expression profile in PTC. In our study, we found that the genes related to our observed differentially expressed lncRNAs and mRNAs were associated with many kinds of cancer. These genes included MT1s, which were significantly down-regulated in PTC and in several other types of cancer. A study found that among 34 human thyroid tumors 31 showed either nuclear and/or cytoplasmic localization of MT, while less than 20% of the normal thyroid gland tissue exhibited diffuse MT distribution (38). Compared to normal thyroid, down-regulated MT-I+II genes have been found both in papillary and follicular thyroid carcinomas by microarray analysis (39). When MT-1s expression was restored to malignant cells via cDNA transfection, cells growth rates and tumorigenicity were suppressed in vivo (39), indicating a suppressive role for MT-1s in thyroid carcinogenesis. This was consistent with a previous study which found that MT-1s was down-regulated in PTC samples relative to normal thyroid tissue owing to the hypermethylation of the MT-1s gene promoter. Furthermore, a negative correlation was detected between MT-1s and poor thyroid cancer prognosis in patients with other poor prognostic factors (40). MT-1s is also vital to many cellular processes, regulating the homeostasis of zinc ions and thus potentially regulating the activity of multiple zinc-dependent transcription factors and enzymes, such as metalloproteinase (41,42). There is abundant evidence that MT-1s plays a role in cell proliferation—MT-1s can be observed in the nuclei of cells in certain stages of the cell cycle, whereas MT-1s is primarily nuclear in quiescent cells (43). Biochemical studies have found that the concentrations of MT-1s are highest in the cytoplasm in the G1 phase, while MT-1s translocates into the nucleus during the G1/S phase (21). Furthermore, in the S and G2 phases the concentration of MT-1s was highest, supporting a role for this protein in cell division (44). Cytoplasmic-nuclear expression of MT-1s has been noted in regenerating kidney (45) and liver tissue (45,46), as well as in proliferating basal and parabasal cells of stratified epithelia (47). The role of MT-1s in the activation of cell proliferation is also supported by studies that found a positive correlation between MT-1s and the expression Ki-67 antigen in various tumors, including in breast cancer (48), squamous and basal cell skin cancers (49), adrenocortical cancers (ACC) (50), non-small cell lung cancers (NSCLC) (51), gastrointestinal stromal tumors (GISTs) (52), and sarcomas (53). Many studies have described the important role of MT-1s in protecting DNA against free radical damage (54) induced by radiation (55) or chemotherapeutic agents, such as cisplatin, etoposide, melfalan, bleomycin, or irinotecan administered during the treatment of breast, ovarian, gastric, or prostate cancer (56,57). Furthermore, studies have also highlighted the stimulatory role of MT-1s during angiogenesis and collaterogenesis under stress conditions such as hypoxia (58). In hypoxic conditions, MT-1s expression increased, and it was able to stabilize the expression of hypoxia-inducible factor-1a (HIF-1a) (59). MT-1s can also induce proangiogenic factors such as fibroblast growth factor (FGF), transforming growth factor b (TGF-b), and vascular endothelial growth factor A (VEGF-A) (60). In spite of these cytoprotective activities, the over-expression of MT-1s is linked with an unfavorable prognosis in several cancers including NSCLC (61), intrahepatic cholangiocarcinoma (ICC) (62), renal cancer (63), ovarian cancer (64), and melanoma (65). However, the prognostic value of MT-1s seems to be limited in some tumors, such as thyroid or laryngeal cancer (66). Some isoforms of MT-1s have been shown to be cancer suppressors. The MT-1G isoform has been found to be a potent oncosuppressor in papillary thyroid (39) and hepatocellular cancers (67), and MT-1F may play a suppressive role in colon cancer (68).

GO annotation and enrichment analysis of twelve differentially expressed transcripts including MT1G, MT1H, MT1F, TPO, uc021yst, SLC34A2, SLC26A4, IYD, n340297, n407932, n342369, and n336033 was implemented using the local FunNet software (28), with results shown in . With respect to molecular functions, 96 GO terms were assigned of which the three most significantly enriched terms were protein binding, RNA binding, and structural constituent of ribosome. With respect to cellular components, 70 GO terms were assigned of which cytoplasm, membrane, and cytosol were the most significantly enriched. With respect to biological processes, 237 GO terms were assigned of which gene expression and cellular protein metabolic process were the top two significantly over-represented terms. A KEGG analysis revealed that protein processing in the endoplasmic reticulum was the most enriched pathway among differentially expressed transcripts, with eight of these transcripts participating in this pathway. In addition, oxidative phosphorylation, Huntington’s disease, and Alzheimer’s disease was also significant enriched in this pathway analysis.

In conclusion, our study characterized the significantly altered expression profiles of lncRNAs and mRNAs in PTC. Abnormal regulation of these altered lncRNAs and mRNA may play vital role in the development, progression, and metastasis of this disease. Further research will be needed to explore whether these alterations are linked with genotype, familial cancer risk, inflammatory states, or a combination of these cancer risk factors. In addition, further study of the molecular mechanisms in which these mRNAs and lncRNAs participate may yield further insights into PTC pathogenesis.