USP13 genetics and expression in a family with thyroid cancer

Papillary thyroid carcinoma (PTC) is the most common type of thyroid carcinoma and its incidence has greatly increased in the last 30 years. Ubiquitin-specific protease 13 (USP13) is a class of deubiquitinating enzymes (DUBs) and plays an important role in cellular functions such as cell cycle regulation, DNA damage repair, and different cell signaling pathways. Studies regarding the role of USP13 in cancer development and progression are divergent and there are no previous data regarding the role of USP13 gene in PTCs. In this study, we investigated the genetic cause of PTC diagnosed in multiple members of a Brazilian family. Whole exome sequencing (WES) was performed to identify the genetic cause of PTC. Cycloheximide chase assay and clonogenic assay were performed to study USP13 stability and function in vitro. WES analysis identified a heterozygous missense variant c.1483G > A (p.V495M) in the USP13 gene that fully segregates with the disease. In silico modeling suggests that this variant may cause protein structural perturbations. USP13 overexpression increased the potential of a single cell to form colonies. The USP13 c.1483G > A variant enhanced the effects seen in USP13 overexpression and preserved protein stability for longer hours compared to the non-mutated USP13 protein. Our study suggests that USP13 overexpression may play a role in tumorigenesis of PTCs; and that the USP13 p.V495M (c.1483G > A) variant enhances USP13 estability.


Introduction
Papillary thyroid carcinoma (PTC) is the most common type of endocrine malignancy, representing approximately 85% of all differentiated thyroid carcinomas [1,2]. The incidence of thyroid carcinomas has greatly increased in the last 30 years with a prevalence in women about four times higher than in men [3]. This increase is due to highresolution imaging techniques improvement that led to early detection. However, studies also suggest that external elements as lifestyle changes, environmental conditions, radiation exposure, and family history have implications in the increase of thyroid cancer incidence, particularly PTCs [1,4]. Despite the high incidence, mortality rates have decreased over the years due to improved management and treatment of the disease [5]. Yet, high resistance to standard chemotherapy, extensive resection, recurrence, and metastasis can lead to poor outcome [6]. In this context, there is still lack of information regarding genetic alterations that lead to PTC onset and progression [7].
In the advent of genetic sequencing techniques, new possibilities for diagnosis of different types of tumors with hereditary origin have emerged. Whole exome sequencing (WES) showed to be a potent tool to identify existing variants and became more accessible to clinicians to investigate different diseases [6] and therefore, improved early diagnosis and more prognosis. Particular driver pathogenic variants were found to be associated with different histological subtypes of thyroid cancers. Regarding PTCs, classification by gene expression profile includes BRAF-like pathogenic variants: BRAF V600E , RET/PTC fusion, NTRK fusions, or RAS -like molecular subtypes: RAS, BRAF K601E , PAX8/ PPARG, PTEN, EIF1AX; with BRAF V600E being the most common pathogenic variant found in PTCs (40-80%) [8]. Although less common, other genes have showed to have an important role in PTC onset such as DICER1, APC, KTM2D, and WDR77 [9][10][11].
Ubiquitin-specific proteases (USPs) are the biggest class of deubiquitinating enzymes (DUBs) and play an important role in cellular functions such as cell cycle regulation, DNA damage repair, and different cell signaling pathways. Therefore, due to interactions with different cellular targets, USPs have been associated with tumor development and progression [12]. USP13 belongs to the USPs and its function in cancer is controversial. A recent study showed that USP13 may act as an oncogene in cervical cancer [13]. In the same direction, it has been shown a pro-tumor role of USP13 in glioblastoma [14]. On the other hand, in some studies, USP13 gene seems to act as a tumor suppressor by interacting with PTEN signaling in human bladder cancer and oral squamous cell carcinoma [15][16][17]. To our knowledge, there are no prior studies in the literature investigating the role of USP13 gene in thyroid cancers.
In the present work, we investigated the genetic cause of thyroid dysfunctions, including diagnosis of PTCs, Hashimoto thyroiditis, and goiter, presented by several individuals from the same family. Whole exome sequencing identified a heterozygous missense variant c.1483G > A (p.V495M) in the USP13 gene that fully segregates with the disease. Since USP13 function in thyroid development and progression is unexplored and little is known about this specific variant in tumor development and progression, we evaluated whether USP13 plays a role in MDA-T32 thyroid cancer cells and whether the variant affects USP13 functionality and stability in these cells.

