Immune-related key gene CLDN10 correlates with lymph node metastasis but predicts favorable prognosis in papillary thyroid carcinoma

Relationship of immune and stromal scores with clinical features

The immune and stromal scores and the clinical information of the TCGA cohort were obtained. The characteristics of the PTC cohort are presented in Table 1. The immune and stromal scores were respectively divided into high and low groups on the basis of their median levels. As shown in Table 2, high immune scores showed a significant association with advanced clinical stages, high T stages, and LNM (P<0.05). Although the comparisons of stromal scores were not significant, a borderline association was shown in the LNM comparison (P=0.07). The data indicated that both immune and stromal scores might have a correlation with LNM in PTC.

Survival curves based on the TCGA cohort were constructed to evaluate the prognostic values of the immune and stromal scores. The results of log-rank tests failed to show that either immune or stromal scores influenced the prognosis of PTC (P>0.05; Figure 1A).

Identification of differentially expressed genes (DEGs) related to immune scores and stromal scores in PTC

Immune and stromal scores were used to stratify patients with PTC into high and low groups, respectively, as mentioned above. The gene expression profiles were compared between the high and low groups. Heatmaps in Figure 1B show the results generated from the comparisons. The comparison based on immune scores (high vs. low) revealed 1324 upregulated and 112 downregulated genes. Similarly, the comparison based on stromal scores (high vs. low) generated 1237 upregulated and 5 downregulated genes.

Moreover, Venn analysis was performed to narrow the scope of the target genes that are associated with immune and stromal cells. The data showed 1003 upregulated and 5 downregulated DEGs in the intersection (Figure 1C). Thus, these genes were chosen for further analysis.

Functional annotation of DEGs in the intersection

Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses were utilized to evaluate the possible functions of the screened DEGs in the intersection.

As shown in Figure 1D, the top five GO terms were immune response, defense response, response to biotic stimulus, response to pest, pathogen or parasite, and organismal physiological process. The top five pathways were cytokine–cytokine receptor interaction, complement and coagulation cascades, Toll-like receptor signaling pathway, ECM–receptor interaction, and focal adhesion.

The results of GO analysis showed that these DEGs might be involved in various cell biological processes, including responses to stimulus and immune response. The data were consistent with our results that the screened DEGs may be correlated with immune responses. Pathway enrichment analysis showed that these DEGs may be enriched in pathways related to cancer development, suggesting that these genes function in immune response and cancer progression.

Prognostic values of individual DEGs in PTC

Univariate cox regression analyses of the TCGA cohort were conducted to explore the prognostic values of the screened DEGs in the intersection. Of the 1008 DEGs that were assessed, 87 predicted favorable or unfavorable prognosis for PTC (P<0.05; Table 3).

Construction of co-expression network and identification of hub genes associated with LNM

A total of 1008 DEGs in the intersection were included in the co-expression network analysis. All samples were selected because no outliers were observed among the samples.

The soft threshold was determined by scale independence and mean connectivity analysis of modules with different power values. As shown in Figure 2A, when the power value was set to 9, the scale independence value achieved 0.9 and lower mean connectivity. Figure 2B shows the cluster dendrogram among the modules. Four modules containing 194 genes in blue, 126 in brown, 15 in grey, and 671 in turquoise were generated with different colors.

The module that had a relationship with LNM was evaluated. As shown in Figure 2C, the blue module may have the closest association (r=0.38, P=4.6e-08). Thus, the genes in the blue module were selected for further analysis.

Key gene selection

The following steps were carried out to screen the immune-related and LNM-related key genes in PTC. First, 87 genes that had prognostic values for PTC, as shown in Table 3, were selected (immune-related). Second, the genes in the blue module (LNM-related), with a GS value more than 0.15, were selected (20 genes). Third, Venn analysis was conducted to obtain the intersection of these genes. As a result, CLDN10 was screened as the key gene. Hence, the key gene might be immune- and LNM-related (Figure 2D).

Assessment of CLDN10 expression in a PTC cohort from the TCGA and gene expression omnibus (GEO) databases

The expression data of CLDN10 in a PTC cohort were extracted from the TCGA database. The expression levels of CLDN10 were classified as high and low on the basis of the median level. Table 4 lists the main results of the association of CLDN10 expression and the clinical features.

