HLA-G Gene Variability Is Associated with Papillary Thyroid Carcinoma Morbidity and the HLA-G Protein Profile

Human leukocyte antigen (HLA)-G is an immune checkpoint molecule that is highly expressed in papillary thyroid carcinoma (PTC). The HLA-G gene presents several functional polymorphisms distributed across the coding and regulatory regions (5′URR: 5′ upstream regulatory region and 3′UTR: 3′ untranslated region) and some of them may impact HLA-G expression and human malignancy. To understand the contribution of the HLA-G genetic background in PTC, we studied the HLA-G gene variability in PTC patients in association with tumor morbidity, HLA-G tissue expression, and plasma soluble (sHLA-G) levels. We evaluated 185 PTC patients and 154 healthy controls. Polymorphic sites defining coding, regulatory and extended haplotypes were characterized by sequencing analyses. HLA-G tissue expression and plasma soluble HLA-G levels were evaluated by immunohistochemistry and ELISA, respectively. Compared to the controls, the G0104a(5′URR)G*01:04:04(coding)UTR-03(3’UTR) extended haplotype was underrepresented in the PTC patients, while G0104a(5′URR)G*01:04:01(coding)UTR-03(3′UTR) was less frequent in patients with metastatic and multifocal tumors. Decreased HLA-G tissue expression and undetectable plasma sHLA-G were associated with the G010102a(5′URR)G*01:01:02:01(coding)UTR-02(3′UTR) extended haplotype. We concluded that the HLA-G variability was associated with PTC development and morbidity, as well as the magnitude of the encoded protein expression at local and systemic levels.


Introduction
Thyroid cancer (TC) is the most common endocrine malignancy, and its incidence has markedly increased during the past few decades [1].Differentiated thyroid carcinomas (DTCs) are the most frequent ones, originating from thyroid follicular cells.Among DTCs, papillary thyroid carcinoma (PTC) comprises about 80% of cases [2].Although the majority of PTC patients present a good prognosis when adequately treated with the standard therapy (i.e., surgery and selective use of radioactive iodine (RAI)) [3], some patients present an aggressive tumor behavior [4].Due to the frequent inconclusive diagnosis of TC, a substantial number of people undergo unnecessary partial or total thyroidectomy, with them being exposed to an avoidable risk [5].
Over the past two decades, considerable progress has been made in the understanding of the molecular mechanisms associated with DTC development.The elucidation of genetic and epigenetic events affecting key signaling pathways, such as the mitogenactivated protein kinase (MAPK) pathway, and the role of these mediators in thyroid carcinogenesis has brought beneficial clinical management to DTC patients [6].More recently, the development of kinase inhibition-based therapies, which allow for the targeting of tumor cells, are considered a novel tool for personalized medicine in TC [7].Despite this progress, there are still limitations concerning the detection of malignant cases and the prediction of unfavorable outcomes [4].Studies on TC have reported several candidate genes, including those associated with the regulation of inflammation and the immune response [8].
Human leucocyte antigen-G (HLA-G) is a unique non-classical immunomodulatory major histocompatibility complex (MHC) class Ib molecule [9] that is physiologically expressed in the thymus, cornea, pancreas, erythroid and endothelial precursors, and in trophoblastic cells at the maternal-fetal interface.In pathological conditions, HLA-G is expressed in a great number of tumors [10], particularly in PTC [11,12].The mechanisms that stimulate HLA-G expression in the context of cancer are complex, depending on the individual's genetic background [13] and various environmental factors such as hypoxia, stress, hormones, and specific cytokines [14].HLA-G expression in tumor cells is deleterious to the host because of its interaction with leukocyte receptors, particularly ILT-2 and ILT-4, which downregulate the cytotoxic activity of CD8 + T and natural killer (NK) cells and inhibit antigen presentation, lymphocyte proliferation, and antibody formation [15,16].Alternative splicing membrane-bound (HLA-G1, -G2, -G3, and -G4) and soluble (HLA-G5, -G6, and -G7) HLA-G isoforms may also present immunomodulatory properties [17].
The main feature that distinguishes HLA-G from classical class I genes is the limited variability at the coding region.The overall structure of the HLA-G alleles has been maintained during the evolution process, in which most of its variable sites are synonymous mutations or are observed in introns [18].Currently, 117 alleles encoding 38 different proteins have been reported by the International Immunogenetics Database (IMGT, 3.52.0)[19][20][21].Nucleotide/amino acid variability at the coding region may produce conformational changes in the molecule, which may modify its principal functions, i.e., interaction with cell receptors, isoform production, modulation of the immune response, and the ability to couple peptides [13].On the other hand, nucleotide variability at the HLA-G regulatory regions (5 upstream regulatory region: 5 URR and 3 untranslated region: 3 UTR) may potentially influence the targeting of transcriptional and post-transcriptional elements, affecting the HLA-G gene and protein expression [22,23].Accordingly, these variation sites along the HLA-G gene have been reported in pre-eclampsia [24], cancer [25,26], transplants [27], and autoimmune [28,29], chronic inflammatory [30], and infectious disorders [31,32].In some of these conditions, the HLA-G gene variability potentially affected the magnitude of HLA-G expression [23,33,34].
Considering that (i) HLA-G is a well-recognized immune checkpoint molecule [9], (ii) HLA-G is expressed in several thyroidal disorders, especially in PTC, but it is not in histologically normal thyroid specimens [12], (iii) the HLA-G gene has several nucleotide variation sites at the regulatory regions that have been associated with the magnitude of molecule production [18], and (iv) the study of the HLA-G variation sites in PTC has mostly focused at the 3 UTR of the gene [35], we evaluated the complete HLA-G gene variability (5 URR, coding, 3 UTR, and extended haplotypes) in PTC patients in association with the clinical, histopathological, and tumor features and with the HLA-G protein expression profile.We identified HLA-G-specific regions and/or extended haplotypes associated with: (i) PTC histological and prognosis features, (ii) HLA-G expression in tumor specimens, and (iii) soluble plasma HLA-G detection before thyroidectomy.

