RNA-seq–based miRNA signature as an independent predictor of relapse in pediatric B-cell acute lymphoblastic leukemia

Recent genomic and transcriptomic analyses of pediatric B-cell precursor acute lymphoblastic leukemia (BCP-ALL) have revealed novel genetic subgroups.¹ These novel genetic subgroups and their driver mutations affect prognosis; therefore, treatment stratification based on risk classification has been proposed to improve treatment outcomes.^2,3 However, ∼20% of pediatric patients with BCP-ALL relapse, and the prognosis remains poor; indeed, <50% of the patients who relapse survive long term,⁴ whereas, approximately half of the patients who relapse have a favorable prognosis; these include patients with hyperdiploidy and ETV6::RUNX1.⁵ Although minimal residual disease (MRD), a marker of response to therapy, is the strongest prognostic factor,⁶ many patients who achieve undetectable MRD experience relapse; therefore, current prognostic factors are not thought to be sufficiently stratified.

Micro RNAs (miRNAs) regulate posttranscriptional gene expression by functioning as messenger RNA (mRNA) translation repressors or degraders.⁷ Aberrant expression of miRNAs by tumor cells is a novel biomarker or therapeutic target for various malignancies.⁸ Recent studies have revealed characteristic expression patterns according to genetic subgroups of pediatric BCP-ALL as well as their prognostic significance.^9,10 In addition, novel prognostic factors identified by integrated analysis, including miRNAs and mRNAs, have been reported in several types of malignant tumors¹¹; however, few studies have conducted an integrated analysis focusing on miRNA profiles in a large cohort of pediatric BCP-ALL.

In this study, we aimed to clarify the expression profile of miRNAs in pediatric BCP-ALL by performing an integrative analysis of miRNA sequencing (miRNA-seq) data. We then used this information to investigate their clinical significance and identify prognostic factors on which we could base a more refined risk stratification approach.

The TARGET ALL-P2 cohorts have been published previously, and data were used with permission from the database of Genotypes and Phenotypes (accession number phs000218).¹ The study included 111 primary samples of BCP-ALL (of which 35 patients had paired sample data at relapse) from the TARGET cohort, a high-risk cohort with clinical information, for whom both miRNA-seq and mRNA sequencing (mRNA-seq) data were available (supplemental Tables 1 and 2). Of these 111 cases, 105 had single nucleotide polymorphism (SNP) array data (Affymetrix SNP Array 6.0), and 53 had methylation data (HpaII tiny fragment enrichment by ligation-mediated PCR [HELP]) Assay with Roche NimbleGen). The samples sequenced in this study were derived from 69 primary samples (including 8 paired samples at relapse) provided by multiple institutions in Japan (supplemental Tables 1 and 3). All samples were obtained with written consent from the patient or their parent/guardian under the protocols approved by the institutional review board of each institution. All studies were conducted in accordance with the international ethical guidelines for biomedical research involving human subjects. The chemotherapy regimens used in the TARGET cohort were AALL0232¹² and AALL0331,¹³ and in the Japanese cohort, were based on JACLS¹⁴ or TCCSG¹⁵ protocols.

miRNA-seq was performed in the Japanese cohort (validation cohort). Total RNA was extracted using the miRNeasy Kit (Qiagen). Libraries for miRNA-seq were generated using the NEXTflex Small RNA Sequencing kit version 3 (PerkinElmer) and sequenced using the Illumina NextSeq 550 using a 75 bp single-end read protocol. miRNA-seq analysis was performed using the nf-core/smrnaseq pipeline version 1.0.0.¹⁶ Briefly, all steps consisted of the removal of 3′ adapter sequences by TrimGalore version 0.6.5, mapping to mature and hairpin miRNAs using miRbase version 22.1, and mapping to the GRCh37 human reference genome using Bowtie version 1.3.0. After alignment and trimming, sorted binary alignment map (BAM) files were used for further analysis using edgeR version 3.32 for reads per million normalization, and mirtop version 0.4.23 for raw counts. Normalized count data obtained from the regularized logarithm function of the R package DESeq2 version 1.26.0 were used for clustering analysis. Cluster stability was determined by consensus clustering with 1000 iterations using the R package ConsensusClusterPlus version 1.50.0. The function removeBatchEffect from the R package limma package version 3.42.2 was used to remove batch effects for merged miRNA expression during unsupervised clustering analysis. Differential expression analysis was performed using the DESeq2. After the rlog gene-expression level was calculated using DESeq2, Rtsne in the R package was used to map the samples to a 2-dimensional t-distributed stochastic neighbor embedding (tSNE) plot of the top 500 most variable miRNAs (on the basis of median absolute deviation), and the tSNE perplexity parameter was set to 30.

