Genetic Analyses in a Cohort of Pediatric Patients with Congenital Hypothyroidism Based on Congenital Hypothyroidism Consensus Guideline Open Access

Congenital hypothyroidism (CH) stands as the most prevalent endocrine disorder affecting newborns, potentially leading to significant neurodevelopmental impairments if not promptly treated. Due to this reason, the thyroid screening program is widely used during the neonatal period, as in our country, throughout the world. The incidence of CH is between 1/2,000 and 1/4,000 [1]. CH can be categorized into primary and central (secondary or tertiary). Primary CH is distinguished by low free T4 and T3 levels and elevated thyroid stimulating hormone (TSH) levels, traditionally classified into two groups known as thyroid dysgenesis (TD) and dyshormonogenesis (TDH). TD accounts for 60–65% of primary cases of CH and encompasses various atypical developmental patterns of the thyroid gland, such as athyreosis, ectopy, orthotopic hypoplasia, and hemiagenesis [2, 3]. In a mere 5% of TD instances, identifiable genetic origins are traced to mutations in the TSHR gene [4] or in genes encoding transcription factors pivotal for thyroid gland development, such as TTF1/NKX2-1, PAX8, FOXE1, NKX2-5, and GLIS3; variants of these genes are mainly inherited in an autosomal dominant manner [5, 6]. Dyshormonogenesis constitutes approximately 35–40% of permanent CH (PCH) [7]. TDH arises from malfunctions in the hormone biosynthesis process and typically manifests with a thyroid that is either normal size (a normal in situ thyroid) or enlarged (goiter). TDH is primarily inherited in an autosomal recessive manner and is associated with genetic abnormalities in several genes, including TPO, TG, SLC5A5, SLC26A7, SLC26A4/PDS, IYD/DEHAL1, DUOX2,and DUOXA2 [8].

In addition to TD and TDH, all neonates with CH should be carefully examined for dysmorphic features suggestive of syndromic CH and congenital malformations [9]. PCH can be isolated or part of a syndrome, necessitating early detection of syndromic indicators. Syndromes such as Bamforth-Lazarus (due to biallelic loss-of-function variants in the FOXE1 gene) [10] and brain-lung-thyroid (also known as NKX2-1-related disorders) are linked to specific genetic mutations and come with distinct features like TD and neurological symptoms [11, 12]. Other conditions, including Alagille and Pendred syndromes, exhibit CH with organ malformations or hearing loss [9]. Congenital cardiac and renal defects are also more prevalent in CH cases, underlining the importance of comprehensive early assessment [9].

Despite genetic advances, identifying CH’s genetic causes in most patients with TD remains challenging, suggesting the potential influence of unknown genes or environmental factors. Recently, other new genes have been discovered in the etiology of CH, which could confirm this condition [13, 14]. Therefore, research is recommended in patient groups where the detection of genetic causes is more likely. Congenital Hypothyroidism: A 2020–2021 Consensus Guidelines Update recommends genetic testing for CH if there is a familial occurrence of dysgenesis, central hypothyroidism, or TDH [9]. Additionally, it is crucial to conduct genetic studies on syndromic associations with CH patients to enhance genetic counseling and uncover novel candidate genes that may elucidate the association. We assessed the genetic analysis of our patients with PCH, as recommended by the Congenital Hypothyroidism: A 2020–2021 Consensus Guidelines Update Guidelines, and also discussed the impact of novel variants on these cases.

The study commenced after receiving approval from the local Ethics Committee following the Helsinki Declaration. Genetic analyses of the cases with PCH included in the study were conducted using blood samples collected during visits between 2018 and early 2024 at a single center. Patient selection was updated in 2021 according to the up-to-date CH consensus. Testing was performed on familial TD, TDH, or syndromic CH cases.

It was not feasible to conduct genetic analyses immediately in cases diagnosed with PCH before 2018. Therefore, in some cases, these assessments were performed later than 3 years old, which is the recommended time for re-evaluation of the diagnosis to determine whether PCH or transient CH. However, in some cases, these evaluations could be conducted shortly after 3 years if the case was diagnosed with PCH between 2018 and early 2024.