As shown in this table, high expression may have a significant correlation with low ages, high T stages, and LNM (P<0.05).

Relevant data from a cohort in a dataset (GSE35570) were downloaded from the gene expression omnibus (GEO) database to evaluate the diagnostic value of CLDN10 expression in PTC [12]. This cohort included 65 PTC cases and 51 normal controls. As shown in Figure 3A-a, the area under the ROC curve (AUC) of the receiver-operating characteristic (ROC) curves achieved 0.954, with a sensitivity of 0.9077 and a specificity of 0.9804, suggesting that CLDN10 can be used as a potential biomarker distinguishing PTC cases from normal controls.

Surprisingly, survival analysis showed that patients with low CLDN10 expression may have a shorter overall survival time than those with high CLDN10 expression (Figure 3A-b). In other words, CLDN10 might act as a preventive factor for PTC.

Protein expression of CLDN10 in a PTC cohort on the basis of a tissue microarray

The association of CLDN10 with LNM was further validated by measuring its protein expression on a tissue chip through an immunohistochemical (IHC) assay.

Specific staining was observed mainly in the membrane of cancer and normal cells, and weak staining was observed in cytoplasm (Figure 3B-ab). In addition, the expression scores of CLDN10 were significantly higher in PTC tissues than in para-carcinoma tissues (P<0.05), as shown in Figure 3B-c.

In this cohort, no cases presented distant metastasis, and the survival information was unavailable. Only the confounding factors, such as clinical stages, lymph node metastasis, and age, can be addressed.

Interestingly, the data indicated that the high expression score of CLDN10 was significantly associated with LNM (P<0.05). No significant associations were observed in the comparisons regarding age, T stage, and clinical stage (P>0.05) (Figure 3C).

Association of CLDN10 expression with immune cell infiltration levels

The TIMER algorithm [13] was used to determine the possible association between CLDN10 expression and immune cell infiltration.

On the basis of the PTC cohort from TCGA, a weak negative correlation was observed between CLDN10 expression and tumor purity (r=-0.133, P<0.05). Conversely, CLDN10 expression had significantly positive correlations with infiltrating levels of B cells (r=0.348, P<0.05), CD4+T cells (r=0.351, P<0.05), neutrophils (r=0.509, P<0.05), and dendritic cells (r=0.52, P<0.05) (Figure 4A). Moreover, weak positive correlations of CLDN10 expression with CD8+T cells (r=0.105, P<0.05) and macrophages (r=0.181, P<0.05), respectively, were observed. These findings indicate that CLDN10 is closely associated with immune cell infiltration in PTC.

Multivariate COX regression analyses were performed using the TIMER tool to learn the potential effects of these immune cells on PTC survival. As shown in Table 5, tumor purity was correlated with poor survival, confirming that CLDN10 expression was negatively correlated with tumor purity. Notably, CD8+T cells (HR=0.000, 95%CI= 0.000–0.032) and macrophages (HR=0.000, 95%CI= 0.000–0.152) might be independent favorable prognostic indicators that lower the cancer risk and prolong the overall survival time of patients.

Association of CLDN10 expression with immune cell markers

Markers of immune cells were considered to investigate further the association of CLDN10 expression with CD8+T cells and macrophages. As shown in Figure 4B and Table 6, CLDN10 expression was correlated with immune markers of macrophages and CD8+T cells. In consideration that macrophages and CD8+T cells might be favorable prognostic indicators for PTC, the data may help explain why low CLDN10 expression is related to poor prognosis in patients with PTC.

Data source

A cohort of PTC was retrieved from the TCGA database [25]. The gene expression profile based on RNA-seq and relevant clinical data were downloaded from the TCGA and UCSC Xena databases, respectively [26].

The stromal and immune scores in the PTC TCGA dataset were calculated based on the ESTIMATE algorithm to predict the level of infiltrating stromal and immune cells in PTC and infer tumor purity in tumor tissue [27].

The data were analyzed by using the package limma of R program (version 3.6.0). The samples were divided into two groups (high and low) on the basis of the immune and stromal scores, and the two groups were compared. The results were downloaded and analyzed, in which the genes that met the cut-off criteria of adjusted P<0.05 and a |log fold-change| of >1.0 were screened.

