RNA-Seq transcriptome analysis reveals Maackia amurensis leukoagglutinin has antitumor activity in human anaplastic thyroid cancer cells.

BACKGROUND: Lectins are carbohydrate-binding molecules that can bind specifically to the sugar residues of glycoconjugates and are found in almost all organisms. Plant lectins subjected to many studies reported exhibiting anti-cancer activity. This study aimed to investigate the possible molecular mechanisms of Maackia amurensis leukoagglutinin II (MAL-II) treated ATCCs. METHODS AND
RESULTS: We tested the effects of MAL-II, which is isolated from Amur seeds, on cancerous features of 8505C human anaplastic thyroid cancer cells (ATCCs) on a large scale using RNA-Seq. Transcriptome analysis was performed using Illumina next-generation sequencing technology by using cDNA libraries obtained from total RNA isolates of ATCCs treated with 0.25 µM MAL-II for 24 h. Gene ontology and pathway enrichment analysis were performed for the systematic analysis of gene functions. Moreover, we validated RNA-Seq findings using qPCR. Our results showed that many cancer-related genes such as TENM4, STIM2, SYT12, PIEZO2, ABCG1, SPNS2, ARRB1, and IRX5 were downregulated and many anticancer genes such as HSPA6, G0S2, TNFAIP3, GEM, GADD45G, RND1, SERPINB2, and IL24 were upregulated. Also, pathway enrichment analysis showed that differentially expressed genes were found to be associated with Ras, p53, and apoptosis signaling pathways, which are some important signal transduction pathways in development, proliferation, stem cell control, and carcinogenesis.
CONCLUSION: Collectively, our results show that MAL-II treatment reveals significant antitumor activity by changing the expression of many cancer-related genes and implies that MAL-II treatment might be a potential candidate molecule to inhibit the malignancy of human anaplastic thyroid cancer.

Introduction

Thyroid cancer is the most common endocrine system cancer worldwide, accounting for approximately 2.1% of all malignant cancers [1]. The incidence is higher in developed countries and 5.12% of women were affected worldwide in 2018. The incidence in women is 3 times higher than in men. The incidence of thyroid cancer has increased 3 fold in the last thirty years, and the incidence-based mortality rate has increased by 1.1% per year [2]. It is histologically classified into four groups: papillary, follicular, medullary, and undifferentiated or anaplastic carcinoma [3]. Papillary, follicular, and anaplastic thyroid carcinomas (ATC) are derived from follicular cells, while medullary carcinomas are derived from para-follicular cells [4]. These different types of thyroid cancer show different features in terms of tumor differentiation, metastatic potential, and aggressiveness. Anaplastic thyroid cancer cells (ATCCs) lose their thyroid-specific gene expression patterns, such as the ability to take up iodine, and production of thyroglobulin during the carcinogenesis process. In particular, the sodium-iodide symporter (NIS) gene, which mediates the uptake of iodide in thyrocytes, is not expressed in ATCCs. Thus, the methods such as radioiodine or TSH suppressor therapies used to treat patients with differentiated thyroid cancer are not often effective in patients with undifferentiated thyroid cancer [5]. This type of tumor also fails to respond to alternative treatments such as external radiation and systemic cancer chemotherapies [6]. Therefore, there is a need for new strategies for treating anaplastic thyroid cancer.

