A reference-guided TILLING by amplicon-sequencing platform supports forward and reverse genetics in barley.

Barley is a diploid species with a genome smaller than those of other members of the Triticeae tribe, making it an attractive model for genetic studies in Triticeae crops. The recent development of barley genomics has created a need for a high-throughput platform to identify genetically uniform mutants for gene function investigations. In this study, we report an ethyl methanesulfonate (EMS)-mutagenized population consisting of 8525 M3 lines in the barley landrace "Hatiexi" (HTX), which we complement with a high-quality de novo assembly of a reference genome for this genotype. The mutation rate within the population ranged from 1.51 to 4.09 mutations per megabase, depending on the treatment dosage of EMS and the mutation discrimination platform used for genotype analysis. We implemented a three-dimensional DNA pooling strategy combined with multiplexed amplicon sequencing to create a highly efficient and cost-effective TILLING (targeting induced locus lesion in genomes) platform in barley. Mutations were successfully identified from 72 mixed amplicons within a DNA pool containing 64 individual mutants and from 56 mixed amplicons within a pool containing 144 individuals. We discovered abundant allelic mutants for dozens of genes, including the barley Green Revolution contributor gene Brassinosteroid insensitive 1 (BRI1). As a proof of concept, we rapidly determined the causal gene responsible for a chlorotic mutant by following the MutMap strategy, demonstrating the value of this resource to support forward and reverse genetic studies in barley.

Introduction

Barley (Hordeum vulgare L.) ranks fourth in yield production among the major cereal crops worldwide and is an important supply for the animal feed, human food, and brewing industries. Barley is a diploid species within the Triticeae tribe, which includes many important polyploid crops such as wheat (Triticum aestivum), oat (Avena sativa), and triticale (a hybrid between wheat and rye [Secale cereale]) (Mascher et al., 2017; Appels et al., 2018; Ayalew et al., 2021; Hu et al., 2021; Li et al., 2021). With its small genome size (∼5.1 Gb) compared with other Triticeae, as well as abundant genetic stocks (over 485 000 germplasm accessions hosted by over 200 repositories around the world) (Knüpffer, 2009), barley was a model species for the cytological, genetic, and biochemical analysis of crop species during the 20th century. Barley genes have typically been identified through classical map-based cloning designed to determine the genes that influence diverse processes such as domestication (Pourkheirandish et al., 2015), plant growth habit (Fu et al., 2005; Turner et al., 2005), disease resistance (Büschges et al., 1997), and seed dormancy (Sato et al., 2016). These studies have offered broad insights with direct applications to polyploid wheat and have greatly advanced our knowledge and understanding of the relationships among cereal crops (Wang et al., 2014; Pourkheirandish et al., 2017; Abe et al., 2019; Fernández-Calleja et al., 2021).

Mutants with phenotypes that differ from those of isogenic wild-type plants play critical roles in genetic analysis (Lundqvist, 2014). Such mutants have enabled the identification of genes involved in plant morphogenesis and development (Mascher et al., 2014; Jost et al., 2016; Li et al., 2019). Induced mutations, caused by either physical (e.g., gamma-rays) or chemical (e.g., ethyl methanesulfonate [EMS], sodium azide [NaN3], or N-methyl-N-nitrosourea) mutagenesis, increase the range and pool of accessible genetic diversity beyond that associated with natural variation (Waugh et al., 2006). For example, very few natural short-culm alleles (i.e., semi-dwarf 1 [sdw1] and uzu) were identified during the first Green Revolution and introduced into commercial barley cultivars (Dockter et al., 2014). However, hundreds of short-culmed mutants have been developed by chemical mutagenesis and can be found in stock collections (Franckowiak and Lundqvist, 2012). Early work attributed these mutants to variations in different genes, inferring that barley plant height could be optimized by combining preferred alleles at separate loci (Dockter et al., 2014). The well-known malting barley cultivar “Golden Promise,” a mainstay of the Scottish whiskey industry and a popular recipient for transformation, was originally isolated as a mutant with short stature and stiff straw after gamma-ray mutagenesis (Schreiber et al., 2020). Although condensed planting was possible in maize (Zea mays) owing to its erect leaf architecture, no natural germplasm of sorghum (Sorghum bicolor) exhibited a similar phenotype. However, multiple mutants with erect leaves were identified from a population of chemically mutagenized sorghum mutants, thus opening the door for condensed planting in sorghum (Jiao et al., 2016).

Along with the rapid development of high-throughput sequencing methodologies, it has become easier and more affordable to clone the causal gene responsible for a phenotype of interest in a particular mutant. For species with small genomes, such as Arabidopsis (Arabidopsis thaliana) and rice (Oryza sativa), whole-genome re-sequencing of segregating populations or bulked mutants enables the fast identification of a causal gene (Schneeberger et al., 2009; Ashelford et al., 2011; Uchida et al., 2011; Abe et al., 2012). In crop species with large, complex genomes, such as allohexaploid wheat (T. aestivum), direct sequencing of sorted chromosomes in multiple allelic mutants facilitates the rapid isolation of causal genes (Sánchez-Martín et al., 2016; Steuernagel et al., 2016). Likewise, high-throughput sequencing of bulked-segregant pools for barley mutants has revealed genes that control plant architecture and chloroplast development (Mascher et al., 2014; Jost et al., 2016; Li et al., 2019).

A population of induced mutants enables the rapid validation of candidate gene function through reverse genetic approaches (Uauy et al., 2017). Targeting induced local lesions in genomes (TILLING) is a technical platform that can identify induced mutations in heteroduplexes of PCR amplicons from pooled samples (McCallum et al., 2000). The heteroduplexed fragments are then discriminated by denaturing high-pressure liquid chromatography, fluorescence-based gel electrophoresis (e.g., LI-COR), or unlabeled capillary electrophoresis of heteroduplexes digested with the endonuclease CelI (Waugh et al., 2006; Yang et al., 2014). Given the low sensitivity of heteroduplex discrimination, DNA pooling is typically limited to 8–16 individuals (Comai and Henikoff, 2006), thus requiring hundreds to thousands of DNA pools for a population of over 10 000 lines. By contrast, high-throughput sequencing of multiplexed PCR amplicons would enable the discrimination of all DNA molecules within larger pools, even those harboring single-nucleotide mutations present in heterozygous individuals within a heavily pooled sample (Tsai et al., 2011). Heterozygous mutants were unambiguously identified from DNA pools of 64–96 individuals in rice (Tsai et al., 2011), 128 individuals in tomato (Solanum lycopersicum) (Gupta et al., 2017), 16 individuals in barley (Schreiber et al., 2019), and 96 individuals from the amphidiploid rapeseed (Brassica napus) (Sashidhar et al., 2020). With the dramatic decrease in sequencing costs, whole-genome re-sequencing of hundreds of mutant lines is now possible for species such as sorghum, with a 730-Mb genome (Jiao et al., 2016).