Patients
Thirteen members of a large Brazilian family were studied. This family was previously screened for genes associated to cancer, including SDHA, SDHB, SDHC, SDHD, BRAF, JUNB, CDKN1B, and PTEN [18]. In this study, only members without thyroid problems and older than 40 years old were treated as unaffected. All patients were recruited under protocols approved by the Institutional Review Board of the participating institutions and written informed consent was obtained from all participants.

Clinical and DNA studies
The sampling of epithelial buccal cells was performed as previously described [19]. After the collection, 10 μl of proteinase K (20 mg/mL) were added to the solution, being left overnight at 65°C, DNA was purified by adding ammonium acetate 10 M, precipitated with isopropanol, and resuspended with 50 μl of Tris 10 mM (pH 7.6) and EDTA 1 mM. Exome sequencing was conducted as we have described elsewhere [20,21]. In brief, exome DNA was captured using a Sur-eSelect Human all exome v6 (Agilent Technologies, USA). The 150-bp library was sequenced paired-end on the Nova-Seq6000 platform (Illumina), and read alignment was done against the human chromosome reference assembly build 37 (GRCh37, NCBI). Filtering and analysis of the variant calling files was performed using the VarSeq v2.2.2 (Golden Helix Inc., Bozeman, MT) using the hereditary gene panel with modifications adjusted to our studied pedigree. In particular, we looked for variants that can be related with the phenotype in our 13 subjects, making sure that the subjects treated as unaffected do not harbor any of these variants. In order to identify potential causative variants, we applied the following filters: we searched for missence variants, deletions, and insertions that are present in gnomAD with a frequency less that 1%. Variants that were predicted as damaging by more than four in silico tools (see in silico analyses section) were selected for further studies. By applying this filtering, we were able to identify 78 variants. However, the only gene variant that could segregate with the phenotype in all affected member of the family was a variant in USP13 gene. The identified USP13 genetic variant was confirmed through Sanger sequencing in all members of the family carrying it based on the exome sequencing findings. biocompute.org.uk/fathmmMKL.htm) to predict the possible impact of the amino acid substitution on the structure and function of the corresponding proteins. Variants that were predicted as damaging by more than four in silico tools and also have a frequency in gnomAD less than 1% were selected for further studies.