Some approaches suppress the tumorigenic and malignant properties of thyroid cancer using lectins. It has been reported that Maackia amurensis leukoagglutinin (MAL-II), a plant lectin obtained from the Amur seeds, has an anti-tumoral effect on ATCCs and is promising for treating ATC [7]. Lectins are a family of proteins that can bind specifically to carbohydrates and are found in various organisms, from viruses to bacteria and from plants to animals [8]. They have been subjected to research in many cancers because of their ability to stimulate or suppress various cellular processes in different ways [7, 9]. Lectins perform anti-tumoral effects by binding to the cell surface carbohydrate residues, which are being altered due to the oncogenic stimuli and so-called abnormal glycosylation in cancer cells [9]. Many investigators have also pointed out that changes in cell surface glycans play an important role in promoting the malignant characteristics of the tumors [10]. Currently, changes in some glycan structures have come forward as tumor markers as well [11]. Studies have shown that glycocalyx is rich in α-2,3 sialic acid residues in many types of cancer, including human thyroid cancer cells [12]. Ochoa-Alvarez et al. [13] stated that MAL-II treatment decreased the invasiveness of melanoma cells by binding increasing α2,3 sialic acid motifs in podoplanin, which is the extracellular domain of the mucin receptor giving rise to an increase invasive trait of melanoma cells. Also, MAL-II treatment reduced cell motility in oral squamous carcinoma cells expressing podoplanin [14]. In another study, it was observed that only α1 unit of α5β1 integrin in melanoma cells carries α-2,3 sialic acid motifs and MAL-II treatment resulted in a large decrease in the motility of the cells [15]. Our previous study revealed that MAL-II treatment significantly reduced the tumorigenic and malignant characteristics of ATCCs. MAL-II treatment significantly reduced tumor-endothelial cell interaction, which is an important step in the metastatic process. Besides, non-toxic doses of MAL-II to normal thyroid epithelial cells showed toxic effects in the 0.1– 1 µM dose range in the 8505C human anaplastic cancer line [7]. However, there is no information about its molecular mechanism of action.

In this context, transcriptome analysis can be one of the potential approaches to comprehending the underlying mechanisms of MAL-II treatment on a large scale. Transcriptome analysis using RNA sequencing is a technology that provides the opportunity for the simultaneous analysis of multiple genes/targets and the identification of mechanisms of action in drug treatment experiments [16]. Detection and grouping of transcripts (mRNAs, non-encoded RNAs) of species, understanding of post-translational modifications, elucidation of splicing mechanisms, and comparison of transcript-expression levels under different conditions [17]. Transcriptome studies are important for understanding the genomic characterization and transcriptomic outcomes of various cancers and their contribution to the implementation of sustainable and effective treatment methods for each patient because of the identification of biomarkers for the diagnosis and course of disease [18]. Although MAL-II has been shown to inhibit the proliferative and invasive properties of ATCCs, the mechanism of action has not been investigated. In this study, it was analyzed the differences or changes in whole genome-wide transcriptome profile in ATCCs under MAL-II treatment using RNA-Seq analysis.

Materials and methods

Cell culture and treatment

The human ATC cell line 8505C was obtained from the European Collection of Authenticated Cell Cultures (ECACC). The cells were cultured in Eagle’s Minimum Essential Medium (EMEM/Gibco) containing 10% fetal bovine serum (FBS/Biowest) and antibiotic (100 U/ml penicillin and 100 µg/ml streptomycin/Gibco) and incubated in the conditions of 37 °C and 5% CO2 at 1 atmospheric pressure. The stock lectin solution was prepared at a concentration of 2 mg/ml in HEPES buffer using lyophilized MAL-II supplied by Vector Lab (L-1260-2). Two different lectin doses were prepared from the stock solution at the concentrations of 0.25 and 0.5 µM in a cell culture medium and treated to ATCCs.

Cell viability assay

Thiazolyl blue tetrazolium bromide (MTT/Biomatik A3338) assay was used to detect the effect of MAL-II treatment on the viability of 8505C cells. For this purpose, 5 × 103 cells per well were seeded in 96 well plates. On the following day,the culture medium was removed from the wells and fresh medium was added to the control group, and medium containing 0.25 and 0.5 µM MAL-II to the experimental group. After 24 and 48 h of incubation with MAL-II, 10 µl of 5 g/L MTT was added to each well and incubated under the condition of 37 °C, 5% CO2 for 4 h. The medium with MTT was withdrawn from the wells and then 100 µl dimethyl sulfoxide (DMSO) was added to each well to dissolve formazan crystals formed by the reduction of MTT. The plate was shaken on an orbital shaker at 150 rpm for 5 min. Optical density was measured at 570 nm test wavelength and 630 nm reference wavelengths.