Several barley populations of induced mutants have been developed, mostly in malting barley varieties (Caldwell et al., 2004; Talamè et al., 2008; Gottwald et al., 2009; Lababidi et al., 2009; Kurowska et al., 2012; Szurman-Zubrzycka et al., 2018; Schreiber et al., 2019). Two of the wild-type cultivars, “Morex” and “Golden Promise,” have high-quality genomes (Jayakodi et al., 2020). Such TILLING mutants have enabled the functional validation of a list of candidate genes resulting from forward genetics and mapping (Komatsuda et al., 2007; Taketa et al., 2008; Yang et al., 2014; van Esse et al., 2017; Poursarebani et al., 2020). In this study, we developed an EMS-induced mutant population in the special barley landrace “Hatiexi” (HTX) and generated a de novo chromosomal-scale assembly of its genome. We established a robust platform for TILLING by amplicon sequencing (amplicon-seq) involving three-dimensional DNA pooling matrices of 64 or 144 individuals for target amplification, along with high-throughput sequencing of multiplexed amplicons. Following bioinformatics analysis, we identified a series of allelic mutants at 119 gene-specific fragments, highlighting the mutant population as a valuable resource for barley reverse genetic studies. Moreover, we applied the MutMap strategy (Abe et al., 2012) to rapidly isolate the causal gene for a chlorosis mutant, illustrating the use of this resource for future forward genetics.

Results

Establishing an EMS-mutagenized population for the barley landrace “Hatiexi”

The barley landrace HTX is a spring-type, two-rowed hulled barley whose entire body (including leaves, stem knots, awns, lemma, and aleurone) turns black when the plants reach full maturity (Figure 1A). This landrace is characterized by higher crude protein and lower starch and crude fat content in grains compared with varieties that are popularly chosen for the malting industry, animal feed, or food processing (Figure 1B). Genetic diversity analysis of nearly 20 000 domesticated barley accessions (Milner et al., 2019) revealed that HTX is distinct from the cultivated barley of either Europe/Near East or East Asia group (Figure 1C).

Figure 1

Development of an EMS-induced mutagenized population in the barley landrace “Hatiexi” (HTX).

(A) Morphological characteristics of HTX. HTX is a two-rowed, spring-type, hulled barley cultivar with black aleurone and black lemma, whose entire plant body turns black when fully mature. The images of a single seed from left to right show the dorsal and ventral face of the mature seed with hulls maintained or removed.

(B) Contents of crude proteins, neutral detergent fibers (NDF), starch, and crude fat in grains of HTX and other popular barley varieties. Hua 30 (H30), Beiqing 7 (BQ7), and Yangsimai 3 (YSM3) are commercial barley varieties used for malting, food, and feed, respectively. The grains were harvested from three biological replicates in Xinxiang City, China, in 2021. Multiple comparisons were performed using Duncan’s test (p < 0.05).

(C) HTX shows a distinct genetic background from either the Europe/Near East or the East Asia barley groups, based on an analysis of genetic diversity among over 20 000 barley germplasms (Milner et al., 2019).

(D) Development of the EMS-mutagenized population. Genomic DNA was extracted from a single M2 plant for each line, and its derived M3 seeds were harvested.

We established a population of mutants in the HTX background by treating seeds with either 22 mM (0.28% v/v) or 32 mM (0.40% v/v) EMS (Table 1) (Gottwald et al., 2009). The germination and survival rates of M1 seedlings were slightly lower in seed batches exposed to 32 mM versus 22 mM EMS (71.40% versus 79.20%). Similarly, fertility (including semi-sterility) of M1 plants was much lower in seeds treated with the higher EMS concentration (67.64% with 22 mM versus 21.83% with 32 mM), resulting in 6764 and 3492 M2 lines, respectively. Most of the 10 256 M1 plants that produced seeds were partially sterile, with ca. 10 seeds per individual. We sowed all the 10 256 M2 lines; 9190 M2 lines grew, and 8525 M2 lines produced M3 seeds (Figure 1D). We collected leaf samples from a single plant that showed growth vigor in each M2 line (Table 1 and Figure 1D) for DNA extraction and thereby obtained a library with 9190 high-quality M2-DNA samples.

Table 1

Overview statistics of the HTX mutagenized population.

Category	22 mM (0.28%, v/v)		32 mM (0.40%,v/v)
	No.	%	No.	%
Statistics of population development
Mutagenized seeds (M1)	10 000	n.a.	16 000	n.a.
Germinated seeds (M1)	7920	79.20a	11 424	71.40a
Harvested plants (M1)	6764	67.64a	3492	21.83a
Planted lines (M2)	6764	67.64a	3492	21.83a
Developed lines (M2)	6603	66.03a	2587	16.17a
Harvested lines (M3)	6288	62.88a	2237	13.98a
Statistics of M2lines with altered phenotypes
Chlorophyll content	475	7.19b	149	5.76b
Plant height	224	3.39b	196	7.58b
Prostrate/semi-prostrate	1	0.02b	17	0.66b
Leaf morphology	6	0.09b	19	0.73b
Tiller	25	0.38b	12	0.46b
Hypersensitivity response mimicry	0	0	3	0.12b
Waxless	16	0.24b	28	1.08b
Brittle stem	4	0.06b	6	0.23b
Row type	14	0.21b	14	0.54b
Spike architecture	31	0.47b	30	1.16b
Pericarp/lemma color	8	0.12b	6	0.23b
Awn	13	0.20b	11	0.43b
Maturation stage	6	0.09b	7	0.27b
Statistics of M3lines with altered phenotypes
Planted lines	52	n.a.	95	n.a.
Phenotype validated	39	75.00c	75	78.95c

n.a., not applicable.

The percentage was calculated based on the number of mutagenized seeds (M1).

The percentage was calculated based on the number of developed lines (M2).

The percentage was calculated based on the number of planted lines (M3).