RNA was isolated from fresh/frozen samples using the NucleoSpin Triprep (Macherey-Nagel). Sequencing libraries were prepared using the NEBNext Ultra II Directional RNA Library Prep Kit for Illumina (New England BioLabs). The library was sequenced using an Illumina HiSeq X platform using a 150-bp paired-end read protocol. mRNA-seq was performed in the Japanese cohort. Sequence alignment and read counting were conducted using the Genomon version 2.6.2 pipeline (http://genomon.readthedocs.io/ja/latest/index.html). Read counting was also normalized to obtain fragments per kilobase million values. Fusion transcripts were detected by the CICERO software (supplemental Table 4).¹⁷ To identify BCR::ABL1–like and ETV6::RUNX1 ALL, 198 genes were identified using the recognition of outliers by sampling ends (ROSE) method,¹⁸ including the top 25 genes in clusters 1 to 8 in the TARGET cohort (supplemental Figure 1). Ward hierarchical clustering method was then used for clustering analysis of the ROSE gene set. Samples in cluster 5 had a signature similar to that of BCR::ABL1 ALL and were, therefore, labeled BCR::ABL1–like, whereas samples in cluster 8 had a signature similar to that of ETV6::RUNX1 ALL and were labeled ETV6::RUNX1–like.

We used mRNA-seq data to examine driver mutations in variant calls, following the Genome Analysis Toolkit (GATK) best practices for calling variants (https://gatk.broadinstitute.org/hc/). Variants were called by the GATK HaplotypeCaller version 4.1 according to the GATK best practices for variant calling in mRNA-seq, as implemented in the CRG-CNAG/CalliNGS-NF pipeline (https://github.com/CRG-CNAG/CalliNGS-NF). Overall, 110 known or putative driver genes and target genes of activation-induced cytidine deaminase in BCP-ALL were filtered (supplemental Table 5).³ The filtered variants were then annotated against the COSMIC database version 92 (a cancer mutation catalog) from the GRCh37 assembly to retain those variants previously identified as cancer mutations.¹⁹ 

Additional details for “Methods” are provided in the supplemental Data.

To characterize the miRNA profile of patients with BCP-ALL having poor outcome, we focused on the TARGET ALL-P2 cohort, which mostly consisted of cases with early bone marrow relapse. We included 111 cases for which miRNA-seq, mRNA-seq data, and clinical information were all available and used them as the discovery cohort (supplemental Tables 1 and 2). First, we performed unsupervised consensus clustering of miRNA expression data obtained from primary samples. The analysis revealed a strong cluster (n = 31) that was not associated with well-known genetic subtypes (cluster 2 in Figure 1A-B; supplemental Figure 2A-B). This cluster included diverse genetic subtypes, such as hyperdiploidy (n = 7) and ETV6::RUNX1 (n = 4) fusions. Differential expression analysis of miRNAs with high average expression across the whole population composing this subgroup demonstrated that they were mainly downregulated (Figure 1C-D; supplemental Table 6). Hereafter, we refer to this subgroup with a distinct miRNA expression profile as the miR-low cluster (MLC). Notable differentially expressed miRNAs in the MLC included tumor suppressive (eg, miR-26b,²⁰ miR-101,²¹ miR-29, miR-142,²² miR-146b,²³ and miR-326²⁴) and oncogenic (eg, miR-130b²⁵ and miR-92a²⁶) miRNAs linked to leukemia. Consecutive 2-step unsupervised consensus clustering of the MLC and non-MLC group identified 4 clusters respectively, which generally clustered together with well-known subtypes (supplemental Figure 2C), further highlighting the distinct features of the aberrant miRNA expression profile identified in the MLC.

Next, we assessed the impact of the miRNA signature of the MLC on the clinical outcome of pediatric patients with BCP-ALL. In the TARGET ALL cohort, 83 patients experienced relapse. We found that the MLC group had significantly shorter event-free survival (EFS) than the non-MLC group (median, 21 vs 33 months, respectively; 95% confidence interval [CI], 11.1-26.8; log-rank P = 3 × 10⁻⁵; Figure 1E). Although not significant, the MLC group also showed a trend toward shorter overall survival (OS) (median, 29 vs 64 months, respectively; 95% CI, 21.2-84.9; log-rank P = .09). In addition, when the analysis was restricted to patients with hyperdiploidy (n = 23), the MLC group (7 of 23 patients) had significantly shorter EFS (median, 28 months; 95% CI, 3.3-39.1; log-rank P = 5 × 10⁻⁴). Although the difference in OS was not statistically significant, there was a trend toward shorter OS in the MLC group of those with hyperdiploidy (log-rank P = .15). Multivariable Cox regression analysis of MLC status and established prognostic factors, including end-induction MRD assessed by multicolor flow cytometry,^27,28 showed that MLC status was an independent predictor of relapse (hazard ratio, 2.59; 95% CI, 1.56-4.3; P < .001; Figure 1F). This was also true in the analysis in which each genetic subtype was included as a covariate (supplemental Figure 3). Taken together, these results suggest that the miRNA signature of MLC could be a novel prognostic factor that predicts relapse better than current prognostic factors, such as MRD.

To validate the poor prognostic miRNA profile of MLC, we evaluated the miRNA profile in the Japanese BCP-ALL cohort. We collected primary samples of 69 patients, 34 of whom relapsed (validation cohort; supplemental Tables 1 and 3). Two-step consecutive unsupervised consensus clustering (supplemental Figure 4A-B) identified a cluster (n = 19) that was independent of known genetic subtypes (cluster 2 in supplemental Figure 4C), as observed for the MLC in the TARGET cohort. To identify patients similar to the MLC in this small Japanese cohort, we applied a statistical regression algorithm based on the least absolute shrinkage selection operator.^19,29 This analysis identified a 19-miRNA signature characterizing the MLC (Figure 2A). Next, the weighted sum of the expression levels of these 19 miRNAs was calculated for each patient in the validation cohort and was expressed as the MLC score. The group identified by the 2-step consecutive unsupervised clustering method described earlier (supplemental Figure 4C) had MLC scores almost exclusively in the highest tertile within the validation cohort; therefore, we referred to this group as patients with MLC-like characteristics (Figure 2B). Compared with the non–MLC-like group, patients with MLC-like characteristicswere more likely to relapse (Figure 2C) and had a significantly shorter EFS (median, 25 months; 95% CI, 16.1-33.0; log-rank P = 1 × 10⁻⁶) and OS (median, 42.6 months; 95% CI, 22.3-78.4; log-rank P = 8 × 10⁻⁴; Figure 2D). Furthermore, evaluation of the area under the curve of the time-dependent receiver-operating characteristic curve, which measures the prognostic ability of the MLC score and miRNAs to predict EFS and OS at 2 years, showed that the MLC score was one of the best predictors of poor prognosis compared with each single miRNA in both the TARGET and Japanese cohorts (Figure 2E; supplemental Tables 7 and 8). In particular, among the elements of the MLC score, high expression of miR-92b-3p and miR-1304-3p were good predictors of poor prognosis (supplemental Figure 5A-B). Taken together, these data confirm the use of the MLC profile for predicting poor prognosis in the validation cohort.

To investigate the relationship between MLC status and previously reported genetic alterations in BCP-ALL, we examined driver mutations in short variants and copy number alterations from the discovery cohort (supplemental Tables 9 and 10). In general, the distribution of these genetic alterations did not differ significantly according to MLC status, although alterations in IKZF1 or TP53 and rearrangement in ZNF384 were significantly mutually exclusive with MLC status (Figure 3A-B; supplemental Figure 6A). Variant call and copy number analysis were also performed using whole-exome sequencing from the validation cohort (supplemental Figure 6B-C; supplemental Tables 11 and 12), but we found no abnormalities characterizing MLC-like cases. Analysis of focal copy number (supplemental Figure 7) and DNA methylation data (supplemental Figure 8) from the TARGET ALL cohort did not identify characteristic abnormalities for MLC. In addition, the poor risk chromosomal profile for hyperdiploidy³⁰ was observed in only 1 case in the MLC (supplemental Table 2). Thus, the MLC miRNA profile does not seem to be associated with any of the well-known high-risk genetic subgroups, driver sequence mutations, copy number alterations, such as IKZF1 or IKZF1^plus, or aberrant DNA methylation, making it a strong candidate for novel risk stratification.

Next, we analyzed paired samples (obtained at diagnosis and relapse) to determine whether the miRNA profile changes at the point of relapse. Unsupervised clustering of 43 paired primary and relapse samples identified 2 clusters (Figure 4A). Cluster 1 contained all 20 primary samples classified as MLC at diagnosis; this was considered to be the group with the MLC profile. Eighteen of these were also categorized into the same cluster at relapse (Figure 4B; supplemental Table 13). By contrast, 9 of 23 patients classified as having a non-MLC profile at diagnosis showed a change in their miRNA profile from cluster 2 at diagnosis to cluster 1 at relapse (Fisher exact test, P = .039). This suggests that patients with non-MLC can acquire an MLC profile at relapse, which is consistent with the fact that MLC scores were significantly higher at relapse than at diagnosis in these 9 patients (Wilcoxon signed-rank test, P = .021; Figure 4C). Furthermore, evaluation of the MLC scores for the paired primary and relapse samples from the 23 non-MLC cases showed that the MLC score increased significantly at relapse (Wilcoxon signed-rank test, P = .014; Figure 4D). These results imply that an MLC profile may confer resistance to treatment and contribute to disease progression and relapse.

To identify genes targeted by deregulated miRNAs in the MLC, we performed an integrated analysis of miRNAs and mRNAs in the TARGET cohort (Figure 5A; supplemental Table 14).¹¹ Gene set enrichment analysis of the upregulated genes in the 2331 pairs showing negative correlation between miRNAs and mRNAs revealed significant enrichment of MYC proto-oncogene, bHLH transcription factor (MYC) target and oxidative phosphorylation genes (Figure 5B) as well as genes related to RNA splicing (Figure 5C). Furthermore, gene set enrichment analysis of genes upregulated in the paired primary and relapse samples from the 9 patients who acquired an MLC profile at relapse also showed significant enrichment of MYC target genes and genes related to oxidative phosphorylation (Figure 5D-E). Collectively, these results suggest that these pathways contribute to the aggressive phenotype of the MLC.

Given the enrichment of RNA splicing terms within the genes upregulated in the MLC, we next focused on splicing alterations in the MLC. We identified that almost all of the alternative-splicing events in MLC cases were reductions in intron retention (IR) (Figure 5F; supplemental Table 15), which was also observed when compared with normal pro–B cells, and in 9 relapsed cases that acquired the MLC-like profile at relapse (supplemental Figure 9A-B). Genes showing decreased IR and increased expression because of a reduction in miRNA expression may be crucial for the MLC; therefore, we looked for genes that met all of the following criteria: (1) genes for which expression increased in MLC compared with that in non-MLC samples (supplemental Table 16); (2) genes for which an miRNA reduction in the MLC may be associated with an increase in mRNA expression (compared with that in non-MLCs); and (3) genes with decreased IR (compared with that in non-MLCs and normal pro–B cells). This resulted in identification of the MTA1 gene, an RNA-binding protein gene that posttranscriptionally activates MYC³¹ (Figure 5G). Notably, miR-30e-5p, a previously identified MTA1 suppressor,³² showed a significant inverse correlation with MTA1 expression in the TARGET cohort (Figure 5H) and may play a role in increasing MYC activity in the MLC. Next, we performed gene-set variation analysis.³³ Based on the B-lineage cluster-specific gene set,³⁴ the gene-expression pattern of MLC leukemic blasts resembled that of pre–B cells; this was not the case for non-MLC patterns (Figure 5I). A previous study shows that MYC-driven transcriptional programs contribute to pre–B-cell transformation and that the pre–B-cell receptor signaling expression pattern could serve as a prognostic marker for high-risk BCP-ALL.^34,35 Similarly for MLC, it is suggested that increased pre–B-cell signaling upon activation of the MYC pathway (caused by miRNA aberrations) may be associated with a poor prognosis.

Finally, to identify the mechanism underlying the downregulation in key miRNA expression in the MLC, we explored abnormalities in the miRNA biogenesis pathway (supplemental Figure 10A). We found that expression of DICER1, also a crucial factor for B-cell development,³⁶ differed significantly between MLC, non-MLC, and normal B-progenitor cells (Kruskal-Wallis test, P = 5.0 × 10⁻⁶); DICER1 was downregulated significantly in the MLC compared with the non-MLC samples (Wilcoxon rank-sum test, Bonferroni-corrected, P = 2.0 × 10⁻⁵; Figure 6A), which was also the case in the patients with MLC-like characteristics of the Japanese cohort (Wilcoxon rank-sum test, P = .0153; Figure 6B). In addition, DICER1 was downregulated significantly at relapse, compared with that in the paired primary sample, in all 9 patients who acquired an MLC profile at relapse (Wilcoxon signed-rank test, P = 7.81 × 10⁻³; Figure 6C). Furthermore, DICER1 expression was significantly lower in patients at relapse than at diagnosis in the other B-ALL cohort³⁷ (supplemental Figure 10B), and patients with low expression of DICER1 were associated with poor prognosis, although there was no significant difference between the intermediate and high expression groups in the TARGET cohort, probably because of the high proportion of relapsed cases (supplemental Figure 10C). Indeed, sequencing of mature miRNAs revealed a marked reduction in the 5p- to 3p-mature miRNA ratio in MLC compared with non-MLC samples (Wilcoxon signed-rank test, P = 6.59 × 10⁻¹⁵; Figure 6D), which is consistent with the finding that DICER1 loss-of-function mutations lead to depletion of 5p miRNAs in other cancers.³⁸ In vitro experiments showed that although DICER1 knockdown of REH cell lines had no effect on cell growth, it did confer resistance to cytarabine, for which oxidative phosphorylation has been reported to be involved in resistance,³⁹ suggesting that DICER1 is crucial for chemoresistance of BCP-ALL (Figure 6E-F). Taken together, these observations suggest that the downregulation of DICER1 is a potential mechanism underlying the miRNA dysregulation that characterizes the MLC. Although reduced expression of DICER1 because of DNA methylation or copy number alterations promotes tumorigenesis in other tumors,⁴⁰ these abnormalities were not observed in MLC (supplemental Figure 10D).

This study describes a novel RNA-seq–based miRNA signature that identifies pediatric patients with BCP-ALL at high risk of relapse and can lead to optimized treatment stratification as a novel risk factor. Unlike mRNA expression profiles that are dependent on existing subtypes,⁴¹ the miRNA profile of the MLC was independent of existing subtypes, driver mutations, and methylation changes, meaning that it cannot be identified in the absence of miRNA-seq data. Several studies report prognostic prediction based on miRNA signatures in pediatric ALL.^42,43 However, these studies were limited to analysis of specific miRNAs by quantitative polymerase chain reaction or microarrays, and did not analyze their association with comprehensive genetic alterations in patients with BCP-ALL. Previous studies of miRNA expression in cases of pediatric ALL are inconsistent with respect to methodology and study design, which reduces reproducibility and prevents adequate meta-analysis; therefore, no definitive prognostic miRNAs have been identified.¹⁰ Among miRNA profiling methods, miRNA-seq is the most accurate, sensitive, and specific, and it can quantify expression over a wide dynamic range.⁴⁴ Thus, comprehensive integrative analysis using miRNA-seq allowed us to identify an aberrant miRNA signature in the MLC.

Predicting relapse is particularly important for standard-risk BCP-ALL. Recently, several novel prognostic factors for ALL have been reported. However, we found that IKZF1^plus45 and a focal 22q11.22 deletion combined with IKZF1 alterations,⁴⁶ both of which are poor prognostic factors, were not associated with the MLC. Other recently identified poor prognostic factors, including the poor risk chromosomal profile for hyperdiploidy³⁰ and the prognostic epigenetic regulator mutations¹ (ie, TBL1XR1 and SETD2), were not associated with the MLC. The miRNA profile observed in the MLC could be a predictor of relapse even in hyperdiploid ALL; indeed, it may be a better predictor of relapse than MRD. More precise identification of patients at high risk for relapse at the time of diagnosis based on miRNA expression would improve the prognosis of pediatric BCP-ALL through risk-stratified treatment, including induction of immunotherapy.⁴⁷ Because evaluation of the miRNA transcriptome is still difficult in clinical practice, it was suggested that high expression of miRNA-92b-3p and miRNA-1304-3p could alternatively predict MLC profile and prognosis, suggesting the potential use of their individual assessment, but further validation is needed.

Our integrative miRNA-mRNA analysis showed that a distinctive reduction in miRNA expression in the MLC may lead to increased expression of MYC target genes and genes related to oxidative phosphorylation, which are associated with a poor prognosis and cytarabine resistance in BCP-ALL.⁴⁸ Given the challenges associated with direct targeting of MYC, exploring the therapeutic potential of drugs targeting the oxidative phosphorylation pathway, such as venetoclax,³⁹ warrants further investigation. Elucidating the functional basis for the association between the MLC miRNA profile and adverse outcomes and the founding genetic lesion underlying DICER1 downregulation is beyond the scope of this study, but further research is warranted.

The MLC was characterized by reduced IR, accompanied by increased expression of RNA processing genes. Aberrant alternative splicing (AS) is relevant to BCP-ALL.⁴⁹ In addition, a recent study of BCP-ALL found that recurrent alterations in the RNA machinery genes¹ as well as RNA splicing alterations⁵⁰ have prognostic implications. IR is a type of AS event that occurs when an intron is retained in the final mature mRNA transcript. Many malignancies are characterized by increased IR; however, we found (unexpectedly) that the MLC had lower levels of IR than the non-MLC and normal B-cell precursors. Decreased IR has been reported in SF3B1-mutated myelodysplastic syndrome⁵¹ and breast cancer⁵², and an inverse relationship between IR and cell proliferation has been observed during B-cell development.⁵³ Reduced IR in the MLC is presumably due to altered expression of spliceosome components through dysregulated miRNAs or to increased MYC activity.⁵⁴ Because mRNA transcripts containing introns are degraded by nonsense-mediated decay of host genes, reduced IR is associated with the upregulation of gene products. In the MLC, reduced IR of the MTA1 gene, an activator of MYC, has been implicated in disease pathogenesis. Overexpression of MTA1 results in spontaneous development of B-cell lymphoma, and PAX5 has been identified as a potential target of MTA1,⁵⁵ which could be associated with increased pre-BCR signaling. Pre-BCR signaling may be a useful therapeutic target because a subset of BCP-ALL depends on it for survival,⁵⁶ and Bruton tyrosine kinase inhibitors, such as ibrutinib are effective.⁵⁷ 

This study has several limitations. First, we did not perform functional analyses of the aberrant miRNAs in MLC. The phenotypic changes conferred by these miRNA alterations in BCP-ALL cells in vitro and in vivo need to be investigated. Second, because this study was a retrospective cohort study with a relatively small sample size and various treatment contexts, caution should be exercised when interpreting the results. Future validation is needed to determine whether the results can be replicated in the modern era when immunotherapy is available for relapsed cases. Because only limited miRNA-seq data for pediatric BCP-ALL are still available, further data accumulation in large prospective cohorts will be required to optimize prognostic models for clinical application. Recently, whole-genome sequencing identified new genetic abnormalities in BCP-ALL,⁵⁸ and future analyses combining transcriptome and whole-genome sequencing are needed to obtain a complete overview of the genomics of miRNAs in BCP-ALL.

In conclusion, our miRNA-mRNA integration analysis comprehensively identified the miRNA profiles associated with pediatric BCP-ALL and clearly identified the novel prognostic miRNA signature in a group of patients at high risk of relapse. Risk stratification based on miRNA expression profiles is expected to improve the treatment outcomes for pediatric patients with BCP-ALL.

The authors thank the patients and their families for their cooperation. The authors also thank K. Kodama for technical assistance.

This work was supported by Japan Society for the Promotion of Science (JSPS) KAKENHI grant numbers JP17H04224, JP18K19467, JP20H00528, and JP21K19405 (J.T.); the Project for Cancer Research and Therapeutic Evolution (grant no. JP19cm0106509h9904), the Project for Promotion of Cancer Research and Therapeutic Evolution (grant no. JP22cm0106xxxh000x and 23ama221505h0002), the Practical Research for Innovative Cancer Control (grant no. JP19ck0106468h0001) from the Japan Agency for Medical Research and Development, the Princess Takamatsu Cancer Research Fund (all to J.T.); the Takeda Hosho Grants for Research Manitoba (J.T.); and Ishizue from Kyoto University Research Administration (J.T.).

Contribution: H.K., H.U., H.H., and J.T. designed the study; M.H., D.H., and T. Imamura provided specimens and collected clinical data; H.K. and K.T. prepared the samples; H.K., N.K., Y.N., and S.O. performed sequencing; H.K., H.U., and T. Isobe. performed bioinformatics analysis; S.S., I.K., and K.U. participated in helpful discussions and commented on the research direction; H.K., H.H., and J.T. wrote the manuscript; and all authors reviewed the manuscript during preparation.

Conflict-of-interest disclosure: The authors declare no conflicting financial interests.

Correspondence: Junko Takita, Department of Pediatrics, Graduate School of Medicine, Kyoto University, 54 Kawahara-cho, Shogoin, Sakyo-ku, Kyoto 606-8507, Japan; email: jtakita@kuhp.kyoto-u.ac.jp; and Hidefumi Hiramatsu, Department of Pediatrics, Graduate School of Medicine, Kyoto University, 54 Kawahara-cho, Shogoin, Sakyo-ku, Kyoto 606-8507, Japan; email: hiramatu@kuhp.kyoto-u.ac.jp.