RNA extraction

8505C cells were inoculated in 6 well plates (3 × 105 cells per well). Following day, 0.25 µM MAL-II was added to each well of the experimental group, whereas a fresh culture medium was added to each well of the control group and incubated for 24 h. The instructions Ambion PureLink RNA isolation kit (Life Technologies) was used to extract RNA and the content, purity, and integrity of RNA were detected using Nanodrop ND-2000 (Thermo Scientific) and Agilent-2100 Bioanalyzer. The ratio of OD260/OD280 was between 1.8 and 2.0 with high purity. RNA integrity number (RIN) was at least above 7 for all samples on a scale from 1 to 10.

RNA-Seq and data analysis

Transcriptome analysis of 8505C cells was performed on a Genome Sequencer Illumina HiSeq next-generation sequencing system in sequence mode HiSeq 4000 50 bp SR. The RNA-Seq reads were aligned to the reference genome (hg19 genome assembly) using Bowtie generating transcriptome alignments. Cufflinks identified and quantified the transcripts from the pre-processed RNA-Seq alignment assembly. After this, Cuffmerge merged the identified transcript pieces into full-length transcripts and annotated the transcripts based on the given annotations. Finally, merged transcripts from two or more samples/conditions were compared using Cuffdiff to determine the differential expression levels at transcript and gene levels, including a measure of significance between samples. Operating on the RNA-Seq alignments and Cufflinks processing, Cuffdiff tracked the mapped reads and determined the fragment per kilobase per million mapped reads (FPKM) for each transcript in the samples. Primary transcripts and gene FPKMs were then computed by adding up the FPKMs of each primary transcript group or gene group.

Validation of RNA-Seq data by qPCR

Total RNA isolated from the cells was converted into cDNA using High-Capacity cDNA Kit. 200 ng of cDNA was used as a template in the qPCR reaction. The following primers were used for the amplification of the desired mRNA transcript. Leucine-rich G protein-bound receptor 5 (LGR5: forward 5′-TCCGGTCGTTGGAAGAACTCTCG-3′, reverse 5′-TCAGTCGCCACT GCACATCCAGC-3′) and prostaglandin F receptor (PTGFR: forward 5′-CTCCGTCTTCTG CTCCTCAGAGA-3′, reverse 5′-GGAATTGTTCATGGACATTGTTGTG-3′), Glipican 5 (GPC5: forward 5′-TCCGGTCGTTGGAAGAACTCTCG-3′, reverse 5′-tcagtc GCCACTGC ACATCCAGC-3′) mRNA transcripts. The qPCR was performed in Applied Biosystems 7500 fast real-time PCR system by using SYBR Green. The cycling parameters were as follows: reverse transcription at 48 °C for 30 min, AmpliTaq activation at 95 °C for 10 min, denaturation at 95 °C for 15 s, and annealing/extension at 60 °C for 1 min (cycle repeated 40 times) according to GeneMark SYBR kit protocol. The relative changes in gene expression were analyzed using 2−ΔΔCt method. The amount of fold change in the target mRNA was normalized to an endogenous control (GAPDH: forward: 5′-GATGACATCAAGAAGGTGGTGA-3′ reverse: 5′-TTCGTTGTCATACCAGG AAATG-3′).

Functional and pathway enrichment analyzes

Gene Ontology enrichment analysis was performed using ShinyGO (v.0.60) to predict the enrichment degree and the potential functions of differentially expressed genes of MAL-II treated 8505C cells in biological processes (BP), cellular components (CC), and molecular functions. Additionally, Panther (16.0) pathway enrichment analysis was used for systematic analysis of gene functions, revealing a statistically significant difference. Differentially expressed genes with p < 0.05; log-fold change cutoff = 0 were used for the enrichment analyses.