We also phenotyped all 9190 M2 lines by visual assessment of 13 traits for plant architecture, morphology, and physiology (Table 1 and Supplemental Figure 1). As expected, we observed an overall higher rate of visible phenotypes among the M2 lines derived from mutagenesis at the higher EMS concentration compared with the lower EMS concentration. In particular, we obtained frequent mutants affected in leaf color (including albino, xantha-type, viridus-type, and striata-type phenotypes) or plant height. We planted 147 M3 lines derived from M2 progenitors with visibly altered phenotypes, and 114 of these lines (ca. 77.55%) showed stable inheritance of their respective mutant phenotypes (Table 1). These validated mutants constitute a valuable germplasm resource for future genetic studies in barley.

Generation of a chromosome-level reference genome

High-quality reference genomes empower genetic studies and crop improvement (Varshney et al., 2021). To facilitate the use of our mutagenized HTX population in barley genetics, we constructed a de novo assembly of the HTX genome using a combination of multiple sequencing strategies (PacBio continuous long-read sequencing, Illumina paired-end sequencing, and Hi-C) (Table 2). First, we used 394.72 Gb of sequence (ca. 77× coverage) derived from 23 559 024 PacBio subreads for genome assembly. We polished the resulting assembly with 327.78 Gb (ca. 65× coverage) of Illumina high-quality short reads and 166 million Hi-C read pairs. The final assembly consisted of 3.83 Gb anchored to the seven barley chromosomes, with 10 359 unassigned contigs (size 223.10 Mb) (Table 2). The lengths of the assembled chromosomes varied from 472.1 Mb to 611.7 Mb (Supplemental Table 1), comparable to the released genome assemblies of barley varieties (Jayakodi et al., 2020). The average size of unanchored contigs was smaller than that of the contigs included in the pseudomolecules (Table 2), as smaller contigs are more difficult to map by Hi-C. To assess the accuracy of the HTX pseudomolecules, we inspected Hi-C intrachromosomal contact matrices and chromosomal alignments between the HTX genome and the reference genome of the barley cultivar “MorexV3” (Mascher et al., 2021). The Hi-C matrices showed the expected diagonal–antidiagonal pattern consistent with the Rabl configuration (Cowan et al., 2001) of interphase nuclei (Figure 2A). The two genomes were highly collinear, with the exception of a segmental inversion at the centromeric region of chromosome 6H (Figure 2B). We undertook a reference-based gene annotation strategy (Jayakodi et al., 2020) to project coding sequences of high-confidence (HC) genes from the latest barley reference genome “MorexV3” (Mascher et al., 2021) to the HTX assembly (Table 2). We annotated 42 010 gene models in the HTX genome. We evaluated the completeness of the HTX genome assembly by calculating its benchmarking universal single-copy orthologs (BUSCO) score against plant datasets. The assembly contained 93.5%  single-copy ortholog groups from the embryophyta dataset, highlighting the high quality of this HTX genome assembly.

Table 2

Summary statistics of the HTX genome assembly.

Genomic feature	HTX
Genome assembly (Mb)	4052.30
Number of contigs	12 664
Contig N50 (Mb)	2.77
Longest contig (Mb)	21.90
Contigs on chromosomes	2305
Size of contigs on chromosomes (Mb)	3829.20
Mean size of assembled contigs (Mb)	1.66
Unassigned contigs	10 359
Size of unassigned contigs (Mb)	223.10
Mean size of unassigned contigs (Mb)	0.02
Number of annotated gene models	42 010
Annotated gene sequence (Mb)	105 (2.59%)
Mean gene length (bp)	2499
Mean coding sequence length (bp)	1022
BUSCOs	1417 complete BUSCOs

Figure 2

De novo assembly of the HTX genome.

(A) Hi-C interaction matrix for each of the seven barley pseudomolecules.

(B) Genome-wide gene-based collinearity between HTX and the barley reference genome “MorexV3.”

Optimizing an amplicon-seq pipeline to identify EMS-induced mutants

To establish a high-throughput TILLING platform to identify the EMS-induced mutants, we used an amplicon-seq pipeline (Tsai et al., 2011) that combined a three-dimensional DNA pooling strategy with high-throughput sequencing of pooled amplicons. Accordingly, we generated two sets of multiplexed DNA samples: one set of 512 M2 individuals distributed across eight plates of 8 × 8 wells (64× DNA pooling) (Figure 3A) and another set of 1728 M2 individuals distributed across 12 plates of 12 × 12 wells (144× DNA pooling). DNA samples were pooled within each plate, set of columns, or set of rows (Figure 3A).

Figure 3

TILLING by amplicon-seq.

(A) A subpanel with 512 M2 individual DNA samples was distributed into eight plates with 8 × 8 wells. All 64 individuals from the same plate (Pn), in the same column of each plate (Cn), or along the same row of each plate (Rn) were pooled.

(B) PCR amplification using each DNA pool as a template and pooling of multiple amplicons. The red-highlighted rectangle on the right gel indicates the fragments purified for sequencing library construction.

(C) Discrimination of mutation sites for the HvBRI1 target fragment. Only G/C to A/T transitions are shown. Normalized depth (≥25% mean depth) was used as the threshold to identify the candidate mutations represented in three distinct libraries.

(D) Decoding of pooled information and identification of the three M2 individuals harboring mutations within the 1054-bp fragment amplified for HvBRI1.

To define the optimal number of amplicons for sequencing, we conducted two independent experiments using the 512-individual subpanel by pooling either 47 or 72 gene-specific amplicons (Figure 3B) and then constructed 24 amplicon-seq libraries for Illumina short-read sequencing (Supplemental Table 2 and Supplemental Table 3). We obtained 1.75 Gb or 2.39 Gb of sequences for the two experiments, yielding an average of 531× or 359× coverage, respectively, of the target amplicons for each M2 individual (Supplemental Table 3). As a previous study (Gupta et al., 2017) indicated that 13.61× to 30.31× coverage of the target amplicons was sufficient to identify a mutation in a diploid genome, we were confident that the level of coverage achieved for the amplicons in this study would allow the detection of most, if not all, induced mutations. We processed all sequencing reads through a standard bioinformatics pipeline and filtered the results for normalized depth (≥25% of mean depths) across libraries (Supplemental Figure 2 and Figure 3C), enabling the identification of 65 mutations (for the 47 gene-specific amplicons) and 136 mutations (for the 72 gene-specific amplicons) in the two experiments (Supplemental Table 3). The estimated average mutation rate was lower for lines derived from 22 mM EMS mutagenesis than from 32 mM EMS (1.51 versus 3.56 mutations/Mb). We then performed a third experiment with the 1728-individual subpanel by pooling equal molar amounts of 56 gene-specific amplicons. We constructed 36 amplicon-seq libraries for Illumina short-read sequencing, which yielded 2.98 Gb of sequence per library, representing 255× coverage of each M2 individual within the 144× DNA pools (Supplemental Table 3). Three hundred nine mutations were identified with a rate of 2.31 mutations/Mb.