HLA-G Gene Variability in PTC Patients and Control Subjects
For HLA-G gene nucleotide variability analyses, we evaluated 151 PTC patients and 154 unrelated healthy controls from the same geographical region.All detected HLA-G polymorphic sites fit the Hardy-Weinberg expectations and present strong linkage disequilibrium (LD) (p < 0.05), indicating that the alleles do not segregate randomly, justifying the haplotype inferences.Association analyses were conducted considering the three major HLA-G gene segments, taken individually (5 URR, coding, and 3'UTR) (Table S1), and as extended (5 URR/coding/3 UTR) haplotypes (Table 1).

HLA-G Gene Variability and Histopathological Features of PTC
To assess whether the HLA-G gene variability is associated with histopathological features related to PTC severity and prognosis, we analyzed the whole group of 185 PTC patients from two different geographical regions of Brazil (Table S2 for the HLA-G gene segments individually; Table 2 and Table S3 for the HLA-G extended haplotypes).

Discussion
Most TCs are differentiated tumors and, therefore, are associated with a low mortality rate [4].However, 20-30% of these patients develop persistent or recurrent disease, most commonly in association with cervical and locoregional lymph nodes (LNs), requiring additional therapy and surgical intervention, which may increase morbidity [37].The identification of patients presenting risk factors associated with aggressive disease can guide tumor management and prevent the overtreatment of patients with low-risk tumors [4].Thus, many studies aiming to understand the role of genetic features in tumorigenesis and their effects on clinical DTC outcomes have so far been conducted to assess pre-operative diagnosis, prognosis, and treatment [38].
HLA-G gene polymorphisms have been studied in several types of human malignancies and have been associated with risk and predictive factors, and most of these studies are limited to evaluating single polymorphisms at the HLA-G regulatory gene segments (i.e., 5 URR and/or 3 UTR) [34].Although many polymorphic sites along the HLA-G gene are in LD, a single polymorphism or a specific coding region haplotype may be associated with distinct regulatory region haplotypes [39], indicating that the typing of the entire gene is more adequate to investigate the HLA-G gene's association with human tumors [23].Indeed, in this study, we evaluated the HLA-G gene regions separately and as extended haplotypes in terms of their association with HLA-G tissue expression, sHLA-G plasma detection, and PTC outcomes.
The neo-expression of HLA-G in human tissues that constitutively do not express the molecule, as is the case of the thyroid, depends on the combination of factors present in the tumor milieu, such as environmental and hormonal factors, as well as genetic factors, including transcription and post-transcriptional factors [14].These transcriptional and post-transcriptional factors may be differentially affected by the polymorphisms located at the 5 URR and 3 UTR of the HLA-G gene, respectively [9], yielding differential HLA-G expression in physiological and pathological conditions [40].The HLA-G coding variability may affect the differential pattern of alternative splicing, leukocyte-receptor interaction, and peptide coupling [41,42].We and other groups have previously reported that HLA-G is commonly expressed in PTC specimens [11,12], while sHLA-G plasma levels are frequently absent or detected at low levels [36,43].The contrasting profile between the high HLA-G expression in the tumor microenvironment and the low plasma sHLA-G levels in PTC patients is still poorly understood.Both local modulatory factors and individual genetic features may account for the increased expression of HLA-G in the tumor microenvironment.
Undoubtedly, the HLA-G 3 UTR gene segment is the most studied one in terms of disease susceptibility and morbidity, particularly with respect to the presence (insertion or "Ins") or absence (deletion or "Del") of a 14-base-pair (14-bp) fragment, in which the 14-bp Del has been associated with increased HLA-G expression, whereas the 14-bp Ins has been associated with decreased HLA-G expression [44].In this series, we reported that the HLA-G UTR-01 haplotype was overrepresented in PTC patients compared to controls.Since UTR-01 contains the 14-bp Del and is associated with increased sHLA-G levels even in healthy individuals [45,46], patients carrying the UTR-01 haplotype may have an additional risk of PTC development.In contrast, the UTR-03 haplotype that also contains the 14-bp Del was associated with protection against PTC development, as well as the G0104a (5 URR) G*01:04:04 (coding) UTR-03 (3'UTR) extended haplotype.Considering that the HLA-G 3 UTR and the 5 URR segments present several elements that have been functionally studied, exhibiting potential roles in HLA-G production [40,45], the evaluation of a single polymorphism inside a specific gene region or a single gene region may only explore partial content of the HLA-G and disease association.
With respect to the association of HLA-G variability with clinical features, we observed that PTC patients carrying the (i) 5 URR G0104a haplotype, (ii) 3 UTR UTR-03 haplotype, and (iii) G0104a (5 URR) G*01:04:01 (coding) UTR-03 (3 UTR) extended haplotype exhibited clinical benefit with respect to PTC development, since their frequencies were underrepresented in patients exhibiting metastatic and multifocal tumors.Cervical lymph node metastases (LNMs) are the main way through which PTC spread cancer cells [47] and present a worse prognosis due to a more frequent recurrence rate [48].Tumor recurrence in the LNs after primary PTC treatment had an independent and significant negative effect on survival in patients aged over 45 years [47].Although PTC metastasis is usually uncommon, it is associated with high mortality [49].Nevertheless, the presence of neck LNM upon diagnosis was one of the most important factors associated with a greater risk of having distant metastases in a Brazilian cohort [50].Additionally, multicentricity also represents a significant risk factor for disease progression in TC and increases the risk of disease recurrence [51].In contrast, UTR-02 was associated with a risk of multicentricity, which is in agreement with our previous finding evaluating a smaller cohort of PTC patients [35].
Regarding the association between HLA-G gene variability with HLA-G tissue expression, we observed that the coding G*01:01:02:01 and the 3 UTR UTR-02 haplotypes individually or assembled as the G010102a (5 URR) G*01:01:02:01 (coding) UTR-02 (3 UTR) extended haplotype exhibited an association with low HLA-G expression in PTC specimens.Individuals carrying UTR-02 have already been described as low sHLA-G producers [46].In contrast, the 5 URR G010101c haplotype and the G010101c (5 URR) G*01:01:01:05 (coding) UTR-04 (3 UTR) extended haplotype were associated with greater HLA-G expression in the PTC specimens.While UTR-04 has been associated with intermediate plasma sHLA-G levels in healthy Brazilian and French individuals [45], functional studies have reported that cells transfected with distinct HLA-G 5 URR yielded differential responses in the presence of HLA-G inducers or repressors.For instance, an HLA-G + cell line (choriocarcinoma JEG-3) transfected with the G010101c haplotype significantly increased the luciferase activity cultured in the presence of interferon (IFN)-β when compared to non-treated transfectants [40].Therefore, as mentioned previously, distinct HLA-G 5 URR, under the influence of diverse transcription factors may differentially induce gene expression.The balance between transcriptional (primarily acting at the 5 URR segment) and post-transcriptional (primarily acting at the 3 UTR segment) factors present in the PTC microenvironment may determine the magnitude of HLA-G tumor expression, corroborating the need to evaluate the entire gene variability.
Concerning sHLA-G detection, we observed that UTR-18 was associated with detectable sHLA-G, whereas the G010102a (5 URR) G*01:01:02:01 (coding) UTR-02 (3 UTR) extended haplotype was associated with undetectable sHLA-G expression.The HLA-G UTR-02 haplotype is composed of the 14-bp Ins and +3142G; both variations were previously associated with low mRNA availability and decreased HLA-G expression, whereas the UTR-18 haplotype is composed of the 14-bp Del and +3142C associated with high HLA-G expression [45,52].Additionally, an HLA-G + cell line (melanoma FON + ) transfected with the G010102a haplotype exhibited the lowest promoter activity among all other 5 URR haplotypes when compared to the empty vector [40].
A recent paper published by our research group characterized the entire HLA-G genetic variability in 4640 individuals from 88 different population samples across the globe from four datasets [53] .These alleles represent 91% of all HLA-G sequences detected worldwide, and their frequencies differ in each biogeographic region [53].All of these HLA-G alleles were found in the present study and represent 96.9% of all detected haplotypes among the Brazilian subjects from the same geographical region.
In our study, three isolated coding region haplotypes were also associated with PTC: (i) the G*01:04:04 haplotype encodes the HLA-G*01:04 molecule and was associated with protection against disease development, (ii) the G*01:01:02:01 (HLA-G*01:01 molecule) and the G*01:04:01 (HLA-G*01:04 molecule) haplotypes were associated with poor biochemical or structural evidence of persistent disease, and (iii) the G*01:01:02:01 (HLA-G*01:01 molecule) haplotype was associated with undetectable sHLA-G.Although the HLA-G molecule may couple peptides in its groove, its major role is the modulation of immune responses through its binding to leukocyte receptors located outside the peptide-binding groove, whose interactive residues are conserved among these molecules [54].The presence of a few proteins frequently observed worldwide corroborates the role of HLA-G as an immune checkpoint molecule rather than as an antigen-presenting molecule [53].Thus, three major reasons may account for these associations: (i) a hitchhiking effect due to the LD between HLA-G coding and regulatory regions, (ii) different coding region haplotypes may generate differential alternative splicing mechanisms, yielding distinct HLA-G isoforms, and (iii) different coding region haplotypes (e.g., the G*01:04:01 and G*01:04:04) encoding the same molecule (HLA-G*01:04) may present distinct regulatory region haplotypes [39].
One limitation of the present study is that a few associations exhibited a very broad 95% CI and/or exhibited a 95% CI that includes the number 1, particularly involving less frequent HLA-G haplotypes.In human genetic studies, sample size plays a critical role in conducting reliable analyses of polymorphisms [55].Thus, the OR and its associated 95% CI are valuable statistical measures to assess the strength and accuracy of associations between variables.These associations suggest that, for those haplotypes, the sample size of the study is not sufficient to produce a tight 95% CI and an OR measure that is truthfully estimated.Consequently, obtaining a sufficiently larger sample size would be necessary for more precise estimations.On the other hand, most of the significant associations presented a narrow 95% CI, a finding that suggests that our sample size is sufficient to produce statistically significant OR measures, supporting the reliability of the present results.
In contrast, using a Brazilian population is an important strength of our study for its remarkable and unique genetic diversity, resulting from miscegenation among Brazilian Native Amerindian, African, and European populations [56].Genetic studies in a diverse population can uncover unique risk or protective variants that might not be apparent in more homogeneous populations.This knowledge is essential for advancing population medicine and/or personalized medicine, including but not limited to: (i) risk prediction, (ii) disease susceptibility, diagnosis, or prognosis, (iii) drug therapy optimization to increase treatment effectiveness and minimize adverse effects, and (iv) public health strategies for disease prevention and health promotion [57,58].It is also important to emphasize the regional variations within Brazil.There are ancestry differences between regions since the degree of European, African, and Amerindian contributions are not homogeneous [59], which can have implications for genetic studies.Taking these geographical variations into account improves the relevance and accuracy of our findings since the risk association analysis for PTC development was performed evaluating only a specific Brazilian region.

Studied Population
The study was conducted on 185 patients with PTC (age at diagnosis 47.2 ± 15.8 years old; 154 female = 83.2%)submitted to thyroidectomy due to or the suspicion of PTC following the ultrasound-guided fine-needle aspiration (FNA) malignant cytology of the thyroid nodule (V and VI categories; 2017 Bethesda system for reporting thyroid cytopathology) [60].A total of 34 cases were recruited from the "Hospital Liga Norte Riograndense Contra o Câncer", Natal-RN, Brazil, and 151 were recruited from "Hospital das Clínicas de Ribeirão Preto", Ribeirão Preto-SP, Brazil.The patients were followed up during a median period of 67 months (interquartile range, 45.7 to 125.3 months).Disease evolution during the follow-up was monitored by clinical examination, ultrasound imaging, and biochemical measurement of serum thyroglobulin (Tg) and anti-Tg antibody (TgAb) levels.
The cytological and histopathological analyses were performed by experienced pathologists, as well as the tumor node metastasis (TNM) staging classification according to the eighth edition of the American Joint Committee on Cancer (AJCC)/Union for International Cancer Control (UICC), ranging from stages I to IV to predict the risk of mortality from DTC.Overall, the lower the TNM stage, the lesser probability of cancer spread and the lower the risk of mortality [61].
Clinical and histopathological data of all patients (i.e., age at diagnosis, tumor size, histological subtype, tumor invasion, metastasis at diagnosis, multicentricity, extrathyroidal extension, Hashimoto's thyroiditis, and TNM staging) as well as treatment details, including surgery approach and RAI therapy), are summarized in Table 5.
The individual response to initial conventional therapy of each patient at the last follow-up was assessed according to the criteria proposed by Momesso and Tuttle 2014 [62], adopted by the 2015 American Thyroid Association (ATA) guidelines [4], taking into account patients treated with thyroid lobectomy or total thyroidectomy, accompanied or not by RAI ablation.This assessment is a risk stratification system to predict recurrence/persistence of the disease based on four categories: (i) excellent response (no clinical, biochemical, or structural/imaging evidence of disease), (ii) indeterminate response (nonspecific biochemical or structural findings that require continuous observation), (iii) incomplete biochemical response (abnormal Tg or rising TgAb levels in the absence of localizable disease), and (iv) incomplete structural response (persistent or newly identified loco-regional or distant metastases, exhibiting abnormal Tg or TgAb or not).
A group of 154 unrelated healthy blood donors (aged 34 ± 11; 84 women = 54.6%) with no personal or family history of thyroid diseases and exhibiting negative anti-peroxidase antibodies was recruited from "Hospital das Clínicas de Ribeirão Preto", Ribeirão Preto-SP, Brazil.All subjects gave written informed consent to participate in the study, which was approved by local human research ethics committees (processes #18022/2014, and #029/029/2015).
The sequences were analyzed and compared to the reference sequence recognized by the IMGT, <http://www.ebi.ac.uk/ipd/imgt/hla/>URL (accessed on 10 November 2020).All polymorphisms detected were individually annotated, and haplotypes of each HLA-G region were initially defined.These haplotypes were further used to construct HLA-G extended haplotypes, as previously described [39,53].

HLA-G Expression
The HLA-G expression in the tumor microenvironment of the PTC specimens was evaluated in the same series of patients using an immunohistochemical assay as previously described [11].In short, we used the primary monoclonal antibody (mAb) anti-HLA-G MEM-G/2 (EXBIO, Vestec, Czech Republic) that recognizes all soluble and membranebound HLA-G molecules, diluted at 1:100 [64].From the 185 PTC patients whose DNA was sequenced for HLA-G gene variability in the present study, the HLA-G expression in the tumor was available for 170.According to the percentage of stained tumor cells, HLA-G expression was stratified into two categories: low expression (≤50% of positive tumor cells in the PTC specimens) and high expression (>50% of positive tumor cells in the PTC specimens).
Plasma sHLA-G levels were previously evaluated in the same series of PTC patients before thyroidectomy to evaluate its relationship with several cytokines [36].We detected sHLA-G levels by using a specific sandwich ELISA [45], using MEM-G/9 mouse antihuman HLA-G mAb (1:100, EXBIO, Praha, Czech Republic) as the capture antibody and rabbit anti-human β2-microglobulin (1:10.000,DAKO, Santa Clara, CA, USA) as the detection antibody.From the 185 PTC studied patients, soluble HLA-G levels were available for 84.Samples with sHLA-G levels below the detection limit were conventionally expressed as 0 ng/mL and categorically stratified into two categories: undetectable (0 ng/mL) or detectable sHLA-G levels.

Statistical Analysis
The LD among polymorphic sites at the HLA-G gene was evaluated by Haploview ® 4.2 the program [65].Given the positive LD between the polymorphic sites, but unknown gametic phase, the PHASE algorithm [66] was used to infer the most likely haplotypes for each sample.Adherences of genotypic proportions to expectations under the Hardy-Weinberg equilibrium were tested by the exact test of Guo and Thompson [67] using ARLEQUIN 3.5 [68].A two-tailed Fisher's exact test was used for all comparisons, and OR and 95% CI were estimated.A 5% level of significance (α = 0.05) was considered for rejection of the null hypothesis.

Conclusions
The immune checkpoint HLA-G molecule is highly expressed in PTC specimens, although the impact of the HLA-G gene variability on PTC morbidity and/or the regulation of HLA-G expression in PTC was not fully understood until now.Therefore, by evaluating the entire HLA-G gene in PTC patients together with HLA-G expression in tumor specimens and sHLA-G plasma detection, we uncovered important associations between the HLA-G gene, taken as gene segments and/or as extended haplotypes, and: (i) the susceptibility to PTC development, (ii) tumor morbidity, particularly metastatic or multifocal tumors and poor response to conventional therapy, (iii) the magnitude of HLA-G expression in tumor specimens, and (iv) the presence or absence of detectable sHLA-G in pre-thyroidectomy plasma.Although several local and systemic (cytokines, hormones, hypoxia, etc.) factors contribute to both PTC development and HLA-G expression modulation, the associations reported in the present study can have significant implications in the context of PTC outcomes, particularly the identification of specific genetic markers that can serve as potential targets for personalized therapies in the future.

Supplementary Materials:
The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/ijms241612858/s1.Informed Consent Statement: Informed consent was obtained from all subjects involved in the study.

Data Availability Statement:
The data can be requested from the corresponding authors.

Table 1 .
HLA-G extended haplotypes in papillary thyroid carcinoma (PTC) patients and control subjects.

Table 2 .
Association between HLA-G extended haplotypes and histopathological features of PTC (n = 185).

Table 3 .
HLA-G extended haplotypes in papillary thyroid carcinoma (PTC) patients according to the expression of HLA-G in the tumor microenvironment [n = 170].

Table 5 .
Clinical and histopathological features of the papillary thyroid carcinoma (PTC) patients.

Table 5 .
Cont.Including minimal (immediate perithyroidal soft tissues or sternothyroid muscle typically detected only microscopically) and extensive (subcutaneous soft tissues, i.e., larynx, trachea, esophagus, or recurrent laryngeal nerve) extrathyroidal extension.b At initial diagnosis.Tumor node metastasis (TNM) system refers to the primary tumor, lymph node metastasis, and distant metastasis.c Classification at the final follow-up. a