Statistical analysis

To determine statistical significance between MAL-II treated and control groups, one-way ANOVA, and Tukey multiple comparison post hock tests were performed. T-test was used for statistical evaluation of qPCR gene validation data. Results are presented as the mean ± SD of the mean. GraphPad Prism (v.5.0.1) software was used for statistical calculations.

Results

The cellular viability of MAL-II treated 8505C anaplastic thyroid cancer cells

MAL-II treatment to 8505C cells significantly reduced cell viability at both doses of 0.25 and 0.5 µM. MAL-II treatment of 0.25 µM decreased approximately 35% and 60% in cell viability after 24 and 48 h respectively (p ≤ 0.001). Also, 0.5 µM MAL-II treatment reduced approximately 30% and 60% in cell viability during the same time span (p ≤ 0.001) (Fig. 1A).

Fig. 1

Effects of MAL-II on cell viability and transcriptomic activity of 8505C cells. A MTT cell viability results of 8505C cells under MAL-II treatment at the concentrations of 0.25 and 0.5 µM for 24 and 48 h. Percent viability is given as the mean ± SD of three independent experiments. ***(p ≤ 0.001). B Scatter plot for MAL-II untreated versus treated cells. C Heatmap plot of 100 most abundant genes overlapped genes between MAL-II untreated and treated cells. D Volcano plot of DEGs for MAL-II untreated versus treated cells. Black points represent insignificant DEGs, red points represent significant DEGs of the MAL-II treated cells. E Profiling of upregulated and downregulated significant DEGs (p ≤ 0.05) in MAL-II treated cells

RNA-Seq analysis in the MAL-II treated cells

RNA-seq was performed in both 24 h, 0.25 µM MAL-II treated and untreated 8505C cells. The number of reads mapped to the reference genome for each sample in the experiment ranged from 34 million to 41 million pairs and the percentage of genomic alignment consisted of 99.4% clean reads. Scatter plots in Fig. 1B indicate the general similarities and specific outliers between MAL-II treated and untreated conditions in the RNA-Seq experiment. To inspect the quality of RNA-Seq data, the 100 most abundant genes were taken from the samples, and heatmaps were generated to observe the relation between treated and untreated samples (Fig. 1C). The distributions of the fold changes and p-values of genes performed on the RNA-Seq data are shown in Fig. 1D as volcano plots. The gene expression profiling showed that there were numerous significantly differently expressed genes (p ≤ 0.05) in MAL-II treated ATCC (Fig. 1E).

Identification of differentially expressed genes in the MAL-II treated cells

We identified 8037 significantly differentially expressed genes. However, to facilitate the evaluation of these genes, the threshold level was set as FDR and p-value below 0.01 and as ± 3-fold change limit. Considering these threshold levels, we found that 27 genes were upregulated and 29 genes downregulated among significantly differentially expressed 56 genes in MAL-treated 8505C cells compared to non-MAL treated cells (Tables 1 and 2).

Table 1

Top downregulated and differentially expressed genes (p < 0.01, FDR < 0.01 and log2ratio ≥ 3)