To validate the detected mutations, we turned to Sanger sequencing of amplicons covering the target mutations in their relevant M2 individuals, as determined by the decoding of the DNA pools (Figure 3D). We thus confirmed 57 of 63 randomly selected mutations (ca. 90.48%) (Supplemental Table 3), referring to 31 distinct amplicons (Supplemental Table 2 and Supplemental Table 4). Collectively, these results demonstrated the robustness of the amplicon-seq-based TILLING platform for uncovering EMS-induced mutations in barley.

Estimation of mutation rate by multiple approaches

In addition to amplicon-seq, we evaluated the mutation rate in the mutagenized population by two other methods. First, we tested the classical TILLING platform with CelI digestion of heteroduplexes and capillary electrophoresis (Gao et al., 2020) for the barley Nud gene (HORVU.MOREX.r3.7HG0719680.1), the determinant of hulled versus naked grain (Taketa et al., 2008). Accordingly, we amplified a 1099-bp genomic fragment covering the complete coding region in 384 DNA pools from a 12× DNA pooling strategy, representing 2662 and 1946 M2 individuals derived from the 22 mM and 32 mM mutagenized populations, respectively. From 23 pools that showed digestion, we confirmed 11 mutations in Nud by Sanger sequencing, corresponding to a mutation rate of 1.67 mutations/Mb for the 22 mM EMS population and 4 mutations/Mb for the 32 mM EMS population (Supplemental Table 5).

Second, we independently sequenced the entire genome of six individual mutants. We mapped the resulting short reads to the HTX genome assembly with 8.69× to 9.90× coverage per individual (Supplemental Table 6). We discarded heterozygous mutations because of genome complexity (>80% transposable elements) and limited sequencing coverage. We identified between 6880 and 27 968 homozygous mutations per mutant genome, corresponding to an average mutation rate of 3.71 mutations (22 mM EMS population) and 4.09 mutations (32 mM population) per megabase. We detected fewer mutations in annotated gene regions relative to the genome-wide average, with 2.06 mutations (22 mM EMS population) and 3.22 mutations (32 mM EMS population) per megabase. The mutation rates estimated by whole-genome re-sequencing were slightly higher than those obtained by classical digestion or amplicon-seq; however, they were comparable to those obtained by exome sequencing in another barley mutagenesis population (Schreiber et al., 2019).

Therefore, the results of both analyses confirmed our successful establishment of a highly mutagenized barley population and implementation of a high-throughput TILLING platform for its exploitation.

Application in reverse genetics: the barley Green Revolution contributor gene BRI1 as an example

We mapped the 510 mutations identified by amplicon-seq (Supplemental Table 3) onto the well-annotated “MorexV3” reference genome and used SnpEff (Cingolani et al., 2012) to predict their effects. We determined that 67.3% of these mutations were located within exons (Figure 4A). Among these, 2.75% introduced a premature stop codon or affected a splicing site, another 40.6% were non-synonymous substitutions, and the remaining 23.9% were synonymous. We obtained the PROVEAN score for each non-synonymous mutation using a web-based online tool (Choi and Chan, 2015), setting a cutoff of less than or equal to −2.5 for higher-confidence mutations that would be deleterious to protein function. Of 207 non-synonymous mutations, 65 were predicted to be deleterious and worthy of future functional investigation.

Figure 4

Application of the HTX mutagenized population for reverse genetics.

(A) Predicted effects of 510 identified mutations within target amplicons on their encoded proteins. The PROVEAN database (Choi and Chan, 2015) was used to score the potential effects of non-synonymous mutations; absolute scores ≤2.5 were considered to indicate neutral mutations, and scores ≥2.5 indicated deleterious mutations.

(B) Target amplicon of HvBRI1 and the mutations detected from 2240 M2 plants. Mutations shown in gray are synonymous, whereas mutations in red result in either non-synonymous substitution or a prematurely terminated (Q822∗) protein (Supplemental Table 7). Homozygous M4 plants were obtained from the mutation lines M6649 (S1098N) and M6945 (R953K), whereas no homozygous mutants were identified from M9396 (Q822∗).

(C and D) Plant height of the M4 mutant plants and the wild-type HTX.

(E and F) The spike length and (F) grain number per spike in the mutants and wild-type HTX. Multiple comparisons were performed using Duncan’s test (p < 0.05).

As one example, we evaluated the effects of mutations in the barley gene HvBRI1 (HORVU.MOREX.r3.3HG0285210.1) (Chono et al., 2003). HvBRI1 encodes a receptor in the brassinosteroid pathway and plays crucial roles in plant growth and architecture (Li and Chory, 1997; Yamamuro et al., 2000). A natural variant in HvBRI1, uzu1, was identified as causing semi-dwarfing and was incorporated into most cultivars to establish the barley Green Revolution in East Asia (Chono et al., 2003; Dockter et al., 2014). We detected 12 mutations within the 1054-bp fragment amplified for this gene in 2240 M2 individuals (Figure 4B and Supplemental Table 7). Homozygous M4 plants were derived from the mutant lines M6649, harboring the S1098N substitution, and M6945, harboring the R953K substitution. Compared with the wild-type HTX, the homozygous mutants of both lines showed a significantly reduced plant height (Figure 4C and 4D) and spike length (Figure 4E). M6649 also had a lower grain number per spike, whereas M6945 was similar to the wild type for that phenotype (Figure 4F). These mutants will complement existing genetic resources for the development of semi-dwarf barley cultivars, confirming the ability of our platform to identify mutants suitable for further investigation and application.

Application in forward genetics: rapid isolation of the causal gene in a chlorotic mutant

