INTRODUCTION {#s1}
============

Papillary thyroid cancer (PTC) is an endocrine malignancy that accounts for almost 90% of all thyroid malignancies (SEER Cancer statistics review;<http://seer.cancer.gov/statfacts/html/thyro.html>) with 64,300 estimated cases and 1,980 deaths in 2016 in the United States. Based on case-control studies PTC has one of the highest family risk ratios of all cancers \[[@R1], [@R2], [@R3]\]. Despite this, no mutations have been implicated to account for more than a tiny fraction of all PTC (see e.g. \[[@R4]\] for a review).

Micro RNAs (miRNAs) are an abundant class of short, noncoding RNAs that are widely expressed in mammalian cells and regulate the translation of protein-coding genes by binding to miR-specific sequences in the 3' untranslated regions. This binding may lead to the degradation of the mRNA thus further downregulating gene expression. Strong overexpression of several miRNAs as well as modest downregulation of other miRNAs in PTC tumor tissue has been reported \[[@R5]\].

Sequence variants in microRNA genes were described to impact miRNA expression and function leading to human disease (for review see \[[@R6]\]). Most miRNAs are initially transcribed from the genome as long primary transcripts (pri-miRNAs), which are processed into precursor miRNAs (pre-miRNAs), and later into mature miRNAs. While it is known that sequence variants in the seed regions of miRNAs can affect binding of the miRNA to the target gene, it has also been shown that variants in the pre-miRNA sequence can affect the processing of the mature miRNA \[[@R7], [@R8]\]. The resulting change in the expression levels of miRNAs causes deregulation of protein-coding genes. As an example, a SNP in miR-146a predisposes to PTC via downregulation of the THRB gene \[[@R9]\].

We reasoned that there might be germline variants in miRNAs that affect the expression of thyroid target genes \[[@R10]\]. In an effort to find miRNA variants that could thus cause PTC we sequenced the genomic regions surrounding 427 selected miRs in PTC patients with a family history of the disease.

RESULTS {#s2}
=======

Identification of miRNA variants in PTC patients {#s2_1}
------------------------------------------------

DNA from 80 individuals affected with papillary thyroid cancer (PTC) from 80 families was subjected to SOLiD sequencing. The sequencing targeted 321 genomic regions that encode 427 miRNAs expressed in thyroid \[[@R11]\]. The GATK variant discovery pipeline was used according to best practices ( <https://www.broadinstitute.org/gatk/guide/best-practices.php>) and the short reads were aligned to hg19. We identified 1978 variants with a frequency of \< 1% in 1000 Genomes Phase 1 database. To narrow down the number of variants subjected to further tests, only those (*n* = 195) that did not occur in the dbSNP database (dbSNP build 138) were considered (variants not contained in dbSNP are rare). After eliminating variants located outside of the pre-miRNA (*n* = 170), only 25 candidates were left.

Sanger sequencing and cosegregation analysis {#s2_2}
--------------------------------------------

We then proceeded with Sanger sequencing of DNA from the individuals used for SOLiD sequencing in order to confirm the variants. Out of the 25 variants, only one in mir-450a-1 was not confirmed (Table [1](#T1){ref-type="table"}).

###### miRNA variants found in PTC families

  Var      Families (\#)   miRNA        Chromosomal position   Ref   Var   dbSNP         MAF (%) (1K/ESP)^b^   position (in pre-miR)       Cosegr. Analysis^c^
  -------- --------------- ------------ ---------------------- ----- ----- ------------- --------------------- --------------------------- ---------------------
  **1**    2               miR-200b     chr1:1102563           G     A     rs72563729    1.9/2.0               stem, not in mature mir     no
  **2**    1               miR-2355     chr2:207974775         T     C     NA            NA                    middle 5p                   no
  **3**    1               miR-15b      chr3:160122421         A     G     rs146020563   0/0.04                loop                        **yes**
  **4**    1               miR-15b      chr3:160122458         A     G     rs369598613   NA/0.1                stem, not in mature mir     NA
  **5**    1               miR138-1     chr3:44155749          T     C     rs374104338   0.1/0.2               loop, next to 3'end of 5p   no
  **6**    1               miR-135a-1   chr3:52328253          C     A     rs201615303   NA/0                  3' end of 3p                **yes**
  **7**    1               miR-449a     chr5:54466362          T     C     rs144949536   0.7/0.3               stem, not in mature mir     no
  **8**    1               miR-4286     chr8:10524525          G     C     rs566317227   0.1/NA                loop                        NA
  **9**    1               miR-320a     chr8:22102519          G     A     rs200301891   0.2/0.3               loop                        **yes**
  **10**   1               miR-181b-2   chr9:127456023         G     \-    NA            NA                    3' end of 5p                **yes**
  **11**   1               miR-199b     chr9:131007001         C     A     NA            NA                    stem, not in mature mir     NA
  **12**   1               miR-495      chr14:101500157        C     T     rs199783139   NA/0.04               3' end of 3p                NA
  **13**   1               miR-154      chr14:101526140        C     T     rs561595554   0/NA                  loop                        NA
  **14**   1               miR-484      chr16:15737177         G     A     rs376603273   0.1/0.1               3' end of 5p                **yes**
  **15**   1               miR-328      chr16:67236292         C     G     rs188892061   NA/0.7                5' end of 5p                NA
  **16**   1               miR-152      chr17:46114572         C     G     rs200114569   0.6/0.4               loop                        NA
  **17**   2               miR-187      chr18:33484792         G     A     rs41274312    1.4/1.3               stem, not in mature mir     no
  **18**   5               miR-27a      chr19:13947296         G     A     rs11671784    2.4/2.1               loop of pre-mir             no
  **19**   1               miR-181d     chr19:13985711         C     T     rs367686185   NA                    stem, not in mature mir     no
  **20**   1               let-7e       chr19:52196093         A     G     NA            NA                    5' end of 3p                **yes**
  **21**   1               miR-499      chr20:33578202         G     A     rs140486571   0.3/0.7               stem, not in mature mir     no
  **22**   2               miR-296      chr20:57392686         G     A     rs117258475   1.2/1.5               3' end of 3p                no
  **23**   2               miR-130b     chr22:22007634         G     T     rs72631822    0/0                   loop of pre-mir             no
  **24**   1               miR-502      chrX:49779214          C     G     rs781862846   0/NA                  stem, not in mature mir     NA
  **25**   1               miR-450a-1   chrX:133674446^a^      C     T     NA            NA                    stem, not in mature mir     NA

^a^Not confirmed by Sanger sequencing.

^b^Minor allele frequency (MAF) as reported by 1K (1000 Genomes Project phase 3) genotyping data in European population and ESP (NHLBI GO Exome Sequencing Project) in European American population.

^c^Cosegregation analysis performed (yes, variant cosegregates; no, variant does not segregate with PTC).

NA Data not available.

All variants chosen by us in the initial filtering had minor allele frequencies of less than 1% in 1000 Genomes (Phase 1 database). Later we found that of 24 variants chosen, 20 occurred in 1000 Genomes phase 3 of the database (Table [1](#T1){ref-type="table"}) with frequencies up to 2.4%.

We selected 16 variants mostly for being present in a family with at least three affected individuals or for being present in more than one family. We tested them for co-segregation in all available individuals for each family which carried the trait ([Supplementary Table S1](#SD2){ref-type="supplementary-material"}).

Among the 16 tested variants, 6 showed co-segregation with the disease (namely miR-320a, miR-484, mir-15b, mir-135a, mir-181b, and let-7e; see [Supplementary Table S1](#SD2){ref-type="supplementary-material"}). Families carrying variants in miR-320a and let-7e had only two affected members but let-7e presented a high number of individuals with other types of cancers. The family carrying the variant in miR-484 was small and the inheritance pattern was not informative. We decided to use the last four miRNAs in further testing (see family pedigrees in [Supplementary Figure S2](#SD1){ref-type="supplementary-material"}). For the position of the variants in pre-miRNA see [Supplementary Figure S1](#SD1){ref-type="supplementary-material"}.

Testing miRNA expression in patient RNA {#s2_3}
---------------------------------------

To focus on the four miRNAs that showed co-segregation, we requested paraffin blocks of thyroid for all available individuals carrying any of the 4 variants. Unfortunately, only samples carrying miR-135a and miR-181b variants were available. Both variants were shown to be heterozygous in tumor tissue as well as in normal thyroid.

Using TaqMan assay, we tested expression of the two miRNAs (miR-135a and miR-181b) in normal and tumor tissue of both available samples (Table [2](#T2){ref-type="table"}) in order to have each sample work as internal expression control for the micro RNA whose mutation is not present in the sample. The results show greatly reduced expression of miR-135a in the non-tumorous thyroid tissue from the patient with the genomic variant inside miR-135a and correspondingly repressed expression of miR-181b in tissue from the patient with a genomic variant in miR-181b. This result prompted us to test the effects of the four variant miRNAs after transfection of cell lines.

###### TaqMan assay results of miRNA expression in tumor (T) and normal (N) thyroid tissue from 2 patients

                    TaqMan Assay                                  
  ----------------- -------------- ---------- ------------ ------ ------
  miR-135a-1        T              2.70       3.51         1.53   1.41
  variant carrier   N              **0.89**   **1.03**            
  miR-181b-2        T              2.07       1.17         0.87   1.88
  variant carrier   N              **1.84**   **--0.04**          

Testing the effect of the variants on miRNA maturation {#s2_4}
------------------------------------------------------

Since the human miR-135a is expressed from 2 genomic regions (miR-135a-1 in chromosome 3 and miR-135a-2 in chromosome 12), we tested for the expression of primary transcript of the two miRs (pri-miRNAs) in four thyroid tissue tumor/normal pairs in order to see if both pri-miRs are expressed. Both primary transcripts were expressed in the tested samples ([Supplementary Figure S3A](#SD1){ref-type="supplementary-material"}).

The same was done for miR-181b that also is expressed from 2 genomic regions (miR-181b-1 in chromosome 1 and miR-181b-2 in chromosome 9). For this microRNA too, the pri-miRNA is expressed from both genomic regions ([Supplementary Figure S3B](#SD1){ref-type="supplementary-material"}).

In order to test the effect of the variants on processing, we cloned a region of about 500 bp surrounding the pre-miRNA into an expression vector (see Materials and Methods). Constructs containing no insert (empty vector, EV), wild-type genomic region (WT), and genomic region with the variant (MUT) were transfected into the HEK293 cell line and the expression of mature miRNAs was assessed using TaqMan assays. Transfection efficiency of the plasmids in this cell line were greater than 50%. The HEK293 cell line displayed a baseline expression of all four micro RNAs (miR-15b, miR-135a, miR-181b, and let-7e) indicated by C~T~ values around 20 for all miRNAs.

When wild-type construct was transfected, we noted a significant increase in expression for all 4 miRNAs (Figure [1](#F1){ref-type="fig"}). The difference in relative expression between empty vector and wild type construct is due to the amount of baseline endogenous expression of each studied miR and the amount of mature miR expressed by the plasmid.

![Effect of variants on miRNA expressions\
We measured miRNA expression of 4 miRNAs after transcfecting miRNA expressing construct without (WT) or with (MUT) the variant into HEK-293 cell line. Expression of miRNAs was normalized using basal level when transfecting with empty vector (EV). (**A**) miR-15b. (**B**) let-7e. (**C**) miR-135a. (**D**) miR-181b.](oncotarget-08-6475-g001){#F1}

When the mutated constructs were transfected into cell lines, there was no difference relative to WT for miR-15b and miR-135a (Figure [1A](#F1){ref-type="fig"} and [1B](#F1){ref-type="fig"}) while a significantly lower expression was measured for let-7e and miR-181b (Figure [1C](#F1){ref-type="fig"} and [1D](#F1){ref-type="fig"}).

We proceeded to test if the effect of the variants on miR maturation was similar in different cell lines. We transfected the constructs into the BCPAP and COS-7 cell lines and measured the expression of mature miR in order to see the effect of the variants on RNA processing. COS-7 was selected because of its ease of transfection and was thought to be a good cell line to test the basic differences in processing of wild-type constructs vs. mutation carriers.

The variant in the loop of pre-miR-15b showed no effect on processing of mature miRNA in either the BCPAP or the COS-7 cell lines ([Supplementary Figure S4](#SD1){ref-type="supplementary-material"}), while the variant in let-7e affected processing when the construct was transfected into the BCPAP or COS-7 cell lines ([Supplementary Figure S5](#SD1){ref-type="supplementary-material"}). The variant in miR-181b affected expression of mature miR-181b in the HEK293 (Figure [1D](#F1){ref-type="fig"}) and the COS-7 cell lines, but not in BCPAP ([Supplementary Figure S6](#SD1){ref-type="supplementary-material"}). The variant in miR-135a did not affect expression of the micro RNA in BCPAP or COS-7 cell lines ([Supplementary Figure S7](#SD1){ref-type="supplementary-material"}).

Genotyping the candidate variants in PTC case/control populations {#s2_5}
-----------------------------------------------------------------

The two variants that we show to be affecting the maturation of miRNA (namely let-7e and miR-181b) are not listed in any available database (i.e. have not been seen previously) and were only found in one out of 80 families studied in our SOLiD sequencing experiment. After observing the effect of the variants on miRNA processing, we genotyped the two variants in available familial cases, sporadic cases and controls, all from Ohio. The two variants did not occur in any of the 124 familial cases tested, nor in 1170 sporadic cases or 1404 controls.

DISCUSSION {#s3}
==========

As briefly mentioned in the Introduction, the genetic basis of PTC\'s high familiality is not well understood. While ionizing radiation is a well-known risk factor for the development of PTC, especially in children, it is not inconceivable that unknown environmental factors are responsible for the familial aggregation of some PTCs. Searching for hitherto unrecognized genetic and environmental factors therefore is a major challenge. In this study only heritable genetic factors are considered. Although it is becoming clear that genetic risk factors for PTC are likely to be variants with high frequency and low penetrance \[[@R12], [@R13], [@R14]\], here we studied familial PTC that can be expected to display low frequency and high penetrance.

Our study detected only 2 novel variants that affect miRNA expression, and both are rare, occurring in just one family each out of a total of 124 families tested. They also did not occur at all in sporadic PTC (*n* = 1170) or unaffected controls (*n* = 1404). This would appear to suggest that pathogenic genomic variants in pre-miRNA sequences are very rare, but it is worth considering the very real possibility that other variants in the same miRNAs might occur with appreciable frequency; on initial screen the SOLiD sequencing only tested 80 cases. Moreover our filtering criteria may have led us to miss some positive variants. In particular, of the 24 verified variants with the desired genomic location, only 4 remained after testing for co-segregation with disease in the families. Full co-segregation (all affected but no unaffected family members having the variant) is expected in high-penetrance genetics but may not occur with low-penetrance variants. Our conclusion that high penetrance causative variants in miR sequences are rare or very rare must be confirmed by further work.

Of the four candidate variants, the two (in miR-15b and miR-135a) in which the presence of the mutation (variant) did not alter the miRNA\'s expression have also been studied by others, but not implicated in thyroid cancer. The other two candidate miRNAs (let-7e and miR-181b), harboring variants that altered the abundance of their transcript, have already been well characterized in a variety of circumstances. Both are known to be implicated in thyroid cancer \[[@R15]\]. Previously, a higher concentration of let-7e in serum was found to be associated with PTC cases compared to controls \[[@R16]\].

A study published in 2014 found down-regulation of miR-181b and noted a link between microRNA down-regulation and apoptosis through targeting the CYLD gene, however, the study was done in a human papillary thyroid cancer cell line which may not adequately represent the *in vivo* situation \[[@R17]\]. We not that a micro RNA binding site prediction tool such as Diana Tools (<http://diana.imis.athena-innovation.gr/DianaTools/>) lists a large number of predicted targets for every micro RNA. This could well lead to different genes being regulated in different tissues/cell lines.

A recent review by Svoronos et al. \[[@R18]\] addresses the inconsistency in reports of various miRNAs being oncogenic or tumor suppressive. The authors conclude that a miRNA can act as an oncomiR or a tumor suppressor depending on the context. We conclude that the two variants found in this study are likely to be pathogenic, causing strong decrease in the expression of the two miRNAs in normal thyroid. Hypothetically, this could lead to the overexpression of a coding gene that might be an oncogene. In our study we see upregulation of miRs in tumor relative to unaffected tissue from the same organ. This is a common phenomenon \[[@R5]\] and could represent a defense mechanism against carcinogenesis events even in the presence of a pathogenic sequence variant in the pre-miR.

The logical follow-up of this study is searching for target genes that could act as oncogenes and definitive proof will require further experiments.

MATERIALS AND METHODS {#s4}
=====================

Ethics statement {#s4_1}
----------------

The experimental protocols were approved by the Institutional Review Board at the Ohio State University, and all subjects gave written informed consent before participation. Experiments were carried out in accordance with the approved guidelines.

Individuals and clinical samples {#s4_2}
--------------------------------

From a total of \~180 families of PTC in our repository, we chose to study 80 families displaying Mendelian-like inheritance in at least two generations. We defined "familial" as the occurrence of at least two first or second degree relatives with PTC. Family history information, pathology reports confirming the diagnosis of thyroid cancer or other thyroid disease, and blood samples were collected from all consenting affected individuals and key unaffected individuals.

To evaluate candidate predisposing genetic factors found by the next generation sequencing (NGS), we Sanger sequenced blood DNA from an additional 44 familial cases, and genotyped 1170 Ohio sporadic cases and 1404 Ohio controls. The Ohio sporadic cases (*n* = 1170) were individuals with PTC enrolled in the Ohio State University Wexner Medical Center\'s (OSUWMC) endocrine neoplasia repository. Individuals were recruited from a multi-disciplinary thyroid tumor clinic at OSUWMC and all cases were histologically confirmed as PTC. Ohio control samples (*n* = 1404), matched to cases by age, gender and race were provided by the OSUWMC\'s Human Genetics Sample Bank. Recruitment took place in OSUWMC primary care and internal medicine clinics. All cases and controls provided written informed consent, completed a questionnaire that included demographic, medical and family history information, and donated a blood sample. Relevant clinicpathological data for cases were extracted from the electronic medical record.

DNA and RNA extraction {#s4_3}
----------------------

Genomic DNA from blood samples was extracted using the phenol-chloroform method and stored at 4°C. DNA concentration was assessed using NanoDrop 2.0 (ThermoScientific). DNA quality was assessed by loading samples on 1% agarose gel.

Total RNA from cell lines or fresh frozen thyroid tissue was extracted by using TRIzol reagent following the manufacturer\'s procedure and stored at --80°C. Total RNA from paraffin block cores was extracted using NucleoSpin totalRNA FFPE extraction kit (Clontech). RNA concentration was assessed using NanoDrop 2.0 (ThermoScientific). RNA quality was analyzed with an Agilent 2100 Bioanalyzer.

Next-generation sequencing {#s4_4}
--------------------------

An Agilent SureSelect Target Enrichment kit custom designed for the 321 selected genomic regions (with mature microRNA centered and 500 bp flanking on both sides) was used. After enrichment DNAs were sonicated, gel fractionated, and the small DNA fraction (with length not exceeding 40 bp) was subjected to hybridization and ligation with Adaptor Mix. Subsequently, the DNAs were sequenced on the SOLiD platform.

Sanger sequencing {#s4_5}
-----------------

To confirm the presence of miRNA variants and their co-segregation with PTC blood genomic DNA from patients and their available family members was Sanger sequenced. PCR primer sequences are shown in [Supplementary Table S2](#SD3){ref-type="supplementary-material"}. PCR was performed according to a standard protocol using AmpliTaq Gold DNA polymerase (Life Technologies) as follows: 10 min at 94°C; followed by 32 cycles of 15 sec at 94°C, 15 sec at 58°C, and 1 min at 72°C; followed by a final extension of 5 min at 72°C. The PCR amplicons were sequenced using an ABI3730 DNA sequencer (Genomics Shared Resource, OSU).

Cloning of pre-miRNA regions {#s4_6}
----------------------------

Pre-miRNA expression vectors were constructed by amplifying a \~0.5-kb DNA fragment encompassing the pre-miRNA region using the human genomic DNA (heterozygous for the variant of interest) as a template. Using PfuUltra high-fidelity DNA polymerase (Agilent, Santa Clara, CA, USA), PCR reactions were performed with the designed primers ([Supplementary Table S2](#SD3){ref-type="supplementary-material"}) by an Applied Biosystems 9700 Thermal cycler with annealing and elongation temperatures of 58°C and 72°C. The amplified fragments were purified with PCR purification kit (Qiagen) and analyzed on 1.5% agarose gels. The resulting fragments were cloned into a lentiviral vector (pCDH-CMV-EF1-Puro-GFP) using *Xba*I and *Bam*HI or *Nhe*I and *Bam*HI restriction enzymes ([Supplementary Table S2](#SD3){ref-type="supplementary-material"}). The plasmids that contained the wild-type or mutant miRNA genes were identified by Sanger sequencing.

Cell lines {#s4_7}
----------

HEK293 and COS-7 cells were cultured in DMEM medium (Gibco) with 10% fetal bovine serum (Gibco) at 37°C in 5% CO~2~. BCPAP cells were cultured in RPMI medium (Gibco) with 10% fetal bovine serum (Gibco) also at 37°C in 5% CO~2~.

Transient transfection {#s4_8}
----------------------

HEK293 and COS-7 cells were seeded onto 12-well plates (Corning) and transfected at 80% cell confluence with 600 ng of empty vector or vector containing wild-type or mutant insert mixed with Lipofectamine 2000 (Invitrogen; DNA to Lipo at 1:5 ratio). BCPAP cells were seeded onto 12-well plates (Corning) and transfected at 80% cell confluence with 500 ng empty vector (as control) or vector carrying wild-type or mutant miRNA region with a 1:5 ratio of DNA to Lipofectamine 2000 (Invitrogen). Transfections were viewed under a fluorescence microscope (Zeiss) and transfections were confirmed to be at least 50% efficient for HEK293 and COS-7 and \~40% efficient for BCPAP.

Quantitative RT-PCR analysis {#s4_9}
----------------------------

To quantify the expression of pri-miRNAs we used one μg of RNA in a reverse transcription reaction using the High Capacity cDNA Reverse Transcription kit including Random Hexamers (Applied Biosystems). The following TaqMan assays were used: miR-135a-1 (assay ID Hs03303129_pri), miR-135a-2 (assay ID Hs03303133_pri), miR-181b-1 (assay ID Hs03302963_pri), and miR-181b-2 (assay ID Hs03303356_pri). GAPDH was used as internal control.

To measure miRNAs expression, cDNA synthesis was performed as described by the manufacturer using TaqMan microRNA reverse transcription kit (Applied Biosystems) with the specific RT primer provided in each TaqMan assay (Applied Biosystems). The following TaqMan assays were used: miR-135a (Cat \# 000460), miR-15b (Cat \# 000390), miR-181b (Cat \# 001098), and let-7e (Cat \# 002406). RNU6B was used as internal control (Cat \# 001093).

The real time PCR reactions were run on the ABI Prism 7900HT Sequence Detection System (Applied Biosystems). Fast TaqMan assay reaction mix was used and the conditions were as follows: 95°C for 5 min followed by 40 cycles at 95°C, 5 sec and 60°C, 5 sec.

Statistical analysis {#s4_10}
--------------------

Log transformed expression data were analyzed by applying analysis of variance (ANOVA), followed by post hoc Tukey\'s Honest Significant Difference (HSD) test adjusting for multiple comparisons. *P*-values \< 0.05 were considered as statistically significant. Results are presented as Mean ± SEM.

SUPPLEMENTARY MATERIALS FIGURES AND TABLES {#s5}
==========================================

The OSUCCC Genomics Shared Resource performed SOLiD sequencing, SNP genotyping, and Sanger sequencing. Genotyping of the variants in Ohio cases and controls was kindly performed by Dr. Julius Gudmundsson at deCODE Inc, Reykjavik, Iceland. The expression vector pCDH-CMV-EF1-Puro-GFP was kindly provided by Dr. Pierluigi Gasparini. Thank you to Blake Bertani and Jason Bice for their technical support.

**CONFLICTS OF INTEREST**

No conflicts to disclose.

**GRANT SUPPORT**

This work was supported by National Cancer Institute Grants P01CA124570 and P30CA16058.