Accession no.	Gene symbol	Gene name	Log2 ratio	p value	FDR
NR_026779	LINC00341	long intergenic nc RNA 341	− 4.81	5E−05	1.97E−04
NM_001098816	TENM4	teneurin transmembrane protein 4	− 4.24	5E−05	1.97E−04
NM_001177880	SYT12	synaptotagmin 12	− 3.90	5E−05	1.97E−04
Unknown	LINC00939	long intergenic nc RNA 939	− 3.88	5E−05	1.97E−04
NM_022068	PIEZO2	piezo type mechanosensitive icc 2	− 3.73	5E−05	1.97E−04
NM_001128159	VPS53	VPS53, GARP complex subunit	− 3.71	1.6E−03	4.73E−03
NM_004041	ARRB1	arrestin beta 1	− 3.71	5E−05	1.97E−04
NM_198427	BCAN	brevican	− 3.51	5E−05	1.97E−04
NM_003179	SYP	synaptophysin	− 3.48	5E−05	1.97E−04
NM_005853	IRX5	iroquois homeobox 5	− 3.45	4.5E−04	1.50E−03
NR_036540	LINC00622	long intergenic nc RNA 622	− 3.41	5.5E−04	1.81E−03
Unknown	FAM229A	Family with sequence similarity 229 M.A	− 3.39	5E−05	1.97E−04
NM_001146041	KRTAP4− 9	keratin associated protein 4–9	− 3.34	1.55E−03	4.60E−03
NM_002587	PCDH1	protocadherin 1	− 3.33	5E−05	1.97E−04
NM_001093772	KIT	proto-oncogene receptor tyrosine kinase	− 3.33	5E−05	1.97E−04
NM_001410	MEGF8	multiple EGF like domains 8	− 3.33	5E−05	1.97E−04
NM_001124758	SPNS2	sphingolipid transporter 2	− 3.30	3.5E−04	1.19E−03
NM_014339	IL17RA	interleukin 17 receptor A	− 3.24	5E−05	1.97E−04
NM_199461	NANOS1	nanos C2HC-type zinc finger 1	− 3.24	5E−05	1.97E−04
NR_036476	TMEM198B	transmembrane protein 198B	− 3.21	7E−04	2.25E−03
NM_017525	CDC42BPG	CDC42 binding protein kinase gamma	− 3.20	3.6E−03	9.70E−03
NM_207627	ABCG1	ATP binding cassette subfamily G M.1	− 3.18	5E−05	1.97E−04
NM_178544	ZNF546	zinc finger protein 546	− 3.17	2.2E−03	6.27E−03
NM_020860	STIM2	stromal interaction molecule 2	− 3.12	5E−05	1.97E−04
NM_001134402	CYB561D1	cytochrome b561 family member D1	− 3.11	5E−05	1.97E−04
Unknown	LINC00896	long intergenic nc RNA 896	− 3.11	1.95E−03	5.64E−03
NR_024490	AS1	GABPB1 antisense RNA 1	− 3.10	5E−05	1.97E−04
NM_002048	GAS1	growth arrest specific 1	− 3.09	5E−05	1.97E−04
NR_024396	LINC00638	long intergenic nc RNA 638	− 3.02	5E−05	1.97E−04

Table 2

Top upregulated and differentially expressed genes (p < 0.01, FDR < 0.01 and log2 ratio ≥ 3)