The CH screening program in Turkey has been implemented since 2006. Before 2006, CH was diagnosed using blood samples taken intravenously. After 2006, most of the diagnosed cases in this study were referred through the newborn screening program. Following the establishment of the screening program, a few cases were diagnosed based on thyroid function tests conducted due to known CH presence in siblings or parents, without waiting for the screening program results, based on demand from the family, or due to the physician’s concern about CH. For these reasons, patients were screened, and if abnormalities in thyroid function tests were found, they were referred to our clinic. The cases referred to us were started on treatment if required according to the CH consensus guidelines and were placed under follow-up, and at the age of 3, a distinction between PCH and TCH was made. The study did not include TCH cases or those with incomplete medical data. Patients were placed in the PCH category if they met any of the following conditions: (i) had TD evident on a scan (if any) or ultrasound (USG) at the time of diagnosis (if any), (ii) required an increase in LT4 dosage over time due to a rise in TSH levels during treatment, or (iii) exhibited low FT4 and/or elevated TSH levels upon re-evaluation after 2–3 years of age following treatment withdrawal (with TSH >15 μU/mL at the 1st month or >10 μU/mL at or after the 3rd month) [15, 16].

In cases of PCH, when thyroid tissue appears normal in situ gland or goiter on thyroid USG, it is categorized as TDH. In these cases, we could not perform a perchlorate test to detect abnormalities of iodide organification. Conversely, if thyroid tissue is not visualized on USG or is observed as hemiagenetic or smaller than normal, it is classified as TD. The patients’ thyroid volumes were determined using the child metrics online calculator program (http://www.ceddcozum.com) based on the three-dimensional measurements obtained from thyroid USG. Thyroid volumes falling below a z score of −2 were classified as hypoplasia, those within the range of −2 to +2 as normal, and volumes exceeding a z score of +2 indicative of goiter [17, 18].

Cases of familial hypothyroidism are defined if PCH is present in a sibling or parent. In cases of sibling presentation where one has been diagnosed with PCH and a mutation has been identified through genetic analysis, even if the other sibling is as young as 2–3 years old, genetic analysis has been conducted for these cases. All cases with TDH were included in the study, while only familial cases were included in the TD group. Additionally, cases with syndromic associations among CH patients were included in the study for genetic analysis.

In evaluating cases with detected mutations and those without, the χ² test was used to compare relationships between categorical variables. For numerical parameters in those groups, the t test was applied when the data showed a normal distribution and nonparametric tests were used when the data did not meet the normality assumption. A p value of less than 0.05 was considered statistically significant.

To investigate the molecular etiology of the diagnosis of CH, genomic DNA was isolated from peripheral blood using a QIAamp DNA Blood Mini QIAcube Kit (Qiagen, Hilden, Germany) according to the manufacturer’s instructions. The patients were evaluated with the next-generation sequencing (NGS) method using a commercial kit including 22 genes (DUOX2, DUOXA2, FOXE1, GNAS, HESX1, IGSF1, NKX2-1, NKX2-5, PAX8, POU1F1, PROP1, SECISBP2, SLC16A2, SLC26A4, SLC5A5, TBL1X, TG, THRA, THRB, TPO, TSHB, and TSHR) known to be associated with hypothyroidism (Blueprint Genetics, Seattle, WA, USA). The NGS process was performed using the Illumina Miseq platform (Illumina, CA, USA), and the raw data were analyzed through the Illumina BaseSpace Variant Interpreter bioinformatics program. The raw data were visualized via an Integrative Genomics Viewer. The following criteria were used when analyzing the variations detected in the NGS data in the index cases: (i) all nonsynonymous, missense, nonsense, frameshift, splice-site, indels, and inframe variants in all protein-coding genes; (ii) synonymous or intronic variants affecting the consensus splice sites; (iii) variants with minor allele frequency <1% in population studies (1000 Genomes, Genome Aggregation Database [gnomAD]). In silico prediction tools (MutationTaster, SIFT, PROVEAN, PolyPhen-2, CADD, Varsome), allele frequencies in population studies (1000 G, gnomAD, ExAC), ClinVar and the Human Gene Mutation Database, and American College of Medical Genetics and Genomics (ACMG) criteria were used to evaluate the pathogenicity of novel variants [19, 20]. The HOPE Web server (http://www3.cmbi.umcn.nl/hope/, accessed on May 12, 2024) was used to analyze the structural effects of the novel variants on the protein sequences.

The study included 61 cases from 45 families meeting the criteria outlined in the Congenital Hypothyroidism: A 2020–2021 Consensus Guidelines Update. The age distribution of the cases ranged from 1.8 to 22 years. 32 of the individuals were female, and 29 were male. Among all families (n = 45), the frequency of consanguineous marriage among all families was found to be 40% (18 families). In TD families, 3 out of 5 (60%) had consanguineous marriages, whereas in TDH families, 15 out of 38 (39.5%) had consanguineous marriages. The majority (>95) of our cases were of Turkish descent. The phenotypic features of the thyroid gland in patients with CH and the distribution of the identified variants are illustrated in Table 1. We did not find statistically significant phenotypic differences between patients with identified genetic variants and those without when considering sex and neonatal/perinatal parameters such as birth weight, history of hospital admission after birth, hypoglycemia, or complications during delivery. Additionally, no statistically significant difference was observed in USG findings used to determine the type of CH and TSH and free T4 levels between these two groups.

Pathogenic variants were found in all cases of syndromic associations with CH patients, in 2 cases of familial dysgenesis, and in 22 cases with a normal gland or goiter. The overall frequency of variant detection is as follows: 61 patients from 45 families were studied. Excluding siblings, the analysis percentages for this study are 17 families with positive variants and 28 families with no detected variants, corresponding to a positive variant rate of 37.7%. Segregation was carried out for all families with positive variants. When examining the rates in subgroups, this was found to be 100% in the group with syndromic associations with CH patients, 20% in the group with familial occurrence of dysgenesis, and 36.6% in the TDH group. Variants in the TPO gene are the most frequently encountered, a situation that was identified in 10 families. Variants followed this in the TSHR gene in 7 families, variants in the DUOX2 gene in 5 families, and two variants in the TG and NKX2-1 genes in 2 families.

While most of the mutations in the case had been previously reported, six variants were identified as novel. According to ACMG classification, three were classified as likely pathogenic and three as variants of uncertain significance (Table 2). They were as follows: c.1600G>A (p.Gly534Ser) and c.2442C>G (p.Cys814Trp) variants in the TPO gene, c.1217T>C (p.Val406Ala) in DUOX2 gene, c.6724C>T (p.Arg2242Cys) and c.6877–2A>C variations in the TG gene, c.713G>A (p.Trp238*) in the NKX2-1 gene were described as novel. Furthermore, in the subset of cases categorized within the familial TD spectrum, characterized by the presence of hypoplastic thyroid tissue, an analysis revealed two siblings harboring triallelic (digenic) mutations – TSHR (homozygous mutation) and TG (heterozygous mutation). When examining thyroid nodules aside from goiter in TDH, only 1 case (a case from F16, No. 25 in Table 2) carrying the TG mutation underwent total thyroidectomy at seven due to multinodular goiter.

As seen in Table 3, in the group of syndromic associations with CH patients, there were 2 cases from 2 families. In these cases, an NKX2-1 mutation was detected. Patient 1 (P1) was diagnosed with hypothyroidism during the newborn period, and treatment with an appropriate dose of l-thyroxine was initiated. She did not develop any lung or neurological issues. The only issue was a brief seizure, lasting 1–2 min, after being in the heat for too long when the patient was 5 years old. Tests did not reveal any underlying problems. The patient was observed without any medication, and no more seizures happened during follow-up. When it comes to patient 2 (P2), the diagnosis of hypothyroidism was made based on heel blood tests. At 16 days old, TSH levels measured in our hospital were 28.64 mIU/L (normal range 0.34–5.60), free T4 was 0.93 ng/dL (normal range 0.61–1.12), and thyroglobulin was over 300 ng/mL (normal range 1.6–60). We retested the patient 1 week later and re-evaluated the need for treatment. The second control value was reported as TSH: 6.51 mIU/L and free T4: 0.95 ng/dL, and the case was placed under observation. No medication for hypothyroidism was administered, only follow-up care. Free T4 levels remained normal during follow-ups without medication, and TSH levels returned to normal. The patient did not suffer any lung infections and had no respiratory problems during the newborn period or any lung problems afterward. The first neurological consultation was due to delayed walking, followed by persistent balance problems, which gradually worsened over time. Subsequently, the patient was monitored with a diagnosis of choreoathetosis, and treatment was initiated.

We used NGS to study 22 genes known to be associated with CH in 61 Turkish patients from 45 families affected by CH. Our patient group consisted of two individuals with syndromic associations with CH, 13 with familial TD, and 46 with TDH. The overall frequency of variant detection was found to be 37.7% (out of 45 families, 17 had a positive mutation), consistent with earlier reports [21‒29]. Those studies have shown that between 33.0% and 65% of CH cases were associated with genetic variants. However, the detection rates reported in the literature vary considerably, likely due to the different patient selection methods, targeting only specific genes for analysis, criteria for determining pathogenicity, and different ethnicities of patients. In addition, even considering the highest detection rate, genetic causes cannot be identified in approximately 35–40% of patients based on those studies. This suggests that novel genetic causes have yet to be identified and must be elucidated. This situation has recently been evidenced by the identification of two genes playing a role in the etiology of CH [6, 7]. Genetic research plays a vital role in discovering new candidate genes and enhancing genetic counseling, especially in syndromic associations, TDH, familial TD, and central hypothyroidism. Thus, it will allow us to improve diagnosis, treatment, and prognosis.

Another factor that may affect the rates mentioned above is determining the type of disease. Indeed, in some cases, it is difficult to assign a disease-type diagnosis, such as TD or TDH, even when complex differential diagnostic tools are used [21]. Additionally, the size of the thyroid gland can be affected by l-thyroxine therapy and the body’s iodine status. Moreover, in USG, the imaging method that is most frequently used is suboptimal for technical reasons. In our study, while TSHR mutations, which were identified as the second most common, were expected to be present in TD cases, there were also TSHR mutation cases in the TDH group. A similar result for TSHR was also demonstrated in a study conducted in our country [29]. This study reported TSHR gene variants in 4 cases in the TD group and 7 cases in the TDH group. While this situation supports the difficulty in defining the disease subgroup, it also implies that the etiology in the nonfamilial TD group might be more diverse. As observed in our study and the study conducted in our country [30], TSHR mutations are likely to be common and may account for more than 5% of TD cases. Additionally, a study provided the prevalence (29%) of TSHR gene mutations in children and adolescents with non-autoimmune subclinical hypothyroidism who were not selected through neonatal screening [31]. This indicates that the rate of TSHR gene mutations in PCH cases could be higher.

Biallelic (homozygous or compound heterozygous) and monoallelic (heterozygous) variants in genes like TSHR, TG, DUOX2, and TPO can be associated with TD or TDH. Monoallelic variants in these genes are linked to milder clinical presentations [32, 33]. In our study, we could not report monoallelic variants except digenic inheritance TG (Het) + TSHR (compound Het) gene mutations detected in two siblings from one family. Although our study investigated 22 genes using NGS, a genotyping method superior to the Sanger technology, studies where CH cases were analyzed for candidate genes using the NGS method reported varying rates of digenic (mutations in two genes) or oligogenic (involves mutations in multiple genes) variants because demonstrating oligogenic inheritance may not always be possible. If oligogenic inheritance can be demonstrated, these cases will likely be more severe than monogenic cases due to a gene dosage effect [21]. It has been claimed that this situation may have arisen from a deep, intronic regulatory region or a large gene deletion in the gene that NGS could not detect [25, 32].

In our study, mutations in the TPO gene were most frequently detected among the cases. Some studies on genetic variants in CH cases have reported that the TPO mutation is the most frequently detected variant, consistent with our findings [34‒36]. In a study of 104 cases (83 Turkish and 21 Pakistani) initially diagnosed with TDH, the TPO mutation was the most common variant in Turkish cases, found in 18 out of 83 [37]. However, in two studies conducted in our country, the frequency of TPO mutation was not most commonly encountered [30, 38]. In one study, the most frequently encountered mutation was in the TSHR gene [30], while in another study, TG, TSHR, and DUOX2 explained the genetic etiology in most of the patients [38]. On the other hand, in these studies, the TDH group’s mutation frequency appears higher than ours. Although our patient selection methodology partially differed, a mutation frequency of 64% was detected in the TDH group in the study [30]. In contrast, in our study, this rate was 36.8% for the same group. However, the overall variant detection frequency was 16.4% in the study [30], while it was 37.7% in our study. This is why we believe that the selection criteria for genetic analysis based on the current guidelines are quite rational, considering both benefits and costs.

We presented 2 patients with NKX2-1 variants; while the clinical presentation of P1 included hypothyroidism, the clinical presentation of P2 included neurological findings such as psychomotor delay, benign hereditary chorea (BHC), and hypotonia. NKX2-1 encodes thyroid transcription factor-1 (TTF-1) that binds to and activates the promoters of thyroid-specific genes such as thyroglobulin, thyroperoxidase, and thyrotropin receptor, playing a crucial role in maintaining the thyroid differentiation phenotype. Additionally, TTF-1 is essential for regulating gene transcription in the brain and lungs. The NKX2-1 gene comprises three exons and encodes a protein of 371 amino acids (https://www.uniprot.org). TTF-1 regulates their transcription by binding to the promoter regions of target genes through the homeodomain region encoded by exon 3 [39‒41]. In P1, a de novo c.964_975del (p.His322_Gln325del) variant was identified in the third exon located after the homeodomain region. While only hypothyroidism was present in P1, no clinical pathological findings related to the central nervous system or the lungs were detected. To the best of our knowledge, this represents the second reported case with thyroid involvement only [39]. In the study by Moya et al. [39], two siblings with clinical findings of hypothyroidism and BHC who were found to have a deletion in the NKX2-1 gene, as well as a mother carrying the same mutation but diagnosed only with mild hypothyroidism, were reported. None of these family members had lung involvement. The frameshift variation identified in this family (c.825del, p.Ala276Argfs*75) is located after the homeodomain in the C-terminal region. It was highlighted in a study comparing mutations between the N-terminal and C-terminal domains that the C-terminal mutant protein retains residual DNA-binding activity but not a proper carboxy-terminal transactivation domain, and the N-terminal mutant protein results in the loss of the DNA-binding domain and is incapable of binding DNA [42]. Our case (P1), considering the mutation located in the C-terminal region [39], suggests that mutations occurring after the homeodomain are likely to be associated with milder phenotypes.

In P2, a novel de novo c.713G>A (p.Trp238*) variation was detected within the third exon, specifically within the homeodomain region. This mutation leads to a premature stop codon, producing a truncated protein. In P2, there were no clinical findings related to hypothyroidism or the lungs, but there was a clinical presentation of BHC, psychomotor delay, and hypotonia. A well-conducted study illustrates the disease distribution, with the following findings: 50% of affected individuals had complete brain-lung-thyroid syndrome, 30% had brain and thyroid involvement only, and 13% had chorea only [40]. NKX2-1 variants can have varying functional effects, even in similar protein regions. Clinical variability is observed among affected families and individuals, including monozygotic twins with the same variant as reported in the literature, suggesting no correlation is evident between the mutated region of the protein and the associated phenotype, highlighting the variable penetrance of the different mutations [39‒43]. However, in another carefully conducted study in 2018, gene deletions, frameshift mutations more frequently present with the complete phenotype, and mutations that generate a stop codon or amino acid substitutions show a weaker association with the full-blown phenotype, which is present only in mutations before or in the homeodomain. Frameshift mutations after the homeodomain that generate a partially aberrant protein can also cause the complete phenotype. Furthermore, when the most recent ClinVar data are reviewed (https://www.ncbi.nlm.nih.gov/clinvar), it is observed that while the variations located after the homeodomain region are predominantly frameshift/nonsense variants, disrupting the reading frame, 36 missense variants located after that region have been reported, all except one of which are classified as variants of uncertain significance, which could result in less pathogenic or likely benign variants. These data show that distinct portions of the protein are endowed with diverse developmental functions. Additionally, this underscores the complexity of NKX2-1 regulation, suggesting that its activity is modulated by various factors, many of which remain unidentified or may vary across different tissues [39‒41].

Thyroid peroxidase (TPO) is a membrane-bound, glycosylated, heme-containing protein that catalyzes iodination and coupling of the hormonogenic tyrosines in thyroglobulin to yield the thyroid hormones T3 and T4 [44]. The TPO protein, composed of 933 amino acids, has a domain between the 149th and 711th amino acids (region encoded by exons 5–12) that is conserved across many species and is responsible for heme peroxidase activity [45]. The TPO gene’s exons 8, 9, and 10 encode the heme-binding domain responsible for the enzyme’s catalytic activity [45]. We detected p.Gly534Ser and p.Cys814Trp variants in the TPO gene as novel variants previously not reported in population studies (1000 Genomes, gnomAD) and ClinVar database. In this study, we identified a biallelic p.Gly534Ser variation (exon 10) located in the heme-binding domain (family 6: P9 and P10), which is found as a highly conserved residue in many different species. This variation results in the loss of the glycine amino acid, which provides flexibility to the protein, and its replacement with serine, a larger amino acid. This situation can lead to conformational changes in the protein, potentially altering its function. The biallelic p.Cys814Trp variation identified in family 17 is located in the EGF-Ca2+ binding domain of the TPO protein. It is known that the 16 different cysteine residues located in the EGF-Ca2+ binding domain (codons: 142, 158, 259, 263, 269, 286, 598, 655, 696, 721, 800, 808, 814, 823, 825, 838) are responsible for the formation of 8 different disulfide bonds within the protein [46]. Previously, it has been reported that missense changes in the cysteine amino acids at codons 800, 808, and 838 lead to conformational changes in the protein or a reduction in enzyme activity [47‒49]. Similarly, the p.Cys814Trp variation identified in this study may also affect protein function.

DUOX2 generates hydrogen peroxide, which is required for thyroid peroxidase/TPO activity. DUOX2 has its own peroxidase activity through its N-terminal peroxidase-like domain [50]. In 1 case, a novel missense variant (p.Val406Ala) was detected in the DUOX2 gene (P11). It was compatible with the patient’s hypothyroidism condition without specifying complete or partial iodide organification defect since the perchlorate discharge test could not be performed to confirm this. This variant (p.Val406Ala) is located in the N-terminal peroxidase-like domain of the protein and is a conserved residue in many species. The mutated residue (p.Val406Ala) is located in this domain that is important for binding other molecules and in contact with residues in a domain that is important for the activity of the protein. The mutation might affect this interaction and disturb signal transfer from the binding domain to the activity domain.

Thyroglobulin (TG) is a substrate for the production of iodinated thyroid hormones T4 and T3. Following TG re-internalization and lysosomal-mediated proteolysis, T3 and T4 are released from the polypeptide backbone, leading to their secretion into the bloodstream [51, 52]. We detected two novel variants (p.Arg2242Cys and c.6877–2A>C) in the TG gene. The p.Arg2242Cys change identified in P25 is located in the acetylcholinesterase homology domain of the TG protein and is found as a highly conserved residue in many different species. This variation (p.Arg2242Cys) results in a mutant residue that, compared to the wild type, has a smaller, neutral, and more hydrophobic structure. The difference in size and hydrophobicity between the amino acids resulting from the variation may affect the correct hydrogen bond formation in the wild-type protein. Therefore, the difference in charge will be able to disturb the ionic interaction made by the original wild-type residue. In family 3, a novel c.6877–2A>C variation has been detected at the splice-acceptor site in intron 39 of the TG gene. This variation could disrupt the splicing of intron 39, potentially leading to the skipping of exon 40 in the mRNA sequence. Typically, splice-site variations can lead to a frameshift in the coding sequence through this mechanism. However, these effects need to be investigated through studies conducted at the RNA and protein levels.

In the current investigation, rare new variations in genes known to be related to CH were discovered, adding to the molecular genetic spectrum. When we compare the overall variant detection frequency, the selection criterion for genetic analysis based on the current guidelines is quite rational, considering the benefits and costs; however, new genes are awaiting discovery. Also, TSHR mutations are likely to be common and may account for more than 5% of TD cases if we include nonfamilial TD.

We would particularly like to thank the patients who agreed to participate in this study.

Subjects were enrolled after the medical Ethics Committee (Ethics Committee of the Erzurum City Hospital, Approval No. 2021/05-91) approved the study protocol and received written informed consent from parents. Written informed consent to participate in the study has been obtained from all adult participants and all underaged participants’ parents/legal guardians/next of kin. We have also obtained written informed consent for the publication of identifying details from all adult participants and the parent/legal guardian/next of kin of all underaged participants in Tables 2 and 3.

The authors declare that they do not have any conflicts of interest.

This research received no specific grant from any public, commercial, or not-for-profit funding agency.

E.K. and A.T. designed the study. E.K., A.S.D., and A.C. recruited and clinically characterized the patients. A.T. and O.Y. performed and analyzed the next-generation sequencing. E.K. and A.T. prepared the draft manuscript. E.K., A.T., and A.C. performed a critical review of the manuscript. All authors contributed to the discussion of results and edited and approved the final manuscript.

All data generated or analyzed during this study are included in this article. Further inquiries can be directed to the corresponding author.