Multiplex PCR Targeted Amplicon Sequencing (MTA-Seq): Simple, Flexible, and Versatile SNP Genotyping by Highly Multiplexed PCR Amplicon Sequencing

Next-generation sequencing (NGS) technologies have enabled genome re-sequencing for exploring genome-wide polymorphisms among individuals, as well as targeted re-sequencing for the rapid and simultaneous detection of polymorphisms in genes associated with various biological functions. Therefore, a simple and robust method for targeted re-sequencing should facilitate genotyping in a wide range of biological fields. In this study, we developed a simple, custom, targeted re-sequencing method, designated “multiplex PCR targeted amplicon sequencing (MTA-seq),” and applied it to the genotyping of the model grass Brachypodium distachyon. To assess the practical usability of MTA-seq, we applied it to the genotyping of genome-wide single-nucleotide polymorphisms (SNPs) identified in natural accessions (Bd1-1, Bd3-1, Bd21-3, Bd30-1, Koz-1, Koz-3, and Koz-4) by comparing the re-sequencing data with that of reference accession Bd21. Examination of SNP-genotyping accuracy in 443 amplicons from eight parental accessions and an F1 progeny derived by crossing of Bd21 and Bd3-1 revealed that ~95% of the SNPs were correctly called. The assessment suggested that the method provided an efficient framework for accurate and robust SNP genotyping. The method described here enables easy design of custom target SNP-marker panels in various organisms, facilitating a wide range of high-throughput genetic applications, such as genetic mapping, population analysis, molecular breeding, and genomic diagnostics.