Accession no.	Gene symbol	Gene name	Log2Ratio	Pvalue	FDR
NM_002155	HSPA6	Heat shock protein family A, member 6	6.86	5E−05	1.96−E04
NM_001130046	CCL20	C–C motif chemokine ligand 20	6.81	5E−05	1.96−E04
NM_002090	CXCL3	C-X-C motif chemokine ligand 3	6.30	5E−05	1.96−E04
NM_000594	TNF	Tumor necrosis factor	6.11	5E−05	1.96−E04
NM_015714	G0S2	G0/G1 switch 2	6.01	5E−05	1.96−E04
NM_002089	CXCL2	C-X-C motif chemokine ligand 2	5.91	5E−05	1.96−E04
NM_002517	NPAS1	Neuronal PAS domain protein 1	4.95	5E−05	1.96−E04
NM_006705	GADD45G	Growth arrest and DNA damage inducible ɣ	4.53	5E−05	1.96−E04
NM_004049	BCL2A1	BCL2 related protein A1	4.33	3.5−E03	9.58−E03
NM_018724	IL20	Interleukin 20	4.32	2.2−E03	6.27−E03
NM_014470	RND1	Rho family GTPase 1	4.30	5E−05	1.96−E04
NM_000584	IL8	Interleukin 8	4.16	5E−05	1.96−E04
NM_002575	SERPINB2	Serpin family B member 2	4.08	5E−05	1.96−E04
NM_021009	UBC	Ubiquitin C	4.03	5E−05	1.96−E04
Unknown	ACKR3	Atypical chemokine receptor 3	3.99	5E−05	1.96−E04
NM_002982	CCL2	C-C motif chemokine ligand 2	3.83	5E−05	1.96−E04
NM_006290	TNFAIP3	TNF alpha induced protein 3	3.70	5E−05	1.96−E04
NM_001185157	IL24	Interleukin 24	3.70	5E−05	1.96−E04
NM_020529	NFKBIA	NFKB inhibitor alpha	3.70	5E−05	1.96−E04
NM_000575	IL1A	Interleukin 1 alpha	3.65	5E−05	1.96−E04
NM_000636	SOD2	Superoxide dismutase 2, mitochondrial	3.47	5E−05	1.96−E04
NM_000576	IL1B	interleukin 1 beta	3.46	5E−05	1.96−E04
NM_005346	HSPA1B	Heat shock protein family A, member 1B	3.23	5E−05	1.96−E04
NM_005261	GEM	GTP binding protein	3.22	5E−05	1.96−E04
NM_004428	EFNA1	Ephrin A1	3.14	5E−05	1.96−E04
NM_001168238	ABL2	ABL proto-oncogene 2	3.06	5E−05	1.96−E04
NM_004466	GPC5	Glipican 5	Infinity	5E−05	1.96−E04

To validate the RNA-Seq results, three genes were selected from among the most abundant genes (GPC5, LGR5, PTGFR) and performed qPCR using extracted RNA from different experimental setups. The correlation was detected between relative expression levels of RNA-seq and q-PCR data (R = 0.96, p = 0.16) (Fig. 2A). The results of qPCR showed that the mRNA level of GPC5 increased significantly (p ≤ 0.01) but LGR5 (p ≤ 0.001) and PTGFR (p ≤ 0.01) significantly decreased compared to the control (Fig. 2B).

Fig. 2

Validation of RNA-Seq data. A Correlation of gene expression ratios obtained from RNA-Seq and qPCR assays. The most abundant three genes (GPC5, LGR5, and PTGFR) obtained from RNA-Seq analysis were validated by qPCR. B The relative expression levels of GPC5, LGR5, and PTGFR were determined by qPCR. Fold change is given as mean ± SD, GPC5 Control: 1.38 ± 0.35 vs. MAL-II: 34.39 ± 10.30 **(p ≤ 0.01). LGR5 Control: 1.00 ± 0.16 vs. MAL-II: 0.08 ± 0.01 ***(p ≤ 0.001). PTGFR Control: 1.00 ± 0.13 vs. MAL-II: 0.31 ± 0.04 **(p ≤ 0.01). The statistical significance between control vs. MAL-II treatment was assessed by Student’s t-test

Functional and pathway enrichment analysis of the MAL-II treated cells

To show the roles of differentially expressed genes, significantly enriched functional and pathway analysis would help to understand the underlying mechanism of the effects occurring in MAL-II treatment in the cell. Functional enrichment analysis of the significantly differentially expressed genes mostly took part in several BP including organelle organization, macromolecule modification, cellular protein modification process, and protein modification process. When differentially expressed genes were taken into account in CC, most of the genes were enriched in nuclear compartments such as the nuclear lumen, nuclear part, nucleoplasm, nucleolus, and ribonucleoprotein complex. Molecular function analysis, however, showed that differentially expressed genes were mostly related to RNA binding, nucleic acid binding, enzyme binding, transcription factor binding, transferase activity, and nucleotide-binding (Supplementary Table 1). Besides, PANTHER pathway analysis (Overrepresentation test version 16.0 released 2020-12-01) was conducted to determine the biological roles of the significantly differentially expressed genes. The pathway analysis showed that there were 10 significantly enriched pathways listed in Table 3 with FDR of < 0.05 (Binominal test, p < 0.05). The most conventional pathways were the cell cycle, Ras pathway, p53, apoptosis, and integrin signaling pathways (Table 3).