We identified the chlorotic mutant M4009 among the M2 plants and found that it was characterized by a yellow color throughout its life cycle (Figure 5A). This mutant showed lower total chlorophyll content and photosynthetic efficiency relative to the wild type (Figure 5B–5D). When we backcrossed M4 mutant plants to HTX to generate a segregating population, the resulting F1 hybrids were all green in appearance, and the chlorotic phenotype segregated in a 1:3 ratio in the F2 population, indicating that a single recessive mutation was responsible for the chlorotic phenotype (Supplemental Figure 3A and 3B). We followed the MutMap strategy (Abe et al., 2012) to map the causal mutation by sequencing genomic DNA extracted from two bulk pools, each representing 60 individual segregants of either green or yellow color from the F2 population, to about 40× coverage each. We observed tight linkage between the mutation of interest and a region surrounding the centromere of chromosome 1H spanning nearly 294 Mb, from 30 Mb to 323 Mb based on the “MorexV3” reference genome (Figure 5E). Of 376 mutations within this interval, only four resided within coding sequences (Supplemental Table 8). One mutation in Protochlorophyllide oxidoreductase B (HvPORB; HORVU.MOREX.r3.1HG0047060.1) resulted in an HvPORB protein that lacked the last five amino acids (Figure 5E). We genotyped 396 F2 individuals with a specific Kompetitive Allele-Specific PCR (KASP) marker designed based on the above HvPORB mutation, and we found perfect co-segregation between yellow F2 seedlings and the mutation (Supplemental Figure 3B and Supplemental Table 9). PORB participates in the chlorophyll biosynthetic pathway in Arabidopsis and rice (Kim et al., 2005; Sakuraba et al., 2013), and a frameshift in PORB results in a chlorotic phenotype in the radiation-induced barley mutant “Nanchong Yellow Barley” (Liu et al., 2008; Yuan et al., 2010). To confirm the identity of the gene affected by the M4009 mutation, we performed virus-induced gene silencing of HvPORB by targeting the 5′ end of the HvPORB transcript in wild-type HTX plants. Newly emerging leaves became chlorotic (Figure 5F), confirming that the modification of HvPORB can result in the chlorotic phenotype, as observed in the M4009 mutant. Together, these results show that following the MutMap strategy, which combines a mutant-derived backcrossing population with high-throughput sequencing, enables effective and fast gene cloning in a complex genome such as barley.

Figure 5

Cloning of the causal gene for chlorotic mutant M4009.

(A) Identification of mutant M4009 in the field-grown M2 population, with its yellowing performance throughout its life cycle.

(B) Mutant M4009 shows dramatically lower chlorophyll fluorescence and maximum quantum efficiency of photosystem II (Fv/Fm).

(C) Total chlorophyll levels (Chla + Chlb) are much lower in the mutant under 1000 μmol·m−2·s−1 high light conditions.

(D) Mutant M4009 exhibits a lower net photosynthetic rate (Pn) relative to the wild-type HTX under 1000 μmol·m−2·s−1 high light conditions.

(E) Identification of the candidate gene responsible for the chlorophyll phenotype in mutant M4009. A physical interval was defined on chromosome 1H, associated with the phenotypic segregation of green and yellow F2 plants. The phenotypic alteration was attributed to a nucleotide substitution resulting in a truncated protein of barley Protochlorophyllide oxidoreductase B (HvPORB; HORVU.MOREX.r3.1HG0047060).

(F) Reducing HvPORB transcript levels by BSMV-induced gene silencing induces leaf chlorosis in wild-type HTX, confirming the mutation at HvPORB as the cause of the chlorotic phenotype in mutant M4009.

Discussion

Cultivated barley was a founder crop of Neolithic agriculture (Weiss and Zohary, 2011); the natural germplasms currently preserved across various repositories maintain vast phenotypic and genotypic diversity (Schulte et al., 2009; Druka et al., 2011; Muñoz-Amatriaín et al., 2014; Russell et al., 2016). However, these accessions of diverse genetic backgrounds complicate the process of associating a given phenotype with one of the many polymorphisms present in each background. By contrast, induced mutants in a given background provide ideal materials with which to elucidate the function of a particular gene (Komatsuda et al., 2007; Taketa et al., 2008; Yang et al., 2014; van Esse et al., 2017; Poursarebani et al., 2020). In this study, we developed an EMS-induced mutant population based on the barley lodging-sensitive landrace HTX, with hulled grains and moderate seed dormancy (Yan et al., 2021). Although HTX was collected from northeast China, phylogenetic estimation places it to a branch distinct from either the Europe/Near East or the East Asia. HTX grains have extraordinarily high protein contents but low levels of starch and crude fat, qualities that may fit the needs of either livestock feed (Högy et al., 2013; Kazemi-Bonchenari et al., 2020) or human food processing (Xie et al., 2021). HTX has a black lemma and aleurone, and the entire plant turns black upon full maturity. The black pigmentation of grains is considered to increase their antioxidant properties (Glagoleva et al., 2017; Long et al., 2018). This mutagenesis population may be useful for investigating malting quality because of its high protein content. Previously established mutant populations are based mainly on malting barley varieties (Caldwell et al., 2004; Talamè et al., 2008; Gottwald et al., 2009; Lababidi et al., 2009; Kurowska et al., 2012; Szurman-Zubrzycka et al., 2018; Schreiber et al., 2019), and these populations are useful and accessible. This HTX mutant population can support theoretical and applied genetic studies to analyze grain quality and meet end-use requirements for livestock feed or food processing.

A highly efficient platform is desirable for detecting mutations within a large mutagenized population comprising thousands of lines. The classical workflow based on CelI digestion and capillary electrophoresis limits pooling to 8–16 individual DNA samples (Caldwell et al., 2004; Till et al., 2007; Gao et al., 2020) and is time and labor intensive when screening a large population (i.e., 8525 lines in this study). High-throughput sequencing enables the screening of multiple genes by amplicon-seq of a pooled library (Tsai et al., 2011) or whole-genome re-sequencing of individuals or pooled samples (Ashelford et al., 2011; Abe et al., 2012; Jiao et al., 2016). We tested all three methods here and found that high-throughput sequencing of pooled PCR amplicons from three-dimensional DNA pools provided the optimal balance between commercial cost and discrimination power for mutations. It was possible to condense the HTX mutant population (8525 M2 lines) into 180 DNA pools, each containing 144 M2 individuals represented three times in each of the pooling dimensions. Rather than sequencing single amplicons, we pooled 60 amplicons with an average length of 1.5 kb; producing 2 Gb of sequence for each of the 180 resulting libraries (360 Gb in total) would result in 154× coverage of the target fragments in each M2 individual, which is sufficient to identify a mutation in a diploid organism (Gupta et al., 2017). We validated 90.48% of the detected mutations by Sanger sequencing of target amplicons in corresponding M2 individuals, demonstrating the reliability of this amplicon-seq-based TILLING workflow. Targeted PCR amplification using gene-specific primers produces a unique fragment, thereby overcoming issues typically associated with polyploid genomes or gene family paralogs (Wang et al., 2014; Appels et al., 2018; Kan et al., 2022) in which highly homologous short reads map to several genomic targets and thus interfere with mutation calling (Krasileva et al., 2017; Sashidhar et al., 2020). Individual mutants have been successfully identified via amplicon-seq in the amphidiploid rapeseed (B. napus) (Sashidhar et al., 2020).