The screened genes regarding immune and stromal scores were processed with Venn analysis. Genes of the intersection were considered as the DEGs.

Functional annotation of the DEGs

GO and KEGG pathway enrichment analyses were conducted by using the Gather database to learn the possible functions of the DEGs [28].

GO is a comprehensive community-based bioinformatics resource that provides information about gene or gene product function using ontologies to represent biological knowledge [29]. In the GO analysis, the gene network was presented according to biological processes [30]. KEGG is a reference knowledge base for linking sequences to biological functions from molecular to higher levels [31], which is an integrated database consisting of three generic categories of systems information, genomic information, and chemical information [32].

A P-value less than 0.05 was considered statistically significant.

Construction of weighted co-expression network (WGCNA) and identification of modules related to LNM

Relevant gene modules were identified using the WGCNA method to screen the possible gene sets associated with LNM [33]. WGCNA is an algorithm widely used in gene co-expression network identification in a number of disorders to find significantly correlated gene modules. Scale independence and mean connectivity analysis of modules with different power values were conducted to assess the soft threshold of module analysis. The power value was determined when the scale independence value was 0.9. Then, with the minimal module size of 30 and the merge cut height of 0.25, the co-expression matrix was calculated. The results of dynamic tree cut and merge were displayed by clustering dendrogram.

Information regarding the clinical features was used to identify significant co-expression modules associated with the clinical traits. The related gene modules concerning clinical features (such as LMN) were identified by calculating their correlations. Module membership (MM) and gene significance (GS) were calculated in an intramodular analysis of module statistically. Significant modules were defined with a P value less than 0.05.

Key gene screening

The expression levels of the DEGs were divided into high and low groups on the basis of their median levels. From the TCGA data, the univariate cox regression analyses were used to evaluate their prognostic values for PTC. The genes with P values less than 0.05 were selected. Second, the modules that had a close association with LNM were chosen. The genes with a GS of more than 0.15 were selected. Third, Venn analysis was conducted to filter the intersection of the genes selected from the above two steps.

PTC specimens

A PTC tissue microarray (HThy-Pap120CS-01) comprising 58 PTC tissues and 58 paired adjacent non-cancer tissues was purchased from Shanghai Outdo Biotech Co., Ltd.. All patients were pathologically diagnosed as PTC and received no extra treatment before surgery.

IHC staining and evaluation

The protein expression levels of the key genes were tested by using the two-step method of IHC as previously described in our study [34]. In brief, the slides were deparaffinized, rehydrated, and treated with 3% hydrogen peroxide for 20 min to inhibit endogenous peroxidase. The sections were rinsed with distilled water and saturated in phosphate buffered saline for 5 min and then incubated with a 1:200 dilution of rabbit anti-polyclonal antibody (primary antibody; Abcam) overnight at 4 °C. The staining was visualized using DAB solution and counterstained with hematoxylin. IHC staining was conducted in accordance with the manufacturer’s instructions.

The IHC stain results were identified by integrated scoring. The results were evaluated and scored independently by two pathologists without knowledge of the clinical parameters of the cases. The staining intensities of the proteins were scored from 0 to 3, where 0 means negative, 1 weak, 2 moderate, and 3 strong. The percentages of positively stained cells were scored in scales of 0 to 4, in which 1 represents (0–25%), 2 (26%–50%), 3 (51%–75%), and 4 (76%–100%). The proportion and intensity scores were then multiplied to gain a total score, with a range from 0 to 12.

Statistical analysis

For continuous variables, differences between the groups were analyzed using ANOVA, a t-test, or a Wilcoxon rank sum test in accordance with the concrete types of the data. Chi-squared test was used to differentiate the rates of different groups. The diagnostic accuracy of the genes was measured by the ROC curves and the AUC. The Kaplan–Meier method was conducted to calculate the overall survival curves. A log-rank test was used to determine differences in the survival rates. These analyses were performed by utilizing MedCalc software (15.2.2; Mariakerke, Belgium). Statistical significance was considered at P < 0.05.

Ethical conduct of research

The authors state that they have obtained appropriate institutional review board approval or have followed the principles outlined in the Declaration of Helsinki for all human or animal experimental investigations. The present study used a commercial tissue microarray that was purchased from Shanghai Outdo Biotech Co., Ltd.