Table 3

PANTHER pathway analysis of differentially expressed genes (p ≤ 0.05)

ID	PANTHER pathways		Fold enrichment	Raw p value	FDR
P00013	Cell cycle		2,32	1.57E−03	2.91E−02
P04393	Ras Pathway		1,89	5.77E−05	4.82E−03
P00060	Ubiquitin proteasome pathway		1,82	5.68E−04	1.36E−02
P00046	Oxidative stress response		1,77	1.70E−03	2.85E−02
P00059	p53 pathway		1,72	2.30E−04	6.41E−03
P00006	Apoptosis signaling pathway		1,68	6.31E−05	2.64E−03
P06959	CCKR signaling map		1,55	6.20E−05	3.45E−03
P00047	PDGF signaling pathway		1,55	2.06E−04	6.89E−03
P00049	Parkinson disease		1,54	2.36E−03	3.58E−02
P00034	Integrin signalling pathway		1,40	1.41E−03	2.95E−02

Discussion

The altered or abnormal glycosylation profile of cancer cells is an important category for potential biomarker, drug, and vaccine development. Also, accumulated evidence has shown that the treatment of cancer cells with the exogenous lectins, specific for abnormal glycan motifs of cancer cells, can be a promising strategy to overcome malignant cancers. MAL-II-specific for α-2,3, which is one of the most abundant abnormal glycosylation motif in anaplastic thyroid cancer, suppresses proliferative-invasive characteristics of anaplastic thyroid cancer and significantly reduces tumor cell-endothelial interaction [7]. The question of this study was how the transcriptome of anaplastic thyroid cancer cells changes in response to MAL-II treatment. It also aimed to find out the mechanism of action of MAL-II.

In this study, MAL-II treatment to the anaplastic thyroid cancer cell line 8505C at concentrations of 0.25 and 0.5 µM for 24 and 48 h was revealed to significantly reduce cell viability. This finding we obtained confirmed the results previously reported by Kaptan et al. [7]. Transcriptome analysis showed that TENM4, STIM2, SYT12, PIEZO2, ABCG1, SPNS2, ARRB1, and IRX5 genes, which are previously reported in cancers such as breast, lung, and colorectal cancers [19-26], were downregulated in the cells treated with MAL-II. However, only SYT12 and ABCG1 are also associated with thyroid cancer [27, 28]. The SYT12 gene is usually over-expression in PTC, it was assumed that SYT12 plays a significant role in tumorigenesis and progression of thyroid cancer. Knocked down the expression of SYT12 observably suppressed the proliferation ability of papillary thyroid cancer cell lines [27]. ABCG1 is a member of the ATP-binding cassette family of proteins and its increased mRNA expression has been shown in anaplastic thyroid cancer and is associated with chemoresistance and malignancy [28]. Thus, our findings indicated that MAL-II treatment downregulated the genes, which affected cancerous phenotype. This can mean that MAL-II is a potential therapeutic agent for anaplastic thyroid cancer. Also, future studies on these genes can provide further insight to understand the underlying molecular mechanism of the malignant phenotype of anaplastic thyroid cancer to determine new diagnostic and therapeutic targets.

When we had a look at upregulated gene list, they were not associated with thyroid cancer. However, many of these have previously been reported in cancer studies. Especially, HSPA6, G0S2, TNFAIP3, GEM, GADD45G, RND1, SERPINB2, and IL24 have come forward with their repressive roles in malignant features of some cancers [29-35]. However, some genes like ABL2, EFNA1, CCL20, CXCL3, BCL2A1, IL20, and IL8, support cancer progression and they have already been highlighted as indicators of poor prognosis or suggested as potential therapeutic targets [36-40]. This finding suggests that upregulated anti-cancerous genes override the effects of upregulated cancerous genes when considering MAL-II previously demonstrated effects [7] on ATCCs. Additionally, in a study by Hamilton et al. [41], mRNA expression of UBC, which is also found in our upregulated gene list, increased in oral squamous cell carcinoma cells treated with Maackia amurensis seed lectin (MASL). The UBC gene encodes polyubiquitin C protein that maintains cellular ubiquitin levels under stress conditions [41]. Thus, upregulation of UBC and some other cancer-related genes resulting in MAL-II exposure may occur in response to stress.