Based on our results from amplicon-seq, we estimated the mutation rates within the two subpopulations derived from 22 mM and 32 mM EMS treatment to be 1.5 mutations/Mb and 3.5 mutations/Mb, respectively. We noticed a slight reduction in the mutation rate in the 1728 subpanel compared to the 512 subpanel, possibly due to the dilution of genomic DNA from each pooled individual (64 in the 512 subpanel versus 144 in the 1728 subpanel) and a lower capacity to identify mutations, or to the lower sequencing depth in the third experiment with the 1728 subpanel (531× and 359× versus 255×). The estimated mutation rate was similar between capillary electrophoresis and amplicon-seq, both involving homozygous and heterozygous mutations. Whole-genome sequencing of randomly chosen M2 individuals revealed an average of 3.9 homozygous mutations/Mb over the whole genome or 2.64 mutations/Mb in HC gene regions. If both homozygous and heterozygous mutations are considered in the estimation, the mutation rate is three times higher than that revealed by the other two approaches. This discrepancy may be due to: (1) false positives in mutation calling by whole-genome sequencing (e.g., mutations per megabase over the whole genome versus HC genes = 3.9 versus 2.64) or (2) overstringency in mutation discrimination by capillary electrophoresis and amplicon-seq. Nevertheless, the mutation rate in the mutagenized HTX population was higher than or comparable to those in most previous studies, with the exception of the recently developed “Golden Promise” population, which was exposed to EMS mutagenesis twice before sowing of M1 seeds (Schreiber et al., 2019). In keeping with the considerable mutation rate, we observed phenotypic alterations in at least one plant of each line in 15% of all the M2 lines when grown in the field. For instance, mutants with some defect in chlorophyll content represented over 5% of the population, higher than that observed in other studies of barley (Gottwald et al., 2009; Lababidi et al., 2009; Szurman-Zubrzycka et al., 2018; Schreiber et al., 2019).

The cost and labor investment needed to identify causal mutations responsible for phenotypes of interest has fallen with the rapid development of genomics-based methodologies (Uauy et al., 2017). By employing a bulked segregant analysis, we rapidly identified an atypical EMS-induced mutation (T to A) in the HvPORB gene that explained the defect in chlorophyll content observed in one EMS mutant. The predicted HvPORB protein sequences are identical among 20 diversified barley accessions (Jayakodi et al., 2020), meaning that it is impossible to examine this gene by exploiting natural variation. By contrast, the induced mutant provided an opportunity to recognize this protein as a crucial factor in plant growth. We anticipate that our approach can be used to quickly isolate more barley genes that control growth and development, as well as end-use qualities. Some traits, such as loss of wax on the stem and leaf surface (Li et al., 2021) or spotty leaves (also referred to as hypersensitivity response mimicking; Wang et al., 2020), are rarely observed in natural germplasm. We obtained mutants with increased spike density (referred to as spikelet compactness), an agriculturally important trait that is often controlled by multiple quantitative trait loci in natural germplasm of barley (Shahinnia et al., 2012; Lu et al., 2021) and wheat (Zhai et al., 2016). However, mutant phenotypes are typically the result of a single mutation (Nie et al., 2021), thus simplifying the isolation of the causal gene.

Gene editing based on clustered regularly interspaced short palindromic repeats (CRISPR) and the CRISPR-associated nuclease Cas9 is a powerful tool for the functional validation of candidate genes and rapid genetic improvement of crop performance (Gao, 2021). Simultaneous knockouts of negative-dominant regulators in the three subgenomes of hexaploid wheat enabled the rapid creation of germplasm with desirable agronomic traits (Wang et al., 2014; Abe et al., 2019). However, genome editing tends to produce complete loss-of-function alleles at the target locus and would therefore not be feasible for housekeeping genes whose loss of function might result in plant lethality. The barley host factor eIF4E, which confers recessive resistance to virus disease (Stein et al., 2005), did not have loss-of-function variants in barley gene pools (Yang et al., 2017). In the case of eIF4E, modulation of host–virus compatibility would entail specific modifications of eIF4E residues rather than the generation of a true knockout. In this study, most EMS-induced mutations had a moderate effect on the encoded protein. Of the 510 mutations identified by amplicon-seq, only 2.75% introduced a premature stop codon or affected a splicing junction leading to a frameshift, and another 40.6% were non-synonymous substitutions that may generate an allelic series at the target gene. In addition, EMS-induced mutants may be preferable to gene-edited crops for breeding under current policies in many countries (Takagi et al., 2015).

In conclusion, the mutagenized population in the barley landrace HTX, in conjunction with an amplicon-seq-based TILLING platform and MutMap-based gene cloning, will provide a fundamental and valuable resource for barley forward and reverse genetics.

Methods

Plant materials

Four cultivated barley genotypes were included in this study: the landrace HTX and the three varieties “Yangsimai 3” (hulled grain), “Hua 30” (hulled grain), and “Beiqing 7” (hull-less grain), which are used for animal feed, in the malting industry, and as human food, respectively. M4 plants of the EMS-induced mutant M4009 were crossed with HTX wild-type plants, and the resulting F1 hybrids were selfed to develop an F2 segregating population following the speed breeding procedure (Watson et al., 2018). A mixed F2 population derived from seven independent F1 plants was obtained, consisting of 396 F2 individuals.

Generation of EMS-mutagenized population

