Unraveling the Genetic Landscape of Langerhans Cell Histiocytosis in Korean Patients: Comprehensive Insights from Mutational Profiles and Clinical Correlations
Article information
Abstract
Purpose
This study aimed to conduct a comprehensive genetic analysis of patients with Langerhans cell histiocytosis (LCH), focusing on the frequency of mitogen-activated protein kinase (MAPK) pathway mutations, detailed mutation profiles of MAPK pathway genes, and their correlation with clinical features and prognosis in Korean LCH patients.
Materials and Methods
We performed targeted next-generation sequencing, capable of capturing exons from 382 cancer-related genes, on genomic DNA extracted from formaldehyde-fixed and paraffin-embedded samples of 45 pathologically confirmed LCH patients.
Results
The majority of patients (91.1%) exhibited single-system disease, with bone being the most common location (84.4%). Initial treatments varied, and no patients died during a median follow-up of 6.8 years. Our genetic assays revealed that all patients had MAPK pathway alterations, including BRAF mutations in 51.2%, MAP2K1 mutations in 42.2%, RAF1 mutations in 4.4%, and a KRAS mutation in 2.2%. These mutations were mutually exclusive. Detailed mutation profiles indicated that among the BRAF mutations, there were 18 point mutations and five in-frame deletions, while most MAP2K1 mutations were in-frame deletions, with only one missense mutation. We detected previously unreported variations of point mutations in BRAF, MAP2K1, KRAS, and the first instance of a RAF1-KLC1 fusion in LCH. MAP2K1 mutations occurred more frequently in older patients, whereas BRAF V600 mutations were commonly associated with unifocal bone disease. Genetic mutations did not correlate with high-risk features or event-free survival.
Conclusion
This study identified mutually exclusive MAPK pathway mutations in every LCH patient through comprehensive genetic analysis, highlighting the importance of inclusive testing in understanding the disease’s genetics.
Introduction
Langerhans cell histiocytosis (LCH) is a complex disorder characterized by the clonal neoplastic proliferation of Langerhans-type cells that express CD1a, langerin, and S100 protein, and exhibit Birbeck granules upon ultrastructural examination [1]. As one of the most prevalent forms of histiocytosis, LCH presents with a broad spectrum of clinical manifestations, ranging from spontaneously regressing single lesions to recurrent reactivation and the potential for permanent long-term disabilities or life-threatening multisystem involvement [2,3].
Molecular genetic studies have identified mutations in mitogen-activated protein kinase (MAPK) pathway genes, including BRAF, MAP2K1, KRAS, and NRAS, as pivotal in LCH pathogenesis [4-8]. The BRAF V600E mutation is common in LCH, indicating the potential for targeted therapies, while BRAF indel and MAP2K1 mutations also contribute to disease development through MAPK pathway activation.
Despite these advances in understanding the molecular genetics of LCH, several unresolved issues and conflicting reports persist. Reported frequencies of BRAF mutations in LCH cases range from 35% to 69% [9-12]. Previous studies have also often relied on targeted sequencing or selective exonal coverage of BRAF mutations, overlooking less common mutations outside of known hotspots. The phenotypic differences and prognostic significance associated with various molecular genetic alterations remain to be fully elucidated. Moreover, the influence of ethnic differences on the prevalence and clinical course of LCH is not yet well understood. Despite recent comprehensive genetic analyses of histiocytosis [12,13], there remains a need for integrated studies that combine clinical, histopathological, and molecular genetic data to provide valuable insights into the complexity of LCH.
This study aimed to address these gaps by investigating the frequency of MAPK mutations, including those in BRAF, MAP2K1, KRAS, NRAS, and other related mutations, as well as their correlation with clinical features and prognosis in Korean patients with LCH. The purpose was to enhance understanding of LCH’s molecular basis in this population and inform improved diagnostic and therapeutic strategies, considering the varied findings in existing literature.
Materials and Methods
1. Diagnosis of LCH
We retrospectively investigated cases of LCH managed at Asan Medical Center (AMC) between 2000 and 2015. We identified patients diagnosed with LCH from the Department of Pathology archives and enrolled patients who consented to somatic mutation analysis and agreed to be registered in the Korean LCH network. From these, we selected records with available paraffin-embedded tissue blocks or unstained slides. Relevant clinical and laboratory data were retrieved from the medical records.
To confirm the diagnosis, routinely prepared hematoxylin and eosin–stained slides for all patients were reviewed. In most cases, we used immunohistochemical methods (CD1a, S-100 protein, or langerin) for further confirmation and the VE-1 antibody to detect cytoplasmic staining indicating the presence of the BRAF V600E mutation in selected cases. The VE-1 antibody has been proven to be highly specific for this mutation [14].
2. DNA extraction and targeted next-generation sequencing
Genomic DNA was extracted from formaldehyde-fixed and paraffin-embedded (FFPE) samples from patients confirmed to have LCH. We obtained 6-μm sections from designated areas of each FFPE tissue block. After de-paraffinization with xylene and ethanol, gDNA was isolated using the NEXprep FFPE Tissue Kit (#NexK-9000, Geneslabs), and quantification was performed using the Qubit dsDNA HS Assay kit (Thermo Fisher Scientific) as previously described [15].
Genetic alterations were assessed through two targeted next-generation sequencing panels at AMC: OncoPanel AMC ver. 1 (OP_AMCv1) and OncoPanel AMC ver. 3 (OP_AMCv3). OP_AMCv1 covered a total of 176 genes, including complete exons from 164 genes, partial intronic regions of 38 genes commonly associated with solid cancer, and an additional set of small single-nucleotide polymorphism (SNP) loci (10,034 bp) for genetic profiling. On the other hand, OP_AMCv3 encompassed approximately 1.2 Mbp and comprised 33,524 probes targeting 382 genes. This panel provided complete exon coverage for 199 genes, included 188 hotspot regions, and targeted partial intronic sequences from eight genes frequently involved in cancer-associated rearrangements. The panel design was based on the Genome Reference Consortium Human Build 37 (GRCh37) as the reference genome and was created using SureDesign (Agilent Technologies). In total, 45 LCH cases (23 with OP_AMCv1 and 28 with OP_AMCv3) were evaluated in the study. We focused profiling on genes commonly included in both target panels and somatic mutations in 29 genes related to the MAPK pathway using the OncoPanel. For our analyses, we considered histiocyte infiltration to be sufficient for oncogenic mutation detection when it exceeded 20%.
3. Bioinformatic analysis
The sequencing data were processed as follows: sequenced reads were aligned to the human reference genome (NCBI build 37) using the Burrows-Wheeler Aligner ver. 0.7.17 with default options. De-multiplexing was carried out with the MarkDuplicates tool from the GATK4 package to remove polymerase chain reaction duplicates. De-duplicated reads were realigned at known indel positions using the GATK4 BaseRecalibrator tool ver. 4.1.3.0. Base quality recalibration was performed using the GATK4 ApplyBQSR tool. Somatic variant calling for single nucleotide variants and short indels in “Only tumor” mode option was conducted with the Mutect2 tool. Germline variants were filtered out from the somatic variant candidates using various databases, including the Single Nucleotide Polymorphism database (dbSNP, build 141, with a frequency threshold of 1%), the Exome Aggregation Consortium database (ExAC; release 0.3.1, with a threshold frequency of 0.001), the Korean Reference Genome (KRG) database, and an in-house panel of normal controls. The final somatic variants were annotated with variant consequences using Variant Effect Predictor (VEP; ver. 86). The annotated variants were then transformed into the mutation annotation format (MAF) using the vcf2maf tool. Furthermore, potential false-positive variants were subjected to manual curation using the Integrative Genomics Viewer (IGV). To evaluate structural variation, copy number variations (CNVs) and rearrangements were assessed using the CNVkit ver. 0.9.6, GISTIC ver. 2.0.23, and BreaKmer ver. 0.0.6 algorithms, respectively. After CNV analysis, the GISTIC algorithms were applied to identify significant focal and arm-level amplifications and deletions in segmented files (.cns) from the CNVkit outputs, with a GISTIC value cut-off set at 0.25 according to the software’s manual. To minimize candidates for germline mutations or false positives related to rearrangement alterations by BreaKmer, we applied an in-house panel of normal controls and manual review.
4. Germline variant study
To identify germline variants, we extracted variants sets with “Germline” or “Panel of Normal” classifications in the FILTER column of the VCF files generated by Mutect2, using only the tumor samples. To identify germline variants, we cross-referenced the data with the Clinvar database (https://github.com/macarthur-lab/clinvar). We compared the 24,863 SNPs—encompassing clinical evidence related to pathogenicity, drug responses, and risk factors—with manually curated information.
5. Statistical analysis
We used the Fisher exact test to investigate the association between mutation status and clinicopathological features such as age, sex, anatomical site, and disease extent. The age at diagnosis was compared between groups using the Mann-Whitney U test. For clinical follow-up, we estimated event-free survival (EFS), which was calculated from the date of diagnosis to the date of relapse, progression, or death from any cause. We assessed the correlation between mutation status and follow-up using Kaplan-Meier analysis and the log-rank test. All p-values were two-sided, and p-values less than 0.05 were considered statistically significant. All statistical analysis was performed using R ver. 4.2.0 (R Foundation for Statistical Computing) and SPSS Statistics for Windows, ver. 23.0 (IBM Corp.).
Results
1. Characteristics of LCH patients
We included 45 LCH patients in our analysis. The median age at the time of diagnosis was 4.0 years, with a wide age range of 0.1 to 55.8 years. Among the patients, 34 (75.6%) were in the pediatric and adolescent age group (younger than 18 years), while 11 (24.4%) were adults. The clinical characteristics are summarized in Table 1.
The most common manifestation of LCH was bone involvement, seen in 38 patients (84.4%), followed by skin involvement in six patients (13.3%) and lung involvement in three patients (6.7%). Of the 45 patients, 41 (91.1%) exhibited single-system disease, and four presented with multisystem disease. The single-system disease group was further classified into 21 patients with unifocal bone involvement, 14 with multifocal bone involvement, four with isolated skin involvement, one with isolated lung involvement, and one with isolated lymph node involvement. The multisystem disease group consisted of two patients with risk organ involvement and two without risk organ involvement.
Initial treatments varied: 31 patients (68.9%) received chemotherapy, nine patients (20.0%) were treated with surgery only, and one patient underwent radiotherapy; four patients were managed with watchful waiting.
With a median follow-up of 6.8 years (range, 0.5 to 26.1 years), no patients died, and 14 patients experienced notable events (10 patients with relapse and four patients with poor responses to initial treatment). Among 10 patients with relapse, three experienced multiple relapse events. Three patients were diagnosed as having diabetes insipidus as a permanent consequence, with one patient diagnosed at the time of LCH diagnosis, and two patients diagnosed during follow-up. No patients in this group developed neurodegenerative disease.
2. Mutation profiles
Biopsy specimens were primarily derived from bone in 35 patients (77.8%), followed by skin (n=4, 8.9%), lung (n=2,4.4%), and other tissue types, such as from the gastrointestinal tract, lymph nodes, pituitary gland, and soft tissue masses (each n=1, 2.2%). OncoPanel assay findings revealed that all the patients had MAPK pathway alterations, including 23 patients (51.2%) with BRAF mutations, 19 (42.2%) with MAP2K1 mutations, two (4.4%) with RAF1 mutations, and one (2.2%) with a KRAS mutation (Fig. 1A and B). Among 23 patients with BRAF mutations, 16 had BRAF V600 mutations, and seven had non-V600E mutations. All detected MAPK pathway mutations were mutually exclusive. Furthermore, no CNVs were detected.

Profiles of genetic mutations detected by targeted next-generation sequencing in Langerhans cell histiocytosis. (A) Distribution of genetic mutations across Langerhans cell histiocytosis (LCH) patient Samples. (B) Proportional representation of mutation types associated with LCH. MS, multisystem; SS, single system; VAF, variant allele frequency. (C) Schematic representation of the RAF1-KLC1 fusion detected in a patient with LCH.
Details of the 23 BRAF mutations were as follows (Fig. 2A): 18 point mutations (15 with V600E mutations, 1 with V600D, 1 with V413M, 1 with P686K) and five indel mutations (2 with P.L485_P490delinsF, 1 with p.N486_P490del, 1 with p.N486_P690del, 1 with P.N486_T491delinsK). All but one mutation at V413M occurred at protein kinase domains. Mutations were mainly clustered at exon 15 (n=16), and exon 12 (n=4), with one large deletion covering exon 12 to exon 17.

Lollipop plots of mutations in key Langerhans cell histiocytosis–associated genes. Frequency and distribution of the identified somatic mutations in BRAF (A) and MAP2K1 (B).
Although point mutations were predominant among the BRAF mutations, most of the MAP2K1 mutations were in-frame deletions (n=18, 94.7%), with one point mutation (Fig. 2B). Among the in-frame deletions, the p.Q58_E62del type was most common, occurring 3 times. Other in-frame deletions included p.F53_Q58delinsL (2 occurrences), p.H100_I103delinsQ (1 occurrence), p.K57_G61del (1 occurrence), p.L101_I103delinsL (2 occurrences), p.L98K_K104delinsM (1 occurrence), p.Q56_G61delinsQ (1 occurrence), p.Q56_G61-delinsR (2 occurrences), p.Q56_V60del (1 occurrence), and p.Q58_E62delinsK (2 occurrences). Additionally, there was a single missense mutation: p.P124S. Five of the MAP2K1 mutations occurred at kinase domains, whereas 14 patients had mutations outside kinase domains. Mutations were clustered at exon 2 (n=14) and exon 3 (n=5).
Regarding the RAF1 mutations, 1 was a point mutation (R143W), and the other was a RAF1-KLC1 fusion (Fig. 1C). The KRAS mutation was a point mutation (G13C).
Regarding the germline mutations, four SNPs were included in our analysis: rs1042522 (TP53), rs12252 (IFITM3), rs180-1133 (MTHFR), and rs751975712 (MPL) (Table 2). Among these, rs1042522 (TP53) had the highest incidence in our LCH cohort, being present in 91% of the patients (n=41). When comparing genotypes and allele frequencies among the KRG panel and the LCH cohort, it was found that all patients with LCH had the CG genotype. In contrast, the CG genotype was present in only 46.5% of the individuals in the KRG panel.
3. Correlation of mutation profiles with clinical phenotypes
Among patients with BRAF point mutations (n=18), 11 had unifocal bone disease, four had multifocal bone disease, two had multisystem disease, and one had single-system skin disease. In patients with BRAF indel mutations (n=5), none had unifocal bone disease, but two had skin disease, one had multifocal bone disease, one had lung disease, and one had multisystem disease without risk organ involvement. In patients with MAP2K1 mutations (n=19), seven had single-system multifocal bone disease, 10 had unifocal bone disease, one had a soft tissue mass, and one had multisystem disease with risk organ involvement. One patient with a KRAS mutation presented with multifocal bone disease. Among the patients with RAF1 mutations, one had multifocal bone disease and one had skin disease.
Among patients with BRAF and MAP2K1 mutations (n=42), distinct age-related patterns were observed in their clinical characteristics. Although the overall age distribution did not differ significantly (p=0.386), subgroup analysis revealed notable trends. In patients younger than 10 years (n=25), BRAF point mutations were predominant, identified in 68.0% of cases (n=15). In contrast, among patients older than 10 years (n=17), MAP2K1 mutations were more common, accounting for 64.7% of this subgroup (n=11) (p=0.037).
In terms of organ involvement, bone involvement according to mutational subtype of BRAF point and indel mutation and MAP2K1 mutation was different statistically (p=0.007). Within the subset of patients with single-system disease, unifocal bone disease was notably more common in those with BRAF point mutations, occurring in 61.1% (11 of 18) of cases. Specifically, 68.8% of patients (11 of 16) with the BRAF V600 mutation exhibited unifocal bone involvement. In contrast, none of the patients with BRAF indel mutations (n=5) had unifocal bone disease, while 52.6% of those with MAP2K1 mutations showed unifocal bone involvement.
Five patients with BRAF indel mutations exhibited distinctive clinical characteristics. The indel mutations were found at exon 12. In two neonatal patients, LCH affected the skin and resolved spontaneously. The other three patients with these mutations were adults. Importantly, none of the patients with BRAF indel mutation exhibited the classic presentation of LCH typically seen in children, characterized by a single bone lesion.
Three patients developed diabetes insipidus in two patients with BRAF V600E mutation, and one with MAP2K1 mutation. EFS was not significantly different between patients with BRAF (66.8%, n=23) and MAP2K1 mutations (58.0%, n=19) (p=0.876) (Fig. 3A). When the analysis was limited to patients with single-system disease, EFS was not significantly different between patients with BRAF (71.8%, n=20) and MAP2K1 (61.3%, n=18) mutations (p=0.804) (Fig. 3B). In the subset of patients with single system unifocal bone disease, EFS was not different between patients with BRAF (58.3%, n=11) and MAP2K1 (80%, n=10) mutations (p=0.269) (Fig. 3C). In patients with single system unifocal bone disease, eight out of 11 patients with BRAF mutation and six out of 10 patients with MAP2K1 mutation received chemotherapy, with no significant difference in EFS (p=0.296).
Discussion
We investigated the clinical, pathological, and comprehensive genetic characteristics of 45 patients with LCH. Our findings revealed that all patients diagnosed with LCH had mutations in the MAPK pathway, with BRAF mutations present in 51.2% and MAP2K1 mutations in 42.2%, alongside a few minor and novel mutations, all mutually exclusive. The study also highlighted intriguing patterns of mutation distribution and organ involvement, such as the predominance of BRAF V600E mutations in association with unifocal bone involvement and the distinctive clinical features associated with BRAF indel mutations. Additionally, no CNVs were detected in our study, aligning with previous findings. This further underscores the importance of extensive genetic analysis in the context of LCH, revealing a universal involvement of the MAPK pathway, regardless of age, organ involvement, or disease presentation, and emphasizing the potential discrepancies and differences in BRAF and MAP2K1 mutation frequencies relative to previous studies.
Previous research has reported varying BRAF mutation rates, ranging from 33% to 69% [9-11]. A recent extensive international study found that 61% of patients had BRAF mutations [12], whereas our study identified a higher prevalence of MAP2K1 mutations than previously reported. This variation might be attributable to several factors. Firstly, our study predominantly consisted of low-risk patients, which could explain the lower frequency of BRAF mutations; prior research has indicated these mutations as more prevalent in high-risk LCH patients [16]. Secondly, ethnic variations may contribute to differences in mutation frequency, with incidence and severity of LCH varying according to ethnicity [17,18]. These disparities can be attributed to underlying molecular pathophysiology, though the distribution of molecular alterations between ethnicities has not been fully elucidated. Finally, the composition of our study group, which included both adult and pediatric patients, could have influenced the mutation frequencies observed [19,20]. This highlights the importance of considering patient demographics and disease risk profiles when evaluating the genetic landscape of LCH, and it suggests that broader, more diverse studies may be necessary to fully explain the mutation spectrum in this disease.
Our study has highlighted the critical importance of comprehensive genetic analysis in the context of LCH. Unlike previous studies that employed stepwise testing approaches, starting with BRAF mutations followed by other MAPK pathway mutations, or that used polymerase chain reaction analysis and targeted sequencing for certain selected genes or hotspots [5,9,21,22], our methods revealed that all patients with LCH harbored mutations in MAPK pathway genes. This comprehensive approach successfully identified mutations in all cases, including those with less common BRAF non-V600E mutations and novel variants. Our findings support the adoption of more inclusive genetic testing protocols in the context of LCH, ensuring thorough mutation detection and potentially informing more accurate diagnosis and tailored treatment strategies. Furthermore, recent advancements in targeted therapies, such as BRAF inhibitors (vemurafenib or dabrafenib) or the combination of dabrafenib and trametinib for relapsed or refractory LCH, emphasize the significance of comprehensive genetic screening [8,23,24]. This approach facilitates precision treatment, tailoring therapeutic options to the genetic profiles of individual patients.
In our study, MAP2K1 and BRAF mutations were mutually exclusive, a finding consistent with prior research [5,6]. MAP2K1 mutations predominantly appeared as in-frame deletions, contrasting with the BRAF point mutations. Consistent with our findings, Chakraborty et al. [6] observed that the majority of MAP2K1 mutations in their study were in-frame deletions and further showed that deletions or point mutations resulted in constitutive ERK1/2 phosphorylation, indicating their active role in LCH pathogenesis.
We identified rare LCH-associated mutations that had not been frequently reported before, thereby broadening our knowledge of the genetic variations associated with the disease. BRAF P686K, V413M, MAP2K1 P124S had not been previously reported in LCH patients. Notably, MAP2K1 P124S is a mutation often found in melanoma [25]. Furthermore, KRAS mutations, seldom reported in LCH due to the predominant focus on BRAF and MAP2K1 mutations in previous studies, were detected [12,13]. Additionally, RAF1 mutations, infrequently observed in LCH, were also identified. Although the MBNL1-RAF1 fusion has been reported previously [26], our study was the first to demonstrate a RAF1-KLC1 fusion in LCH, a mutation previously associated with melanoma [27], suggesting its pathogenic relevance. These findings underscore the necessity of comprehensive genetic analyses in the context of LCH to uncover such rare but potentially impactful mutations, offering new insights into the disease’s molecular pathology and potential treatment targets.
Studies on LCH have yielded inconsistent results regarding the correlation between clinicopathologic features and genetic mutations. While some studies report no significant correlations, others indicate a higher prevalence of BRAF mutations in high-risk patients [12,16,22,28]. Earlier research suggested that BRAF mutations were associated with high-risk features and poor outcomes [16]. A recent large-scale international study with stratified analyses has refined our understanding [12]. This study indicates that specific mutation profiles, such as BRAF V600E, are associated with high-risk features, younger age at diagnosis, multisystem involvement, and central nervous system–risk bone lesions. However, these mutations do not independently predict inferior clinical outcomes when adjusted for disease extent.
In our cohort, composed primarily of low-risk patients, BRAF mutations—particularly V600E—were common in younger children under 10 years old with bone disease but were not necessarily associated with high-risk features. Our findings showed that mutations do not significantly affect outcomes in single-system bone disease, which is consistent with a recent international study. However, due to the small number of high-risk patients in our study, the analysis of specific mutations and their association with high-risk features was limited. Recent efforts to collect large-scale data from combined international cohorts covering various ages and clinical presentations have started to clarify the precise relationship between clinical characteristics and molecular findings in LCH [29]. Additionally, future prospective studies with consecutively registered patients are required to comprehensively assess the significance of genetic mutations in LCH and mitigate the risk of potential selection bias inherent in retrospective studies.
Our study had some limitations, including a relatively small sample size and a predominance of low-risk patients. The diversity of clinical presentations within our cohort also posed challenges for conducting detailed subgroup analyses. However, a key strength of our research lies in the comprehensive genetic evaluations we performed. Through these extensive analyses, we were able to identify several novel or rare mutations, underscoring the critical importance of thorough genetic evaluations in the context of LCH. This approach not only enhances our understanding of the disease’s genetic landscape but also highlights the potential for discovering new diagnostic and therapeutic targets through detailed genetic investigations.
In conclusion, our study corroborates earlier findings that BRAF and MAP2K1 mutations are common and mutually exclusive LCH, while also identifying BRAF non-V600E and other MAPK pathway mutations. Through comprehensive genetic analyses of 45 patients, we demonstrated the universal presence of MAPK pathway mutations in LCH, emphasizing the importance of inclusive genetic testing. Our study also uncovered rare point mutations and reported the first RAF1-KLC1 fusion associated with LCH, contributing to a deeper understanding of the genetic intricacies of the disease. Furthermore, our results underscore the potential for genetic testing to guide precision medicine in LCH. The identification of specific mutations has implications for therapeutic decision-making, particularly in the use of targeted therapies. Comprehensive genetic screening can facilitate the selection of these therapies, enabling tailored treatment strategies that align with each patient’s molecular profile.
Notes
Ethical Statement
The study was conducted after approval by the relevant institutional review board of Asan Medical Center (2021-1833). The requirement for informed consent was waived due to the retrospective nature of the study.
Author Contributions
Conceived and designed the analysis: Koh KN, Chun SM.
Collected the data: Koh KN, Jun HR, Kim JY, Yoon SH, Koh YK, Kang SH, Kim H, Im HJ, Chun SM.
Contributed data or analysis tools: Koh KN, Jun HR, Lee JY, Yoon SH, Koh YK, Kang SH, Kim H, Im HJ, Chun SM.
Performed the analysis: Koh KN, Jun HR, Chun SM.
Wrote the paper: Koh KN, Jun HR, Chun SM.
Conflicts of Interest
Conflict of interest relevant to this article was not reported.
Funding
This study was supported by grants from the Korea Disease Control and Prevention Agency (2019ER690301, 2022ER050200) and the Asan Institute for Life Sciences, Asan Medical Center, Seoul, Korea (2022IF0007).
Acknowledgments
We also extend our gratitude to Lee Changik, Park Yerang, and Lee Sangmin for their generous support for this research.