In silico analyses
Ab initio modeling I-Tasser was used for the ab initio modeling of the native (wild type) and the mutated USP13 deduced amino acid sequences of p.V495M (http://zhanglab.ccmb.med.umich. edu/I-TASSER) [9][10][11]. Predicted models were evaluated for discrete optimized protein energy and root mean square distance (RMSD) score (for the measurement of the average distance between the backbones of the superimposed proteins), and the TM score (to assess topological similarity) was calculated using TM-Align [22][23][24]. The chosen 3-dimensional structures of the native and the mutated structures were further analyzed using PyMOL molecular graphics system (DeLano Scientific; http:// pymol.sourceforge.net/).

USP13 siRNA
USP13 siRNA was performed in MDA-T32 cells following manufacturer's protocol. We used Human USP13 on-target plus smart pool siRNA (L-006064-00-0010, Dharmacon, CO, USA) and non-targeting on-target plus control pool (D-001810-10-05, Dharmacon, CO, USA). Briefly, 3 × 10 4 MDA-T32 cells were plated in 12 well plates. The next day, cells were transfected with 1 µM of siRNA or control and incubated for 24 h. Medium was replaced by fresh medium for more 24 h and cells were collected for protein extraction.

Site-directed mutagenesis
The human USP13 wild type (NM_003940) coding sequence was cloned into the pCMV6-Entry vector with C-terminal Myc-DDK tag (RC202190, Origene, Rockville, MD, USA). The p.V495M variant was introduced into the human USP13 wild type template using the QuikChange Lightning Sitedirected Mutagenesis Kit (210518-5, Agilent Technologies, Santa Clara, CA, USA), following the manufacturer's protocol. Insertion of the variant was confirmed by Sanger sequencing. The following mutagenic primers were used: USP13-Val495Met MUT _F: GCATCAGGTAATCC ATCCTCTCCGTGTAGCG USP13-Val495Met MUT_R: CGCTACACGGAGAGG ATGGATTACCTGATGC Cycloheximide chase assay MDA-T32 cells were seeded into 12-well plates at a density of 3 × 10 4 cells per well. After 24 h of incubation, cells were transfected with Effectene Transfection Reagent (301425, Qiagen, Germany) following manufacturer's protocol. 0.5 µg of each vector (USP13 wild type and USP13 p. V495M variant) was used for the transfection. The empty pCMV6-Entry vector was used as a negative control. After 24 h of transfection, medium was removed, cells were treated with 25 µM of cycloheximide (MilliporeSigma, 01810) and incubated for 3, 8, and 24 h. After each incubation time, cells were washed with PBS, and protein was extracted. Experiment was performed in triplicate. System (Licor, USA) were used to acquire the signal of the bands. Densitometric quantification was performed using ImageJ software and a ratio between the protein of interest and GAPDH was calculated and plotted on GraphPad.

Analysis of protein expression
Clonogenic survival assay 1 × 10 3 MDA-T32 cells/well were seeded in six well plates for 24 h and transfected with empty vector (p.CMV6), USP13 wild type, and USP13 p.V495M, as described above, for 48 h. After incubation time, the cell culture medium was removed, and fresh medium was added. Cells were then incubated for 7 days and medium was changed every two days. In the end point of the experiment, medium was removed, and colonies were fixed with a solution containing 0.05% w/v of crystal violet, 1% of formaldehyde, and 1% methanol in PBS. Fixed colonies were photographed at a magnification of 4x. One milliliter of acetic acid 10% was added to each well and plates were shaken for 15 min to discolor the plaque. 200 μL of the violet stain solution removed from the plate was transferred to a 96-well plate and read on a spectrophotometer at an absorbance of 590 nm. The intensity of staining is proportional to the number of colonies established after treatment. The experiment was performed in triplicate.

Cell viability assay
MDA-T32 cells were seeded into 96-well plates at a density of 4 × 10 3 cells per well and incubated for 24 h. Cells were then transfected with the non-target control or siRNA as described above and incubated for 0, 24, 48, and 72 h. At each time point, MTT (3-(4,5-dimethylthiazol-2-yl)-2,5diphenyltetrazolium salt) was added to a final concentration of 0.5 mg/mL into the culture media. After 3 h of incubation with MTT, medium was removed and 200 μL of isopropanol containing acetic acid 0.04 M was added into the wells, and plate was gently shaken to solubilize the product reduced by MTT. Absorbance was measured at 570 nm.

Statistical analysis
Statistical analysis was performed using Student's t-test or One-way ANOVA with multiple comparisons tool. Data are expressed as mean ± SEM. Differences between mean values were considered significant when p < 0.05.

Genetic analysis
Filtering of the variant calling files, performed in the Var-Seq v2.2.2 (Golden Helix Inc., MT, USA) using the hereditary gene panel, revealed a heterozygous missense variant in the USP13 gene c.1483G > A (p.V495M). Initially, we included into the genetic analysis the following patients (Affected: I-1, I-2, II-3, II-4, II-6, II-9; Unaffected: I-3, II-2, II-8, III-2) (Fig. 1). This selection was based on the stringiest inclusion criteria based on the phenotype. The results of this analysis picked the USP13 genetic variant c.1483G > A (p.V495M) out of 78 genetic variants as the only one that segregates with the disease in this family.
Predicted damaging variant in USP13 was found segregating with the disease A heterozygous missense variant was detected in USP13 gene c.1483G > A (p.V495M) that fully segregates with the disease. The genealogic tree of the family shows that  (Fig. 1A). This variant has an allele frequency of 0.000302 (gnomAD Exomes v2.1.1) and it is predicted to be damaging in four out of six in silico analysis tools used (Polyphen-2, MutationTaster, MutationAssessor, and FATHMM MKL). This missense variant is located in the USP domain of the USP13 protein and the mutated amino acid is evolutionarily highly conserved among species ( Fig. 2A).

In silico modeling suggests that USP13 c.1483G > A (p.V495M) variant may cause protein structural perturbations
The TM score of the alignment between that native (wild type) and the mutant protein (p.V495M) was 0.93687 and the RMSD value was 2.98, indicating potential pathogenicity for the identified missense variant. RMSD values greater than 0.15 are considered as an indicator for significant structural perturbations, which could have functional implications for the protein function [22]. The 3-dimensional effect of the p.V495M variant that is located at the USP domain was defined using ab initio modeling. As shown in Fig. 2B, the alignment of the ab initio native and mutant model reveals that the variant at the position 495 may change the folding of the protein as the new amino acid Met contains a sulfur that is predicted to interact with the aromatic residue tyrosine in position 497 as indicated by the dotted line. Moreover, this variant is located in the USP domain, a critical domain for DUBs, and therefore may have an impact in the protein function.

USP13 does not affect PTEN stabilization in MDA-T32 thyroid cancer cells
Studies in different types of tumors already demonstrated that USP13 stabilizes PTEN, the second most frequent mutated gene in human cancers [15-17, 25, 26]. Since there are no records in the literature regarding signaling of USP13 in PTCs, we checked whether USP13 would play a role in PTEN stabilization. We evaluated PTEN expression in USP13 siRNA MDA-T32 cell lines. Our results show that the absence of USP13 does not affect expression of PTEN in this cell lineage (Fig. 3A-C).

Silencing USP13 in MDA-T32 cells decreases cell proliferation and colony formation
The ability of cells to constantly proliferate without an external stimulus is one of the main aspects of cancer cells. The capacity of a single cell to migrate and establish colonies in a distant site of its previous microenvironment is also another important feature acquired by cancer cells [27]. We investigated whether USP13 would have a function on proliferation and clonogenic capacity of thyroid cancer cells. USP13 siRNA was performed in MDA-T32 cells and our results show that reduced expression of USP13 decreases the capacity of cells to form colonies (Fig. 3D, E) and leads to decreased cell proliferation rates after 72 h of USP13 siRNA transfection (Fig. 3F The cycloheximide chase assay enables to observe degradation kinetics of proteins and therefore is a good tool to study protein stability [28]. By treating wild type and USP13 c.1483G > A mutated cells with cycloheximide, an inhibitor of protein synthesis, we were able to evaluate the dynamics of USP13 degradation when MDA-T32 thyroid cells present the variant. After 8 h of cycloheximide treatment, expression of USP13 in wild type cells is significantly decreased, while in cells transfected with the USP13 c.1483G > A variant, this effect is seen after 24 h of cycloheximide treatment (Fig. 4B, C). These results suggest that the presence of the variant have a significant impact on the USP13 stability and consequently on its targets, as USP13 will be active in the cells for longer periods.

USP13 c.1483G > A variant increases MDA-T32 thyroid cancer cells colony formation
As studies available in the literature regarding the role of USP13 in cancer development and progression are controversial, with evidence of USP13 acting as an oncogene or as a tumor suppressor, depending on the type of cells, it is hard to predict whether this variant would be, in fact, contributing to the survival of thyroid cancer cells. To answer this question, we evaluated the effect of the variant in the capacity of colony formation of MDA-T32 cells. Our results show that the overexpression of USP13 gene increased the capacity of the cells to form colonies. Interestingly, the presence of the variant, USP13 c.1483G > A (p.V495M), drastically increased this effect (Fig. 4D, E), suggesting that, in thyroid cancer cells, USP13 might act as an oncogene and, in the presence of the variant, the clonogenic capacity of the cells are enhanced.

Discussion
Attempting to investigate the genetic cause of PTC in multiple members of the same family, we found an association of the USP13 c.1483G > A; p.V495M variant with PTC development. Individuals that carry the variant either developed PTC or other thyroid dysfunctions such as Hashimoto thyroiditis or goiter, suggesting that the onset of these thyroid dysfunctions could be related to genetic factors [29]. It is important to mention that carriers of the variant that did not developed PTC may develop it later, since the age of onset may not have been reached yet.
Ubiquitination is an essential post-translational modification (PTM) event that regulates protein function, localization, and turnover, contributing to sustain cellular homeostasis. The ubiquitination process can be reversed by USPs. Thus, USPs play an important role in cell signaling, control of protein degradation, and tumorigenesis [12,30]. USP13 is a subtype of USP and its role in tumor development and progression is controversial. Qi and co-workers performed WES to study the somatic mutation profile of PTCs in a Chinese cohort. They found mutations in relevant genes such as BRAF, BCR, TP53, between others. Insertion of USP13 was found in the cohort, but it was not further explored in the study [31]. Another WES study, performed in a cohort of 158 patients diagnosed with PTC, found APOBEC SBS13 mutational signature as predictor of severity of PTCs [32]. In a study published by Zhao and coworkers, WES analysis was performed in individuals of the same family diagnosed with PTCs. They identified a variant in the WDR77 gene that led to predisposition to PTCs. Based on the analysis performed, nothing was shown regarding USP13 regulation [11]. In our study, WES analysis identified a variant in the USP13 gene (c.1483G > A; p.V495M) that was present in multiple members of the same family diagnosed with PTC.
Researchers attributed a tumor suppressor role to USP13 gene by showing that USP13 mediates PTEN signaling in different types of cancers [15][16][17]. In our study, silencing USP13 gene in MDA-T32 thyroid cancer cell lines did not interfere on PTEN expression (Fig. 3A-C). Therefore, other pathways could be involved in the USP13 signaling of MDA-T32 cells. On the other hand, decreased expression of USP13 in MDA-T32 thyroid cancer cells led to a less severe phenotype, such as cell viability weakening and decreased clonogenic capacity, suggesting that USP13 may play a protumor role in PTCs (Fig. 3D-F).
Due to the segregation of the USP13 c.1483G > A; p. V495M variant within individuals that presented thyroid disfunctions reported in this study, and the putative protumor role of USP13 seen in vitro, we hypothesized that it could be a gain-of-function variant. In fact, we found that the presence of the variant increased USP13 half-life compared to the wild-type protein in MDA-T32 cells (Fig. 4B, C). Based on available data in the literature, the role of USP13 in tumorigenesis is ambiguous. A few studies suggest a tumor suppressor role for USP13: USP13 showed to stabilize PTEN, a potent tumor suppressor, in breast, bladder cancers, and oral squamous cell carcinomas [15][16][17]. On an opposite direction, other studies show that USP13 may act as an oncogene. For instance, USP13 stabilizes c-Myc that leads to an increased proliferation of glioma cells [14] and Mcl-1, which seems to increase proliferation of lung, ovarian and cervical cancer cells [13,33]. Moreover, silencing USP13 decreased cell proliferation in lung metastasis of hepatocellular carcinoma [34] and its overexpression is correlated with poor prognosis and chemoresistance of ovarian cancer [35]. USP13 was also showed to interact with aurora kinase B and disrupt cell cycle regulation [36].
Our in vitro study showed that the USP13 c.1483G > A; p.V495M variant preserves protein stability in MDA-32 thyroid cancer cells compared to wild type USP13. Supporting the pro-tumor role of USP13 previously described in other types of cancers, we demonstrated that USP13 overexpression increases the capability of a single cell to form colonies and that the c.1483G > A; p.V495M variant enhances the effects seen in the USP13 overexpression, causing an even more severe phenotype compared to the wild type protein (Fig. 4D, E) goiter. 2. In vitro assays suggest that the homozygous USP13 c.1483G > A can favor tumor growth. By performing the in vitro assays, we were able to better understand the role of USP13 in PTCs, while seeing heterozygous patients developing the disease opens a discussion on the importance of genetic screening for early diagnosis of PTCs.
A limitation of the study is the availability of tissues to analyse the prevalence of USP13 variants in individuals with PTCs, as well as compare USP13 expression in neoplastic tissues with contralateral non-neoplastic samples. Despite our in silico and in vitro analysis led us to the believe that USP13 might be acting as an oncogene and it goes in the same direction of the data found in patients, it is important to emphasize that the sample size of this study is small and requires further investigation using in vivo models and different cohorts, to ensure the oncogenic role of USP13 in PTCs and the putative use of USP13 as a marker for the disease. Because the incidence of PTCs is constantly increasing [3], this study may have an impact on early diagnosis and prognosis of the disease, opening new avenues for follow up of individuals that carry USP13 genetic variants in regard to thyroid dysfunction. Taken together, this study brings the first insight that USP13 might be an oncogene in the context of PTCs and USP13 aberrant expression may also cause thyroid dysfunctions such as Hashimoto thyroidits and goiter.

Compliance with ethical standards
Conflict of interest Dr. Stratakis holds patents on the function of the PRKAR1A, PDE11A, and GPR101 genes and related issues; his laboratory has also received research funding on GPR101 and its involvement in acromegaly and/or gigantism, abnormal growth hormone secretion and its treatment by Pfizer, Inc.; Dr. Faucz holds patent on the GPR101 gene and/or its function; The other authors have nothing to disclose.
Publisher's note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.