The barley landrace HTX was used to generate a population of EMS-induced mutants. Mutagenesis induced by EMS (liquid, 99% purity; Macklin, Shanghai) treatment was performed as described previously (Gottwald et al., 2009). In brief, batches of ca. 5000 seeds were presoaked in ddH2O for 4 h at room temperature and then soaked and gently shaken for 16 h in 1 l of EMS solution (22 mM [0.28%, v/v] or 32 mM [0.4%, v/v]). The treatment was terminated by washing the seeds twice in 200 mM sodium thiosulfate solution and then five times in distilled water and drying the seed coats overnight at 4°C. M1 seeds were then sown in the field at a spacing of 10 cm between seeds. The germination and fertility rates (including semi-fertility) of M1 plants were recorded. Up to 10 M2 seeds per individual M1 plant were grown in a row under field conditions. For each row, the leaves of a single M2 seedling were sampled for genomic DNA extraction following a standard procedure (Yang et al., 2013), and the M3 seeds from the sampled M2 plant were harvested. Scoring for phenotypes was performed at the seedling and flowering stages in M2 populations and M3 families.

Library preparation, sequencing, and genome assembly

High-molecular-weight genomic DNA from HTX was subjected to genome sequencing via multiple approaches. Single-molecule real-time SMRTbell genomic libraries (insert size 40 kb) were generated and sequenced by Novogene (Beijing, China) using a PacBio Sequel II instrument. Paired-end libraries were generated for whole-genome sequencing on the Illumina NovaSeq 6000 platform. Chromosome-conformation capture sequencing (Hi-C) libraries were generated and sequenced on the Illumina NovaSeq 6000 platform at IPK Gatersleben in Germany (Padmarasu et al., 2019).