INTRODUCTION
Recent advances in high-throughput sequencing technologies have accelerated research and developments in the field of genomics (Metzker, 2010;Schneeberger, 2014). The dramatic increase in sequenced genomes provides substantial benefits to genomics-based approaches. Currently (as of April 12, 2017), the whole-genome data for 755 animal and 247 plant species are publicly available (http://www.ncbi.nlm.nih.gov/genome/browse/). These datasets provide valuable information for developing genetic tools to discover useful genes for applications in molecular breeding and diagnostics.
In forward-genetics approaches, genotyping at a genome, as well as a population scale, is the primary step to elucidate associations between genetic polymorphisms and traits of interest, which is essential to perform genetic mapping in quantitative trait locus (QTL) analyses and genomewide association studies (Schneeberger and Weigel, 2011;Schneeberger, 2014). Single-nucleotide polymorphisms (SNPs) are the most frequently observed mutations and are useful as genetic markers. Because of their co-dominant segregation patterns, SNPs are often more informative than the dominant markers. A previous review showed that pooled next-generation sequencing (NGS) with barcoding is an efficient approach to genotype selected candidate genes for several hundreds of individuals (Marroni et al., 2012). This approach has been applied with excellent results in plants, such as wheat, rice, and poplar (Kharabian-Masouleh et al., 2011;Marroni et al., 2011;Tsai et al., 2011). Recently, several time-and cost-effective methods for SNP genotyping through high-throughput sequencing, such as restriction-site-associated DNA sequencing and genotypingby-sequencing, have been widely used to elucidate genetic diversity (Davey et al., 2011;Onda and Mochida, 2016). However, because these methods do not target specific polymorphisms, the frequency of commonly called SNP sites is expected to decline along with increasing population size (Marchini and Howie, 2010;Rutkoski et al., 2013;Fu, 2014;Fu et al., 2016), and, therefore, they might not be suitable for genotyping panels.
The sequencing of multiplex PCR-based amplicons is an ideal approach to assaying a polymorphic panel of several hundreds of target SNP markers at the population scale and using a wide range of genetic applications, such as phylogenetic analysis, population structure analysis, and QTL mapping (Brachi et al., 2011;Andrews et al., 2016). Currently, commercial kits for amplicon sequencing are available, including the TruSeq Custom Amplicon developed by Illumina (Csernák et al., 2016) and the Ion AmpliSeq Custom DNA Panels marketed by Thermo Fisher Scientific (Millat et al., 2014;Glotov et al., 2015). Moreover, targeted amplicon sequencing effectively allows genotyping of small amounts of DNA derived from samples, such as dried herbarium and fossil specimens (Gugerli et al., 2005;Délye et al., 2013;Beck and Semple, 2015). Therefore, SNP marker panels optimized for multiplexed amplicon sequencing are useful in various research fields, including evolutionary, ecological, agricultural, and population genetics (Savolainen et al., 2013;Tiffin and Ross-Ibarra, 2014;Kumar et al., 2015).
Brachypodium distachyon is an annual grass species that belongs to the Pooideae subfamily and is phylogenetically related to economically important temperate cereals, such as wheat and barley (Draper et al., 2001;Bevan et al., 2010;Mur et al., 2011). B. distachyon is proposed as a model plant for the study of biological phenomena in temperate cereals and biofuel crops due to its tractable features, such as small body size, simple growth requirements, self-fertility, a short life cycle, and a small, diploid genome (Bevan et al., 2010;Girin et al., 2014;Kellogg, 2015). It is expected that the genetic diversity found in natural accessions of B. distachyon is useful for examining genetic polymorphisms associated with traits involved in adaptation to local environments (Filiz et al., 2009;Vogel et al., 2009;Mur et al., 2011;Gordon et al., 2014). To characterize the natural and genotypic diversity in Turkish populations of B. distachyon, Vogel et al. (2009) developed 398 simple sequence repeat markers to determine phylogenetic relations among 187 natural accessions, which indicated considerable genetic diversity. Eight independent inbred lines developed from accessions from Kozluk ("Koz accessions") showed remarkable phenotypic differences, such as the existence or non-existence of hairs on the lemma  and differential disease responses to infection with barley stripe mosaic virus (BSMV) (Cui et al., 2012), although the Koz accessions showed modest genetic diversity and were on the same branch of the phylogenetic tree based on the 187 natural accessions . Additionally, Bsr1, the BSMV-resistance gene in B. distachyon, was mapped using 768 SNP markers designed for the Illumina GoldenGate assay (discontinued by Illumina in 2014) (Huo et al., 2011). A genome-scale SNP panel together with a high-throughput and cost-effective genotyping method without specific platform dependencies would allow gathering polymorphic data from natural accessions and accelerate the mining of genes useful for molecular breeding of related crops, such as wheat and barley.
In this study, we report a useful SNP-genotyping method that we termed "multiplex PCR targeted amplicon sequencing (MTAseq), " which enabled us to genotype >400 SNPs simultaneously (Figure 1). We used the method for high-throughput genotyping of 443 SNPs in eight natural accessions (Bd21, Bd1-1, Bd3-1, Bd21-3, Bd30-1, Koz-1, Koz-3, and Koz-4) of B. distachyon, as well as an F 1 progeny derived by crossing of Bd21 and Bd3-1. We assessed the accuracy and reliability of our method based on genotyping data for homozygous and heterozygous alleles.

MATERIALS AND METHODS
Plant Materials, DNA Purification, and Whole-Genome Re-sequencing Dry seeds of diploid accessions of B. distachyon, Bd21, Bd1-1, Bd3-1, Bd21-3, Bd30-1, Koz-1, Koz-3, and Koz-4, were originally provided by the National Plant Germplasm System of USDA-ARS and were kindly gifted by Dr. David Garvin and Dr. John Vogel. Plants were grown as described previously (Onda et al., 2015). Genomic DNA was extracted from leaves flash-frozen with liquid nitrogen using the DNeasy plant mini kit (Qiagen K.K, Tokyo, Japan). Quality and quantity of the extracted DNA were assessed with a NanoDrop 3300 fluorospectrometer (Thermo Fisher Scientific K. K., Yokohama, Japan). Whole-genome re-sequencing of Bd21, Koz-1, and Koz-4 was conducted with the Ion Proton System (Life Technologies Japan Ltd., Tokyo, Japan). Briefly, DNA libraries for the Ion Proton System were constructed with an Ion Xpress Plus fragment library kit (Life Technologies Japan) and were prepared using the Ion OneTouch 2 System (Life Technologies Japan) for semiconductor sequencing according to manufacturer instructions. The dataset of the Ion Proton System reads has FIGURE 1 | Schematic overview of the MTA-seq workflow used in this study. Primer pairs (443) were designed to generate amplicons harboring one known SNP. Highly multiplexed PCR was conducted with the TaKaRa multiplex PCR amplification kit version 2. Amplicons were end-repaired, adapter-ligated, and nick-repaired with the SPARK DNA sample prep kit for Ion Torrent. To identify individual samples, the Ion Xpress Barcode Adapters 1-96 kit was adopted to add index adapters. The barcoded library was sequenced on an Ion Proton sequencer. SNPs were called by using the "mpileup2cns" command of the VarScan software.

Production of F 1 Progeny by Artificial Crossing
Artificial crossing was performed according to the Garvin Lab method (https://www.ars.usda.gov/ARSUserFiles/1931/ BrachypodiumCrossing.pdf) and the Vogel Lab method (http:// 1ofdmq2n8tc36m6i46scovo2e.wpengine.netdna-cdn.com/wpcontent/uploads/2015/05/Vogel-lab-Crossing-Brachypodium-2-3-2010.pdf). Anthers were pulled out from the female parent by inserting tweezers into the palea to prevent self-fertility, and the lemma was pushed back to close the floret until pollination. After emasculation, all flowers, except the one targeted for crossing, were removed. Next, flowers with very mature anthers were selected. Appropriate anthers with a yellowish-white color were selected, transferred onto a glass slide, and incubated at 37 • C for 5 min to dehisce efficiently. The dehisced anthers were picked up with tweezers and stacked into a female floret to apply the pollen to the stigma. After several days, a developing endosperm was observed if successfully pollinated.

SNP Data Purification Using Whole-Genome Re-sequencing Data of Bd21
To reduce unreliable SNPs in the candidate SNP set, we mapped our Bd21 whole-genome re-sequencing reads (DRR092702) to the Bd21 reference genome (Bdistachyon_192_v1.2) with the Torrent Mapping Alignment Program on the Torrent Server (Thermo Fisher Scientific K. K.). Unreliable SNPs having >5% of different nucleotides between our Bd21 reads and the reference genome were removed from the candidate SNP dataset.

Primer Design
Primer3 (v2.2.3) (Koressaar and Remm, 2007;Untergasser et al., 2012) was used to design primer pairs, each of which may be used to amplify a 150-200-bp amplicon that contains one SNP loci. Some options of Primer3 were modified to improve primer specificity; these were as follows: the optimum size of primer was set to 25 bases, ranging from 18 to 32 bases. The optimum melting temperature was set to 65 • C, ranging from 60 to 72 • C. The range of GC content was set from 45 to 65%. To design primers for the identical genomic regions across the all B. distachyon accessions analyzed, we excluded genomic regions, including SNPs, discovered by at least one of the five SNP callers from the input genomic sequences for the primer design.

Preparation of Amplicon Libraries, Sequencing, and SNP Calling
The workflow of the MTA-seq method is shown in Figure 1. Multiplex PCR was performed with the Multiplex PCR assay kit version 2 (TaKaRa, Kusatsu, Japan) according to manufacturer instruction, with modifications. Briefly, 20 ng of genomic DNA (5 ng/µL) was used as a template for multiplex PCR using a Veriti Thermal Cycler (Thermo Fisher Scientific K.K.) and the 443 primer pairs, with 100 µM of each primer used (Table S2). The thermo-cycling conditions were 94 • C for 1 min, 30 cycles of 94 • C for 30 s and 60 • C for 4 min, and a final extension at 72 • C for 10 min. The multiplex PCR products were purified using the Agencourt AMPure XP system (Beckman Coulter, Brea, CA, USA). For library preparation, we used the SPARK DNA sample prep kit for Ion Torrent (Enzymatics, Beverly, MA, USA) with an Ion Xpress Barcode Adapters 1-96 kit (Thermo Fisher Scientific K.K.), and size selection was performed using E-Gel SizeSelect (Thermo Fisher Scientific K.K.) according to manufacturer instructions. The SPARK DNA sample prep kit is also available for the Illumina platform. Finally, the library was amplified using PrimeSTAR GXL DNA polymerase (TaKaRa) with the custom forward primer (5 ′ -ATCTCATCCCTGCGTGTCTCC-3 ′ ) and reverse primer (5 ′ -TCCGCTTTCCTCTCTATGGGC-3 ′ ). The thermo-cycling conditions were 95 • C for 5 min, and six cycles of denaturation at 95 • C for 15 s, annealing at 58 • C for 15 s, and extension at 70 • C for 1 min. Amplified products were purified with the Agencourt AMPure XP system and quality checked with the High Sensitivity DNA kit (Agilent Technologies, Santa Clara, CA, USA). Emulsion PCR and template-positive ion-sphere particle enrichment were performed using the Ion OneTouch 2 system and Ion OneTouch ES system, respectively, with the Ion PI Hi-Q OT2 200 kit (Life Technologies Japan). The enriched ion-sphere particles were sequenced with the Ion Proton System with Ion PI Chip using the Ion PI Hi-Q Sequencing 200 kit (Life Technologies Japan). Sequence data were downloaded from the Torrent Server to our in-house server and mapped to the Bd21 reference genome with the Torrent Mapping Alignment Program (Thermo Fisher Scientific). SNPs were called by using the "mpileup2cns" command (http://varscan.sourceforge.net/ using-varscan.html#v2.3_mpileup2cns) in the VarScan software (v2.3.7) (Koboldt et al., 2009) using default settings. We verified the results of our SNP genotyping of Bd3-1 and Bd21-3 by wholegenome re-sequencing using the Ion Proton System (the reads has been submitted to the DDBJ Sequence Read Archive under accession numbers DRR124340 and DRR124339, respectively).

SNP Call Simulation by Random Sampling With Decreasing Read Number
The total reads obtained for each accession in this study were randomly sampled to generate 0.01, 0.02, 0.04, 0.07, 0.1, 0.15, 0.2, 0.4, 0.7, 1, 2, and 3 million-read datasets. SNPs were called from the random sampling data with the "mpileup2cns" command of VarScan (v2.3.7) (Koboldt et al., 2009). We performed five independent simulations and calculated the average and standard deviation of the results.

Development of a Set of Even-Spaced SNP Markers in B. distachyon
We identified SNPs between the reference accession Bd21 and eight natural accessions (Bd1-1, Bd3-1, Bd21-3, Bd30-1, BdTR12c, Koz-1, Koz-3, and Koz-4) of B. distachyon, from which we selected a set of even-spaced SNP markers covering the B. distachyon genome as outlined in Figure 2. We performed whole-genome re-sequencing of Koz-1 and Koz-4, which showed remarkable phenotypic differences as compared to Koz-3, with ∼21.3 × and ∼18.3 × coverages, respectively. For the other six accessions, we retrieved whole-genome sequencing datasets from the DDBJ Sequence Read Archive. By mapping the reads of all eight accessions to the Bd21 genome, we identified between 368,211 and 1,406,967 SNPs between Bd21 and the other accessions (Table 1 and Data Sheets 1-8). The number of SNPs was remarkably higher in Bd1-1 and Bd30-1 than in the other accessions (Table 1 and Figure 3). The genetic diversity of Bd1-1 and Bd30-1 was consistent with their phylogenetic relationships reported previously Gordon et al., 2014), whereas Koz-1 and Koz-4 were rather less divergent (Table 1 and Figure 3). As expected, Bd21-3 showed the fewest SNPs, because Bd21-3 and Bd21 were derived from the same source accession (PI 254867) (Garvin et al., 2008;Vogel and Hill, 2008). We then selected 443 biallelic SNPs distributed almost evenly over the five B. distachyon chromosomes. Each of the biallelic SNPs showed a common nucleotide in each of the eight accessions and a polymorphic nucleotide in Bd21. As a result, we selected 124, 94, 101, 78, and 46 SNP markers on chromosomes 1, 2, 3, 4, and 5, respectively. The average distance between adjacent SNPs was ∼600 kb on each chromosome ( Figure S1), indicating that the two closest markers sandwiched 58.5 ± 21.6 genes annotated in the Bd21 genome on average.

Simultaneous Genotyping of Amplicons
To assess the feasibility of our MTA-seq method for SNP genotyping, we designed PCR primers (Table S2 and Figure 2)  targeting each SNP and assessed their accuracy rate in genotyping. For this assessment, we used eight B. distachyon accessions, Bd21, Bd1-1, Bd3-1, Bd21-3, Bd30-1, Koz-1, Koz-3, and Koz-4, as well as an F 1 hybrid derived from crossing of Bd21 and Bd3-1 to assess the genotyping accuracy of MTA-seq for heterozygous alleles. We performed MTA-seq using an Ion Torrent semiconductor sequencer (Ion Proton) and obtained 5.5 million to 7.5 million reads for each accession ( Table 2). Among the total reads obtained, >99.5% of the total reads were mapped onto the Bd21 reference genome, and 2.1 million to 3.6 million reads (39.1-50.9% of obtained reads) were mapped onto the genomic regions of the 443 SNPs in each sample ( Table 2). The accuracy rate of SNP calling by MTA-seq for each accession was between 95.3 and 97.5% (Table 2), which was also in agreement with our whole-genome re-sequencing data of Bd3-1 and Bd21-3, suggesting that MTA-seq is a feasible method for generating amplicons covering >400 SNP markers in one tube, as well as for the simultaneous genotyping of these amplicons using highthroughput sequencers. Over 90% of the SNP markers (400 of 443 SNPs, indicated by asterisks in Table S3) were adequately amplified and accurately genotyped in all of the nine tested samples.

Assessment of Robustness and Scalability of MTA-Seq
To assess MTA-seq robustness and scalability, we conducted a simulation study with the 443 SNP markers designed in this study for natural accessions of B. distachyon and estimated the sequencing-depth threshold for accurate SNP calling. We randomly sampled 10,000 to 3,000,000 reads from the read dataset obtained in this study for each accession, followed by validation of SNP-calling accuracy for increasing numbers of sampled reads. The SNP-calling accuracy rate for an individual accession reached ∼50% when using only 20,000 reads and increased to > 90% when using 400,000 reads, becoming almost saturated when using ≥700,000 (Table S4 and Figure 4). These simulation results verified the high accuracy rate of MTAseq for a small number of sequencing reads and suggested possible scalability for genotyping of a large number of samples simultaneously.

DISCUSSION
In this study, we developed a high-throughput genotyping method, MTA-seq, using NGS technology and examined its reliability using 443 SNP markers designed for B. distachyon. Moreover, we demonstrated that MTA-seq can provide accurate genotyping data for > 400 SNP markers. MTA-seq has the following three main advantages.

Low-Quantity and/or Low-Quality DNA Samples Can be Processed
MTA-seq is developed based on multiplex PCR, and, therefore, is applicable to very limited amounts of input DNA derived, for example, by laser-capture microdissection from formalinfixed paraffin-embedded tissues (Beadling et al., 2013;Singh et al., 2013). It is also applicable to low-quantity and lowquality input DNA samples, such as those derived from dried herbarium and fossil specimens (Gugerli et al., 2005;Délye et al., 2013;Beck and Semple, 2015), which generally contain heavily degraded DNA. MTA-seq requires amplicons of only 150-200 bp for genotyping, circumventing the need for intact DNA. The genotyping of ancient DNA samples by MTA-seq will accelerate the historical study of agricultural species domestication, wildspecies distributions, and plant taxa extinction with spatial and temporal resolution.

MTA-Seq Is Applicable to Large Sample Numbers
For genotyping applications, such as QTL mapping and phylogeographic analysis of wild populations, methodologies that allow high-throughput processing of a large number of samples are advantageous. In this study, we adopted the Ion Xpress Barcode Adapters 1-96 kit, allowing simultaneous processing of 96 samples for library preparation. Moreover, the theoretical number of reads per run that the Ion Proton sequencer can produce is from 60 million to 80 million. Based on our study, 400,000 reads were sufficient for accurate genotyping of > 400 SNPs (90% of the 443 SNP markers) in an individual by MTA-seq (Table S4 and Figure 4). Therefore, the sequencing capacity of the  Table 2). Although we did not investigate the issue in detail, this finding may be attributed to primer dimer formation or mispriming during the highly multiplexed PCR. A previous study reported that fragmented ends of sheared DNA may form incomplete primer sites, potentially resulting in the production of non-target sequences (Campbell et al., 2015). However, FIGURE 4 | Simulation study of genotyping robustness with increasing numbers of randomly sampled reads. The 0.01, 0.02, 0.04, 0.07, 0.1, 0.15, 0.2, 0.4, 0.7, 1, 2, and 3 million NGS reads were randomly sampled from the total NGS reads obtained in this study for each accession. Data are the mean value of five independent random samplings.
non-target-sequence production could be significantly improved by exonuclease I and shrimp alkaline phosphatase treatments prior to library preparation (Campbell et al., 2015). Therefore, if these problems can be solved within the MTA-seq format, this method will allow even higher throughput.

Target SNPs Can be Added Flexibly
One of the strongest advantages of multiplex PCR-based genotyping methods, including MTA-seq, is that they flexibly allow the addition of target SNPs of interest. For example, after QTL identification using genome-wide SNP markers as reported in this study, locus-specific SNP markers accumulated in a QTL region can easily be added to the MTA-seq format to narrow the location of the QTL and identify the genes underlying the QTL. Importantly, MTA-seq can be applied to any SNP in any organism if a reference genome is available. Conversely, a limitation of the MTA-seq method is that the DNA sequence around the target SNPs must be known in order to design primers.

Comparison of MTA-Seq With Other Non-commercial SNP-Genotyping Methods
Recently, other non-commercial SNP-genotyping methods combining NGS technology with multiplex PCR, such as a three-round multiplex PCR-based SNP genotyping method (Chen et al., 2016), and Genotyping-in-Thousands by sequencing (GT-seq) (Campbell et al., 2015), have been reported. Here, we performed a comparison between these two methods and MTA-seq with respect to factors that are important for amplicon sequencing-based genotyping.
(1) With regard to the amount of genomic DNA used as the multiplex PCR template, MTA-seq and the other two methods were found to require comparable amounts of input DNA (less than 50 ng of genomic DNA; e.g., 10-50 ng for GT-seq, 20 ng for the three-round multiplex PCR-based method, and 50 ng for MTA-seq).
(2) With regard to the size of the amplicon for genotyping, GT-seq was performed with less than 100 bp of amplicon, whereas the size of amplicon for three-round multiplex PCR-based method and MTA-seq was 107-160 and 150-200 bp, respectively. Shorter amplicon size may be more suitable for samples in which the genomic DNA is heavily degraded.
(3) The throughput and accuracy rate of MTA-seq SNP calling are comparable to, or slightly superior than, those of the two above mentioned methods. The three-round multiplex PCRbased method allowed genotyping of 37 SNPs in 757 human genomic DNA samples simultaneously, with a SNP-calling accuracy rate of 90.5%, whereas GT-seq allowed genotyping of 192 SNPs in 2,068 steelhead trout genomic DNA samples, with a SNP-calling accuracy rate of 96.4%. Using MTA-seq, we genotyped 443 SNPs in nine B. distachyon genomic DNA samples simultaneously, with 95.3-97.5% accuracy. Although the number of samples capable of being handled can be increased to 384 by using a greater number of adapters (e.g., NEXTflex-HT Barcodes), the relatively small sample number is the weakest point of MTA-seq, as compared with the other two methods.
(4) With respect to the length of primers, primers used in MTA-seq range from 18-31 bases in length, whereas the other two methods used much longer primers (e.g., 50-74 bp for GT-seq, and 37-38 bp for the three-round multiplex PCRbased method). The reason for the longer primer length used by these two methods is the addition of adapters by twoor three-round multiplex PCR, unlike the present MTA-seq method in which adapters are added by ligation. Generally, shorter primers are beneficial as their synthesis is relatively inexpensive.
NGS offers a powerful tool that enables users to obtain several gigabases of nucleotide sequence, with reduced cost and time. Whole-genome re-sequencing enables us to determine genetic diversity at a high resolution in an unbiased manner (James et al., 2013;Michael and VanBuren, 2015). However, several challenges arise when dealing with a large number of samples from organisms with complex and large genomes. In these cases, MTA-seq allows for high-throughput genotyping in a cost-and time-effective manner.

CONCLUSION
The simple, flexible, and versatile SNP-genotyping method reported in this study (MTA-seq) will accelerate forwardgenetics approaches, not only in B. distachyon but also in other species, including economically important temperate cereals, such as wheat and barley. The SNP marker panel for MTAseq of B. distachyon and the genotype information reported in this study should be useful in various research fields, including evolutionary, ecological, agricultural, and population genetics.

AUTHOR CONTRIBUTIONS
YO and KM: Designed the experiments and wrote the manuscript; YO and MS: Performed the experiments; YO, KT, and KI: Analyzed the NGS data; All authors read and approved the final manuscript for publication.

FUNDING
This work was supported by a Grant-in-Aid for Challenging Exploratory Research (grant No. 15K14634) and a Grant-in-Aid for Young Scientists (B) (grant No. 17K15214) from the Japan Society for the Promotion of Science to YO. This work was partially supported by JST CREST (grant No. JPMJCR16O4). This work was also supported by Subsidies to pay for a research support person for research staff with family responsibilities from RIKEN to YO.