The pathway analysis displayed that differentially expressed genes significantly enriched 10 pathways. One of the most enriched pathways was Ras, also known as MAPK/Erk pathway. The MAPK pathway is the focal point for therapeutic approaches in thyroid cancer due to the high rate of mutations in the MAPK pathway associated with survival and proliferation [42]. It has been suggested that the simultaneous targeting of MAPK signaling pathways using clinically approved inhibitors can be an effective therapeutic strategy for aggressive and recurrent thyroid tumors [43]. In this respect, the decreased survival rate with MAL-II treatment may be related to the MAPK pathway. Deregulation of the p53 pathway is an important step for anaplastic thyroid cancer progression [44]. Therefore, MAL-II can be considered as a bioactive molecule that affects both pathways. Plant lectins affect both apoptosis and autophagy by modulating representative signaling pathways in the Bcl-2 family, caspase family, p53, PI3K/Akt, ERK, BNIP3, Ras-Raf, and ATG families in cancer [45]. Comparative studies using various lectins as cytotoxic and anti-tumoral agents revealed different results depending on the source of the lectin and the type of cancer. Concanavalin A (Con A) is an α-mannose-specific lectin that downregulates different signaling pathways mediated through NF-kβ, Erk, and JNK, and stimulates p53 for the induction of apoptosis in different cancer cells [46]. Polygonum cyrtonema lectin is also a mannose-specific lectin that inhibits the activation of Ras-Raf and PI3K-Akt signaling pathways and induces apoptosis in L929 murine fibrosarcoma [47]. Moreover, it has been shown that PCL promotes apoptosis through the activation of MAPK and NF-kβ pathways via ROS generation in A549 lung cancer cells [48]. Additionally, galactoside-specific Mistletoe lectin exhibits anticancer activity via the regulation of Erk, p38, Akt, and SAPK/JNK pathways [49]. Polygonatum odoratum lectin is another mannose-specific lectin that affects the EGFR-mediated Ras-Raf-MEK-ERK signaling pathway and stimulates MCF-7 cell apoptosis and autophagy [50]. Hence, our findings regarding the most enriched pathways might indicate mechanisms of MAL-II-dependent cellular toxicity.

In summary, the large amount of data obtained by RNAseq analysis provided a broad framework for evaluating the effect of MAL-II Maackia amurensis leukoagglutinin II (MAL-II) treatment on anaplastic thyroid cancer cells. We detected many genes that were up and downregulated genes when cells were treated with MAL-II. Many of the downregulated genes were associated with the functions supporting malignancy and poor prognosis in most widely common cancers. The functional and pathway analysis of these data leads us to the conclusion that some biological processes like cell cycle, proliferation, and some molecular functions such as transcription factor binding and RNA binding in anaplastic thyroid cancer cells were affected by MAL-II treatment. Ras, p53, and the apoptosis pathway were affected as well. These results indicate some of the cancer-related genes, which are previously known, might also play role in cancerous features of anaplastic thyroid cancer. Thus, further studies are needed to confirm the validity of these genes for adoption as therapeutic targets and as functional players in anaplastic thyroid cancer. Our results suggest that MAL-II might be a potential bioactive agent for the treatment of anaplastic thyroid cancer, which is resistant to radioiodine therapy.

Below is the link to the electronic supplementary material.

Supplementary material 1 (DOCX 15.2 kb)