Genome assembly was performed using the software wtdbg2 (v.2.4) (Ruan and Li, 2020). FASTQ files were extracted from the raw PacBio subread BAM files with bam2fastq (v.1.1.0). The program kbm in wtdbg2 was used for self-alignment and correction of long reads (parameters: kbm2 -z 15 -c -t 60 -S 2 -k 0 -p 21 -K 1000.05). The corrected reads were aligned against all-versus-all using the program kbm2 in wtdbg2 with default parameters. After sorting of alignments by mapping quality, only the top 500 alignments for each long read as reference were retained for genome assembly using the assembly modules (program wtdbg2 -g 5g; wtpoa-cns) with the default parameters. The resulting contigs were polished with the wtpoa-cns program with one round of alignment with long reads and one round of alignment with paired-end Illumina short reads. Hi-C datasets were used to generate pseudomolecules with the TRITEX pipeline (Monat et al., 2019). Chromosomal collinearity between the “MorexV3” and HTX genome assemblies was plotted based on coding sequences (CDS) alignments. BUSCO (v.3) (Simão et al., 2015) was used to assess assembly completeness with genome assembly mode. The plant database embryophyta_odb9 was used (https://busco-archive.ezlab.org/frame_plants.html).

Gene projections of representative HC coding sequences from barley morex v3 were detected as described in Jayakodi et al. (2020). In brief, blastn and exonerate alignments of coding sequences were filtered for an initial minimal coverage of 50%. Next, a stepwise approach integrated alignment matches into the final annotation, prioritizing orthologous matches, uniqueness, match score, and completeness; further details are provided in Jayakodi et al. (2020).

TILLING by whole-genome re-sequencing

Genomic DNA of the M2 population was subjected to library construction and sequencing using the Illumina NovaSeq platform. Adapters were trimmed from raw reads; clean reads were then mapped to both the de novo assembled HTX genome and the reference “MorexV3” genome (Mascher et al., 2021) with minimap2 (v.2.17) (Li, 2018). Only reads with a mapping quality ≥20 were used for variant calling with BCFtools (v.1.9) (Li, 2011). SNP calls were subjected to the following filtration steps: (1) only biallelic variants were accepted; (2) genotype calls were considered valid if read depth was ≥20 and ≤100 for HTX (ca. 65× coverage) or ≥2 and ≤30 for other samples (ca. 9× coverage). Two SNP matrices, against either HTX or “MorexV3,” were obtained.

The rate of EMS-induced mutation was calculated based on the two SNP matrices. Only SNP sites meeting the following criteria were acceptable: (1) the genotype calls for HTX were homozygous and identical to the HTX genome reference (genotype call 0/0 in vcf format), (2) the genotype calls in the mutant samples were homozygous and corresponded to the alternate allele (1/1), and (3) the alternate allele was present in a single analyzed mutant sample, whereas the remaining mutant samples had the reference genotype or missing information. The homozygous mutation rate for each mutant was computed as the number of SNP sites with homozygous mutations divided by the cumulative size of genomic sequence with sufficient coverage (read depth ≥2 and ≤30). Heterozygous mutations were not considered. The mutation ratio was separately analyzed across the whole genome and HC gene regions.

CelI digestion and capillary electrophoresis

The classical TILLING platform that combines CelI digestion of PCR products in multiplexed DNA templates with capillary electrophoresis was used to detect mutations (Gao et al., 2020). In brief, pooling 12 M2 individual DNAs in one sample resulted in 384 DNA pools representing 4608 M2 individuals. A target fragment of the Nud gene that controls hulled/naked grain in barley (Taketa et al., 2008) was amplified by specific primers (Supplemental Table 2). The detected candidate mutations were validated by Sanger sequencing of the amplicon obtained from each of the 12 individual M2 DNAs in a multiplexed sample. The mutation rate was calculated following established methods (Gottwald et al., 2009).

TILLING by amplicon-sequencing

A three-dimensional pooling strategy was used for DNA pooling (Tsai et al., 2011; Gupta et al., 2017; Sashidhar et al., 2020). The DNA stocks (200 ng/μl) of M2 plants were pipetted into 8 × 8-well plates, with sets of eight plates comprising a subpanel of 512 individuals. Within each subpanel, the DNA samples were pooled along three dimensions: Pn, pooling of all 64 samples within plate n in equal amounts; Cn, pooling of the 64 samples from a single column n for every set of eight plates; and Rn, pooling of the 64 samples from a single row n for every set of eight plates. Alternatively, a subpanel was developed consisting of 1728 M2 individual DNA samples distributed into 12 plates with 12 × 12 wells each; the DNA pools Pn, Cn, and Rn each represented 144 individual DNAs. The pooled DNA was diluted (20 ng/μl) for PCR amplification with gene-specific primers (Supplemental Table 2). PCR products were quality controlled and quantified by electrophoresis on 1% (w/v) agarose gels, and multiple amplicons derived from the same DNA pool were mixed in equal molar amounts to generate a library. After purification and quantification using an agarose gel extraction kit and a Qubit 3.0 instrument, respectively, the libraries from pooled amplicons were subjected to amplicon sequencing following standard procedures (Sashidhar et al., 2020). The sequencing depth for each library was adjusted based on the total size of the included amplicons.

Adapters, low-quality reads, and the three border bases were trimmed from raw reads using Trimmomatic v.0.39 (Bolger et al., 2014). Trimmed reads were mapped to the amplicon sequences in the HTX reference genome by BWA-MEM v.0.7.17 (Li and Durbin, 2009). The mapped reads were sorted with SAMtools v.1.10 (Li et al., 2009), and duplicate reads were marked with the MarkDuplicates function of GATK v.4.1.9.0 (McKenna et al., 2010). Reads with ≥3 base mismatches were discarded using Perl scripts. SNP calling was performed with the HaplotypeCaller function of GATK v.4.1.9.0. Only G/C to A/T transitions were accepted, and read depths were extracted using VCFtools v.0.1.17 (Danecek et al., 2011). The depths of the mutated sites were normalized with the equation: normalized depth = Alt/(Ref + Alt) × mean depth. Alt and Ref indicate the number of reads representing the alternate allele and the reference allele, respectively. Mean depth was calculated as the number of overall reads divided by the product of the number of libraries and the mutated sites. SNPs with normalized depth ≥ 0.25 of mean depth and depth ≥5 were considered to be mutations (Supplemental Figure 2). SnpEff v.5.0e (Cingolani et al., 2012) was used to annotate and predict the effects of the identified mutation sites. The effect severity of identified mutations on the encoded proteins was predicted using the web-based tool PROVEAN (Choi and Chan, 2015).

Measurement of chlorophyll contents and photosynthetic efficiency

Germinated seeds (3 days after soaking) were incubated in the dark for 7 days at 22°C in a culture room and then exposed to light (16-h day, light intensity of 1000 μmol·m−2·s−1) (Yuan et al., 2010). Approximately 3-cm sections of the first leaf were sampled from five individuals per genotype at four time points (0 h, 24 h, 48 h, and 72 h post onset of illumination), and the content of total chlorophyll was determined for each sample using the acetone extraction method (Wei et al., 2021). Chlorophyll fluorescence and F0 and Fv values were captured with a FluorCam instrument (Photon Systems Instruments, Czech Republic). Plants at the six-leaf stage grown under normal conditions (22°C/16-h day, light intensity 600 μmol·m−2·s−1) were used to measure the net photosynthetic rate with an LI-6400XT photosynthesis system (LI-COR, USA) under a light intensity of 1000 μmol·m−2·s−1 and 500 μmol·s−1 CO2 flow conditions.

Measurement of grain storage contents

Mash from barley grains harvested at the research station in Xinxiang City, China, in 2021 was subjected to quantification of multiple constituents using the well-established near-infrared approach at the China branch (Tangshan City, China) of Cumberland Valley Analytical Services (Waynesboro, PA, USA). The standard curves were estimated based on measurements using Kjeldahl methods.

Bulked segregant analysis and genetic mapping

Two micrograms of high-quality genomic DNA from each of 60 wild-type segregants (green) and 60 chlorotic segregants in the F2 population were separately pooled to generate the wild-type and mutant bulk samples. DNA samples of the two bulk pools and the mutant parent M4009 were subjected to whole-genome sequencing and data processing following an established pipeline (Mascher et al., 2014). Genotype calls in both bulk pools were considered valid when the read depth was ≥20 and ≤100. Only the SNP sites for which the genotype call in the wild-type HTX was homozygous and identical to the reference allele were accepted. Allele frequency, which was calculated as the number of mutated reads divided by the total number of reads at a given SNP site, was calculated with a custom AWK script (Mascher et al., 2014) and visualized using standard functions of the R statistical environment (http://www.r-project.org). SNPs and indels detected in the mutant bulk for the interval 30 Mb to 323 Mb on chromosome 1H according to the reference genome “MorexV3” were annotated. KASP primers targeting the SNPs within candidate genes were generated to increase the genetic resolution by genotyping a larger F2 population.

Virus-induced gene silencing and qRT-PCR analysis

Barley stripe mosaic virus (BSMV)-induced gene silencing was carried out following an established procedure (Yuan et al., 2011). A 284-bp transcript-specific fragment was cloned into the BSMV:γb vector and transformed into Agrobacterium tumefaciens strain GV3101. The resulting Agrobacterium cultures were infiltrated into the leaves of 4-week-old Nicotiana benthamiana plants. The infiltrated leaves were collected at 7 days post inoculation to mechanically inoculate 2-week-old barley seedlings of the wild-type HTX. The second youngest leaves of HTX plants were sampled at 12 days post inoculation for total RNA extraction using the TRIzol reagent (Thermo Fisher, USA), followed by first-strand cDNA synthesis using HiScript III RT SuperMix (Vazyme, Nanjing, China). qRT-PCR was performed using the ChamQ Universal SYBR qPCR Master Mix (Vazyme, Nanjing, China) and gene-specific primers (Supplemental Table 9).

Data and mutant availability

The genome assembly of the barley landrace HTX was deposited at the European Nucleotide Archive (ENA; https://www.ebi.ac.uk/ena) under project ID PRJEB47552. The whole-genome re-sequencing data of individual mutants (ENA project PRJEB47553) or bulked segregants (ENA project PRJEB47554) have been deposited in the NCBI database. The raw data batches generated from this study have also been deposited in database hosted by the Beijing Institute of Genomics (BIG; https://ngdc.cncb.ac.cn) under project ID PRJCA008597. Queries and orders for HTX and derived mutants can be made by contacting the corresponding author P.Y.

Funding

This work was funded by grants from the (2018YFD1000702/2018YFD1000700) and the Agricultural Science and Technology Innovation Program of the (CAAS), China.

Author contributions

P.Y., M.M., Z.F., and J.Z. designed the research; C.J., M.L., G.G., L.S., Y.J., and Y.C. performed the experiments; X.Y. and C.L. contributed TILLING by capillary electrophoresis; A.H. and N.S. contributed the Hi-C data; G.H. and M.S. performed HTX genome annotation; F.D. performed the physiological characterization of the chlorotic mutant; C.J., Y.G., G.G., P.Y., and M.L. analyzed the data; S.Z., Q.H., L.L., and J.L. contributed the amplicon-seq dataset analysis; P.Y., C.J., and M.M. wrote the article.