The genotype-dependent phenotypic landscape of quinoa in salt tolerance and key growth traits.

Cultivation of quinoa (Chenopodium quinoa), an annual pseudocereal crop that originated in the Andes, is spreading globally. Because quinoa is highly nutritious and resistant to multiple abiotic stresses, it is emerging as a valuable crop to provide food and nutrition security worldwide. However, molecular analyses have been hindered by the genetic heterogeneity resulting from partial outcrossing. In this study, we generated 136 inbred quinoa lines as a basis for the molecular identification and characterization of gene functions in quinoa through genotyping and phenotyping. Following genotyping-by-sequencing analysis of the inbred lines, we selected 5,753 single-nucleotide polymorphisms (SNPs) in the quinoa genome. Based on these SNPs, we show that our quinoa inbred lines fall into three genetic sub-populations. Moreover, we measured phenotypes, such as salt tolerance and key growth traits in the inbred quinoa lines and generated a heatmap that provides a succinct overview of the genotype-phenotype relationship between inbred quinoa lines. We also demonstrate that, in contrast to northern highland lines, most lowland and southern highland lines can germinate even under high salinity conditions. These findings provide a basis for the molecular elucidation and genetic improvement of quinoa and improve our understanding of the evolutionary process underlying quinoa domestication.

1. Introduction

The decline in global crop diversity constitutes a potential threat to food and nutrition security. Recently, the increasing global vulnerability of the food supply caused by climate change and rapid population growth has made underutilized, orphan, and neglected crops potentially very important for future food and nutrition security. One such potential underutilized crop is quinoa (Chenopodium quinoa Willd.), an annual pseudocereal of the Amaranthaceae family that also includes spinach (Spinacia oleracea L.) and sugar beet (Beta vulgaris L.). Quinoa seeds and leaves have exceptional nutritional profiles, as they contain a constellation of minerals, vitamins, dietary fibre, fats, and high-quality proteins with a balanced amino acid profile. Quinoa also thrives in a wide range of marginal environments in terms of altitude, annual precipitation, temperature, and soil conditions and is thus resistant to multiple abiotic stresses such as high salinity, drought, low temperatures, and frost.

Quinoa was originally domesticated in the Andean region of South America at least 7,500 years ago and served as a staple food for indigenous people of the Andes in Pre-Columbian times. However, Spanish conquerors in the 1500s forbade the cultivation and consumption of quinoa because it was worshiped as a sacred ‘mother grain’ used in indigenous ceremonies in non-Christian religion. Therefore, in contrast to potato (Solanum tuberosum L.) and tomato (Solanum lycopersicum L.), which also originated in the Andes but expanded globally after their discovery by the Conquistadors, quinoa remained a neglected crop in South America for about 500 years until the 20th century. In 1975, the National Academy of Sciences of the United States (NAS) selected quinoa as one of the 36 underexploited plants with promising economic value. In 1993, the US National Aeronautics and Space Administration (NASA) took an interest in quinoa because it possesses desirable food qualities for astronauts on long-duration space missions. The United Nations declared 2013 as the International Year of Quinoa to raise public awareness of the nutritional benefits of this durable plant in terms of food and nutrition security worldwide. Nowadays, although quinoa is cultivated mostly in Bolivia and Peru, quinoa cultivation is being attempted in over 95 countries.

High salinity stress poses a major environmental threat to agriculture, impairing crop production in ∼20% of global irrigated land. Seed germination and subsequent seedling establishment are crucial steps in the life cycle of plants, especially under harsh environments such as salt-affected lands. Since quinoa has been recognized as a salt-tolerant crop and a potential candidate for introduction into marginal areas, many studies have focussed on salinity tolerance during germination in quinoa. However, as few studies have made use of genotyped quinoa accessions, little is known about the molecular basis of salt tolerance, a prerequisite for the improvement of crop yield on saline soils.

Quinoa is an allotetraploid (2n = 4x = 36) with an estimated genome size of 1.5 Gbp, composed of the two sub-genomes A and B. The underlying genome complexity caused by its allotetraploidy have hampered molecular analysis of quinoa genes over the years. However, the quinoa genome was recently sequenced by three independent research groups,,, opening the door to revealing the molecular mechanisms governing the development of the mysterious quinoa plant. Moreover, quinoa plants exhibit partial outcrossing due to having both hermaphrodite and female flowers on the same plant,, resulting in genetic heterogeneity that has also impeded the rigorous molecular analysis of quinoa.

To overcome this substantial limitation, we generated 136 inbred quinoa lines as a basis for molecular studies aimed at identifying genes and characterizing gene function in quinoa through large-scale genotyping and phenotyping efforts. Armed with the inbred lines, we performed a genotyping-by-sequencing (GBS) analysis and thus detected 5,753 single-nucleotide polymorphisms (SNPs) in the quinoa genome. Based on the SNPs, we show that our quinoa inbred lines are genetically divided into three sub-populations and that the highlands of South America near Bolivia and Peru are the centre of diversity of cultivated quinoa. Further, we analysed the phenotypes of the inbred quinoa lines, which we combined with our genetic data into a heatmap that allows a succinct visual overview of the genotype–phenotype relationship among the inbred quinoa lines. We also established that, unlike the northern highland lines, the seeds of the lowland and southern highland lines can germinate under high salinity conditions (600 mM NaCl). These findings will contribute to the molecular basis needed to compensate for the lost breeding opportunities and studies of quinoa over the past 500 years.

2. Materials and methods

Plant materials, generation of quinoa inbred lines, and growth conditions

Quinoa seeds (except Kd and Iw lines) were obtained from the Germplasm Resources Information Network (GRIN) of the US Department of Agriculture (USDA) (Supplementary Table S1). The sequenced and model inbred quinoa line Kd was established at Japan International Research Center for Agricultural Sciences (JIRCAS) by single-seed propagation in the absence of all other quinoa plants for over 20 years at Kyoto University, Japan. Another quinoa line used for plant virus research, Iw, was kindly provided by Dr Nobuyuki Yoshikawa, Iwate University, Japan. Inbred quinoa lines used in this study were generated by self-crossing through single-seed descent of the quinoa lines collected at JIRCAS in Tsukuba, Japan, as described previously. The seeds of inbred quinoa lines derived from the USDA were named J001–J144, excluding J062 and J070, because a USDA accession often contains seeds with multiple distinct phenotypes in terms of colour and size and none of the J062 and J070 seeds could germinate (Supplementary Table S1). For the bulking of inbred quinoa seeds, all inflorescences of each plant were covered with non-woven pollination bags (Rizo, Tsukuba, Japan) to prevent cross-pollination in phytotrons or a greenhouse.

Quinoa plants were grown on soil in pots as described previously, with minor modifications. Quinoa seeds were sown in a peat moss mixture (Jiffy Mix, Sakata Seeds, Yokohama, Japan) in a cell tray and were allowed to germinate in a growth chamber at 25°C under a 12-h-light/12-h-dark cycle (75 ± 25 µmol photons/m2/s). After 13–14 days, the seedlings were transferred to a standard potting mix (Tsuchitaro, Sumitomo Forestry, Tokyo, Japan) in 1.2-l plant pots and were grown under ambient light in temperature-controlled phytotrons at JIRCAS with a temperature of 25 ± 5°C and relative humidity of 55 ± 25%. For the bulking of inbred quinoa seeds for genotyping and phenotyping, quinoa plants were grown under ambient light in a greenhouse with a temperature of 27 ± 5°C and relative humidity of 55 ± 40%. For the phenotypic evaluation of the quinoa inbred lines at Tottori University, Japan, quinoa seeds were sown in sandy soil from the Tottori sand dunes in a cell tray and were allowed to germinate under ambient light in a greenhouse with a temperature of 18 ± 13°C. After 7 days, the seedlings were transferred to sandy soil in 1.5-l plant pots and were grown in the greenhouse.

GBS and SNP detection

We used 136 of the 142 lines for GBS analysis, since we needed to arrange the total number of samples to be a multiple of 48, including 8 control replicates and in accordance with the combination of the adapters. Total genomic DNA was extracted from fresh quinoa leaves as described previously, using the AP1 buffer from the DNAeasy Plant Mini Kit (Qiagen) instead of the lysis buffer used by Yabe et al. GBS reactions were carried out as described previously with minor modifications. EcoRI and PstI were used as restriction enzymes in this study. Barcode adaptors are listed in Supplementary Table S1. Barcode-labelled amplicons were size-selected (range 200 − 500 bp) with AMPure XP beads (Beckman Coulter) and were sequenced on an Illumina HiSeq2500 System at Macrogen Japan (Kyoto, Japan) as paired-end reads. The reads generated in this study are available from the DDBJ/EMBL/NCBI under the accession numbers listed in Supplementary Table S1. In addition to GBS data, we used published whole genome sequencing (WGS) data of quinoa, for SNP detection; the relevant accessions are listed in Supplementary Table S2. Low quality reads and adaptors were trimmed using Trimmomatic-0.36 with the options ‘SLIDINGWINDOW:4:25’ and ‘MINLEN:40’. The adaptor sequences used were ‘CACGACGCTCTTCCGATCT’ and ‘ACCGCTCTTCCGATCTGTAA’. The trimmed paired-end reads were mapped against the quinoa reference genome (Cq_PI614886_genome_V1_pseudomolecule.fa) using BWA 0.7.15 with -L 10 and -B 10 options. Single-end reads were filtered out from BAM files using samtools 1.3.1. We then carried out a local realignment procedure using RealignerTargetCreator and IndelRealigner in GATK 3.7. SNPs were detected using UnifiedGenotyper in GATK 3.7 with the option -out_mode EMIT_ALL_CONFIDENT_SITES. Heterozygous variants and nucleotide sites supported by fewer than three reads were converted to missing data using Tassel 5.2.38 and VCFtools 0.1.13, respectively. Nucleotide sites with a minimum mean depth of 10 and with proportion of max-missing data over 0.8 were filtered using VCFtools 0.1.13.

Construction of the phylogenetic tree and structure analysis

Since the depth of WGS reads was much lower than those of GBS reads obtained in this study, the WGS data were largely eliminated if the uniform filter conditions were used in all analyses. Therefore, we used different filter conditions and plant materials for each analysis (Supplementary Table S3). For STRUCTURE analysis, we used 149 quinoa lines (136 GBS and 13 WGS samples). Heterozygous variants were converted to missing data using Tassel 5.2.38. Minimum coverage depth was set to 3, and nucleotide sites with a proportion of max-missing data over 0.8 were filtered using VCFtools 0.1.13. We performed a population structure analysis using a model-based clustering method implemented in the STRUCTURE 2.3.4 software. The algorithm was run with a burn-in length of 100,000 followed by 100,000 Markov Chain Monte Carlo simulations for estimating parameters. At least 10 runs were performed to estimate the mean likelihood for a number of populations (K) from 1 to 10 using the admixture model and correlated allele frequencies. An average log-likelihood value, LnP(D), and ΔK was calculated across all runs for each K to infer the best fit to the data, as measured by the likelihood score. Lines with membership probabilities of <80% for all groups were assigned to ‘admixed’ lines. For construction of the phylogenetic tree and principal component analysis (PCA), we used GBS data from 136 quinoa lines and WGS data from 13 quinoa accessions and 2 accessions of C. hircinum (a wild relative of cultivated quinoa). The resulting dataset, consisting of 1.54 Mbp derived from 7,586 SNPs, was used for Neighbour-joining (NJ) analysis. In the NJ analysis, we used SEQBOOT with 100 replicates, followed by DNAdist (with Kimura 2-parameter method), Neighbor and Consensus programmes from the PHYLIP package 3.6 (http://www.evolution.genetics.washington.edu/phylip.htlm (last accessed 21 September, 2020)). NJ tree was rooted using midpoint rooting and visualized by FigTree 1.4.2 (http://tree.bio.ed.ac.uk/software/figtree/ (last accessed 21 September, 2020)). PCA based on covariance was performed by Tassel 5.2.38 using 7,586 SNPs.

Genetic diversity analysis and linkage disequilibrium

Nucleotide diversity (π) and Weir and Cockerham mean F-statistics (Fst) were calculated by VCFtools 0.1.13. Sliding window analysis of π, SNP frequency and number of SNPs were represented in a circular plot using the Circos 0.69. For SNP detection, we used 136 GBS samples in addition to Kd, because the sequencing depth of earlier WGS data was insufficient for these analyses (Supplementary Table S3). Pairwise linkage disequilibrium (LD) (r2) between the SNPs showing >0.05 minor allele frequency (MAF) was calculated using VCFtools 0.1.13. Average r2 values were calculated in non-overlapping 100-kbp intervals.

Germination test

Seeds were germinated in 9 cm diameter Petri dishes containing triple sheets of filter paper (Toyo Roshi Kaisha, Ltd, Tokyo, Japan), soaked with 10 ml of water or saline solution (0, 300, and 600 mM NaCl). We supplemented the 600mM NaCl saline solution with 0.2% of plant preservative mixture (Plant Cell Technology, Inc., Washington DC, USA), a broad-spectrum biocide/fungicide for plant tissue culture. Germination assays were carried out at 25 ± 1°C in the dark. We measured total hypocotyl and root length (THRL) of the germinated seedlings in 0, 300, and 600 mM saline solutions after 24 h, 4, and 14 days.

Phenotypic evaluation

Seeds from 136 lines were harvested at maturity, at which point we phenotyped six traits: 1,000-grain weight (TGW), plant height (PH), stem diameter (SD), leaf dry weight (LDW), seed yield per plant (SYP), and days to flowering (DF). TGW was calculated from the weight and the number of seeds counted by a WAVER IC-VAi seed counter (IDEX, Nagoya, Japan). SD was measured 1.0 cm above ground with a vernier caliper. LDW was measured after drying all leaves from individual plants harvested at maturity at 72°C for at least 72 h. SYP was measured on individual plants after husking seeds.

Statistical analysis of phenotypic data

A statistical analysis was performed using the R package version 3.3.3. A Shapiro−Wilk normality test assessed whether each phenotype distribution fit within the normal distribution. Statistically significant differences between phenotypic groups were determined by one-way ANOVA and Tukey−Kramer’s HSD tests (P < 0.05). We constructed the heatmap with Microsoft Excel for Mac. PCA was performed with the ‘prcomp’ function in R. The lines with missing data were eliminated prior to PCA.

Genome-wide association study

We performed genome-wide association study (GWAS) on 12 phenotypic traits (6 traits mentioned above and THRL under different conditions) using TASSEL 5.3.38. SNPs showing ≥0.05 MAF and >0.2 missing data across the 137 inbred lines were filtered from the 5,753 SNPs (see Section 2.3). Finally, we used 3,156 SNPs for general linear model (Q, GLM) and mixed linear model (Q + K, MLM) methods. Q-matrix derived from the STRUCTURE program was used for the GLM method. We performed 10,000 permutation runs to obtain the permutation-based significance in GLM analysis. The Q and kinsip (K) matrixes were used for the MLM method. K matrix was calculated based on Centered-IBS method by Tassel 5.2.38. Significant marker-trait associations (MTAs) were identified with a Benjamini−Hochberg false discovery rate (FDR) at the 5% level to correct for multiple testing.

3. Results and discussion

Generation of 136 inbred quinoa lines

Although quinoa is essentially self-pollinating, it does display genetic heterogeneity due to partial outcrossing due to the fact that individual plants carry both hermaphrodite and female flowers., Heterozygous alleles at a given genetic locus within an individual accession were seen in 32% of quinoa accessions examined. Nevertheless, to date, inbred quinoa lines have not been used in most phenotyping and genotyping studies of quinoa. Consequently, we first wished to develop inbred quinoa lines that would be a valuable resource for molecular genetics and biological studies of quinoa. We generated 136 inbred quinoa lines by self-crossing and single-seed descent of Kd, Iw, and all quinoa lines obtained from the USDA-GRIN quinoa germplasm collection and then named these inbred lines J lines (J001–J144, excluding J062 and J070) since one USDA accession often contains seeds with multiple distinct phenotypes in terms of colour and size and none of the J062 and J070 seeds could germinate (Supplementary Table S1). We used pollination bags to prevent cross-pollination when plants set seeds. The ratio of heterozygous sites to all SNPs sites ranged from 0.0024 to 0.019% across all lines, with an average value of 0.0075% (Supplementary Table S1, also see below), indicating that our inbred lines were purified enough for genotyping and phenotyping analyses.

SNP discovery by GBS

To evaluate genetic diversity in quinoa, we used a GBS approach to identify SNPs in all inbred quinoa lines. We sequenced DNA libraries derived from EcoRI/PstI double digests of genomic DNA from our 136 quinoa inbred lines on the HiSeq2500 platform. After primary quality filtering steps, we generated a total of 560 million clean reads (covering 56.6 Gbp), with an average of 8,723,466 reads (corresponding to 880.3 Mbp) per line (Supplementary Table S1). On average, 88.0% of trimmed paired-end reads correctly mapped to the reference genome. The average depth for all nucleotide sites in the GBS samples was 148.1. After filtering the nucleotide sites detected by GBS for all 136 samples and by WGS for Kd, we retained 1,709,689 nucleotide sites, which correspond to 0.14% of the reference genome (Supplementary Table S4). From these 1,709,689 nucleotide sites, we detected 5,753 SNP sites, ranging from 88 (Chr13) to 719 (Chr07) per chromosome. The number of SNPs per chromosome showed a high correlation with the length of each assembled chromosome (R2 =0.93), as expected. These results suggest that the GBS reads are properly mapped to the reference sequence, since the low mis-mapping rate of GBS reads is presumably due to the high nucleotide divergence between the A and B genomes, (pair-wise nucleotide divergence at synonymous site = 0.08−0.1). SNPs were distributed over the entire genome (Fig. 1). The average SNP density across all chromosomes was 0.0034 (i.e., one SNP per 297 bp) and ranged from 0.0025 (one SNP per 407 bp, Chr17) to 0.0041 (one SNP per 243 bp, Chr04) depending on the chromosome. These data indicate that the 5,753 SNPs we obtained by GBS can be used as genome-wide markers for population and phylogenetic analyses, as described below.

Figure 1

Distribution of SNPs and nucleotide diversity in the quinoa genome. (A) Distribution of SNP frequency, summarized in non-overlapping 1-Mbp windows along the quinoa genome. Lines along the y-axis are shown in increments of 0.005. (B) Distribution of SNP number, summarized in non-overlapping 1-Mbp windows. Lines along the y-axis are shown in increments of 5. (C) Comparison of nucleotide diversity (π) between highland and lowland lines. Yellow, red, and blue indicate all, lowland, and highland lines, respectively. (D) Comparison of nucleotide diversity (π) between the three populations. Red, green, and blue indicate lowland, northern highland, and southern highland lines, respectively. Mean π was calculated for partially overlapping windows of 2.5 Mbp with a 1-Mbp step. Lines along the y-axis are shown in increments of 0.0005.

Population structure, genetic diversity, and evolutionary history of domestication in quinoa

To investigate the population structure within C. quinoa, we conducted multiple analyses (STRUCTURE, NJ phylogenetic, and PCA) using 5,529 SNPs obtained from both GBS data for 136 lines and WGS data for 13 accessions. Assignment analysis using the program STRUCTURE indicated that the highest ΔK value was observed for K = 2 (15,466) followed by K = 3 (1,490), whereas all larger Ks (K ≥ 4) had substantially lower ΔK values (Figs. 2A and B). With K = 2 (Fig. 2C), one population was mainly composed of lines from Bolivia and Peru, whereas the other population largely contained lines from Chile and a set of US accessions that lack reliable passport data (Supplementary Table S5). These two populations corresponded to those named highland and lowland (Coastal), respectively, in previous reports,,,, (Fig. 2D). Many of the admixed lines (lines with membership probabilities of <80%) were lowland lines (Fig. 2C). This might indicate low and/or recent gene flow had occurred from highland to lowland population. With K = 3 (Fig. 2C), the highland population was further divided into two sub-populations, corresponding to southern and northern highland sub-populations.,, The southern highland sub-population was mainly composed of lines from Bolivia, whereas the northern highland sub-population contained lines from Peru. Notably, although we identified a number of admixtures consisting of lowland and southern highland genetic background, we failed to detect any cases of admixtures between lowland and northern highland genetic material (Fig. 2C), suggesting a closer relationship between southern highland and lowland populations than between northern highland and lowland populations, as might be expected due to physical isolation (Fig. 2D). Considering that these data are consistent with previous reports,,, it would be more plausible that recent gene flow mainly occurred from the Bolivian population to the Chilean population (Fig. 2D).

Figure 2

Population structure of 149 quinoa lines. (A) Plots of mean likelihood at each K value. (B) Plots of ΔK at each K value. (C) Model-based ancestry for each quinoa line. Each individual is shown as a vertical bar. The models with K = 2 and K = 3 show the high ΔK value. At K = 2, red and blue indicate lowland and highland lines, respectively. At K = 3, red, blue, and green indicate lowland, southern highland, and northern highland lines, respectively. (D) Distribution of quinoa in South America. Estimated areas of northern/southern highland and lowland populations are indicated by green, blue, and red, based on the collection sites (circles) of the accessions with the passport data associated with the inbred lines from this study. Red arrows indicate the estimated evolutionary process of domestication in quinoa (see Section 3 for details).

We also conducted NJ analysis and PCA to examine the phylogenetic relationships among C. quinoa lines, using two C. hirucinum accessions (a wild relative of domesticated quinoa) as outliers. These two different approaches strengthen our view that quinoa lines are divided into three clusters: Lowland, northern highland, and southern highland (Figs 2D and 3; Supplementary Fig. S1). Four inbred lines (J010, J027, J032, and J040) assigned to southern highland were included in the lowland cluster of the phylogenetic NJ tree (Fig. 3), supporting the notion that the lowland cluster is more closely related to southern highland lines than it is to the northern highland cluster. One quinoa line (O654) was positioned in an intermediate position between the C. hircinum accessions BYU566 and BYU1101 in the phylogenetic tree (Fig. 3). In addition, line O654 had much higher heterozygosity than the other lines (0.6496 vs 0.0075% for the across-line average; Supplementary Table S2), suggesting that O654, reported in the previous study, may be a hybrid or a hybrid-descendant from a cross with the other species. It is noteworthy that the phylogenetic analysis placed both C. hircinum accessions BYU566 and BYU1101 at the base of all quinoa lines analysed in this study (Fig. 3). Neither C. hircinum accessions were included in the three sub-populations revealed by the PCA (Supplementary Fig. S1). These findings suggest that all quinoa lines may have originated from a single domestication event (Fig. 2D).

Figure 3

NJ tree of 149 quinoa lines and two C. hircinum accessions based on 1.54 Mbp of sequence data corresponding to 7,586 global SNP positions. Red, blue, and green indicate lowland, southern highland, and northern highland lines, respectively. Numbers above branches show bootstrap values based on 100 replicates (not shown if <80%). Red asterisks indicate bootstrap values of 100%. The scale bar corresponds to 2×10−4 substitutions per nucleotide site.

Next, we evaluated population differentiation using Weir and Cockerham’s Fst value. Fst between highland and lowland populations was 0.209, compared with 0.111 between northern and southern highland populations (Supplementary Table S6). These numbers indicate that highland and lowland lines are clearly differentiated, corresponding to the highest ΔK (ΔK = 15,466) at K = 2 obtained by STRUCTURE. For K = 3, Fst between southern highland and lowland populations (0.241) was lower than that between the northern highland and lowland populations (0.291), supporting the view that the lowland population is more closely related to the southern highland population than to the lowland population. We also estimated the genome-wide nucleotide diversity within each population: the nucleotide diversity (π) within the highland population (4.73 × 10−4 for northern highland; 4.86 × 10−4 for southern highland) was higher than across the lowland population (3.90 × 10−4) (Supplementary Table S7). While previous studies had estimated higher genetic diversity in the highland quinoa population based on limited sequence comparison and isozyme banding patterns, our genome-wide SNPs data clearly reinforce the notion that South American highland, including Lake Titicaca in the High Andes on the border between Bolivia and Peru, is the centre of diversity from which cultivated quinoa originated.

Together, the results derived from our genetic analyses based on GBS and WGS data support the single domestication hypothesis:,, quinoa was most likely first domesticated near Lake Titicaca, which straddles the border between the northern and southern highland, and then spread to the northern and southern highlands, finally reaching the Chilean lowland region via the southern highland (Fig. 2D). However, a new question arises from our observations: how did the lowland population undergo high genetic differentiation relative to the southern highland population? Our findings may support the possibility that lowland quinoa populations might be derived from crosses between quinoa and wild species, such as C. hircinum, in the Chilean lowland region. In rice (Oryza sativa L.), two diversified subspecies of cultivated rice (japonica and indica) arose from a single domestication event (the japonica subspecies), and indica was later derived from crosses between japonica and wild rice species. However, earlier phylogenetic data obtained by WGS indicated that BYU566 clustered with the lowland population, suggesting that quinoa may have independently been domesticated in the highland and lowland. There are differences in both the number of quinoa lines (149 in this study vs 16 in the previous report) and in the methods (GBS + WGS vs WGS) between this study and the previous report, but it remains unclear whether modern-day quinoa originated from a single or two independent domestication events. Therefore, further identification of domestication-related genes such as shattering genes, and/or a genome-wide population analysis using unbiased large sample of quinoa and C. hircinum will be required to definitively determine the evolutionary process of domestication in quinoa.

Phenotypic analyses

3.4.1. Salt tolerance during the germination stage

During the plant life cycle, germination is an early and critical growth stage that determines whether plants will survive and thrive at that location. Here, we therefore focussed on salt tolerance during the germination stage as an initial evaluation of salt tolerance in quinoa lines. In total, we assessed 137 inbred quinoa lines for tolerance to high salt concentrations during germination (Supplementary Tables S8 and S9). We observed that southern highland lines had a much longer mean THRL when compared with lowland and northern highland lines after 24 h under control conditions without added NaCl (Supplementary Fig. S3), indicating that southern highland lines grew more rapidly than all other lines. This may suggest that rapid germination ability is important determinant whether plants can grow in arid regions such as Bolivian Southern Altiplano around Uyuni Salt Flats, where the frequency of rainfall occurrence is very low. However, we failed to detect significant differences in mean THRL between the three population groups after 4 days of growth in the presence of 300 mM NaCl (Supplementary Fig. S3). In contrast, we noticed that most northern highland lines did not germinate after 14 days of treatment with 600 mM NaCl, whereas lowland and southern highland lines did, suggesting that these lines are tolerant to germination under high salinity conditions (Supplementary Fig. S3). Interestingly, we saw no correlation between germination potential at 300 and 600 mM NaCl treatments (R2=0.01) (Supplementary Fig. S4). These results strongly suggest that germination tests to determine salt tolerance require high salinity conditions (600 mM NaCl) in order to differentiate quinoa lines on the basis of their salt tolerance.

3.4.2. Growth traits

We also assessed 137 inbred quinoa lines for important growth traits such as TGW, PH, SD, LDW, SYP, and DF (Supplementary Tables S8 and S9). The average TGW from the southern highland and lowland lines was higher than for the northern highland lines regardless of the location and season, suggesting that northern highland lines produced smaller seeds and that TGW is a comparatively conserved trait in quinoa (Supplementary Fig. S3). For the other traits examined in this study, mean values for PH, SD, LDW, SYP, and DF were higher in the lowland lines than in the other lines (Supplementary Fig. S3), consistent with the fact that lowland accessions originated in the Chilean coastal region, which is roughly located on the same latitude as Japan conducted this study.

Phenotypic diversity

To examine the relationship between population structure and phenotypic diversity, we performed a PCA using data for the phenotypic traits. The individual plotting factors for the first two PCs are listed in Supplementary Table S10. PC1 and PC2 accounted for 34.0 and 28.8% of the total phenotypic variance, respectively, and the following variables showed substantial loadings (loading coefficients ≥0.3): THRL (600 mM NaCl at 14 days) and TGW for PC1; four variables of plant size, namely PH, SD, LDW, and DF, for PC2. The PCA plot exhibited partially overlapping distributions for the three sub-groups (Fig. 4): PC1 separated the northern highland lines from the others, whereas PC2 separated the lines into highland and lowland lines, therefore following altitude as a discriminating factor.

Figure 4

PCA of quinoa lines based on data for the phenotypic traits. Scatterplot of the first two components (PC1 and PC2) from the PCA using 123 quinoa lines with complete phenotypic data. The proportion of variance explained by each component is given in parentheses along each axis. TGW, thousand-grain weight.

The PCA plot and associated loading coefficients (Supplementary Table S10) clearly suggest that northern highland lines are characterized by lower salt tolerance and smaller grain size, whereas most southern highland and lowland lines are highly salt-tolerant during germination and produce larger grains (Fig. 4). In addition, most lowland lines displayed larger plant size and late flowering phenotypes, whereas southern highland lines had smaller plant size and early flowering phenotypes under our experimental conditions in Japan (Fig. 4). These data therefore demonstrate that quinoa lines should be classified into three, rather than two populations in terms of both genotype and phenotype. Notably, combining the NJ tree and the heatmap illustrating the phenotypes examined here offered a visual confirmation of the separation of quinoa lines and associated phenotypes into the three phylogenetic sub-groups indicated by the PCA (Fig. 5). Furthermore, this figure provides quick access to the genotype–phenotype state of each quinoa line (Fig. 5).

Figure 5

Heatmap of the phenotypic traits in quinoa lines. The columns of the heatmap are ordered based on the NJ tree, which is shown on the left with its branch lengths transformed to be proportional. Missing values in the heatmap are shown in white. Red, green, and blue dotted lines represent lowland, northern highland, and southern highland lines, respectively. THRL, the average total hypocotyl and root length (mm) of seeds treated with 0, 300, or 600 mM of NaCl for each indicated time; TGW, thousand grain weight (g); PH, plant height (cm); SD, stem diameter (mm); LDW, leaf dry weight (g); SYP, seed yield per plant (g); DF, days to flowering (see Supplementary Tables S8 and S9 for details).

Based on our genotype data, and assuming a single domestication event, (Fig. 2D), the phenotypic landscape of quinoa lines in terms of THRL in the presence of 600 mM NaCl indicates that cultivated quinoa originated near Lake Titicaca, and then gained high salt tolerance during germination in the saline soils around the Uyuni Salt Flats in Bolivia, and then spread to lowland regions in Chile, where it retained its high salt tolerance during germination (Fig. 2D). Therefore, the combination of genotype and phenotype data allows us to easily imagine the evolutionary history of quinoa domestication. However, since the sample size for experiments performed under strict and uniform conditions was small due to limitations of our facilities, future analyses with larger sample sizes and fewer lines are needed to confirm the results obtained here and to elucidate the mechanisms underlying salt tolerance and other key traits in quinoa. Overall, the visual data of the heatmap, combined with the phylogenetic tree, aids our understanding of the phenotypic landscape in quinoa for important growth traits (Fig. 5).

Association mapping of the phenotypic traits

Although association mapping, which is based on the concept of LD, has been used for detecting quantitative trait loci (QTL) in a wide range of organisms, the rate of LD decay is unknown for quinoa. Therefore, before investigating the association between phenotype and genotype by GWAS, we first evaluated LD decay/extension. We calculated the pairwise r2 correlation value between SNPs with an MAF of at least 0.05. Lowland lines showed a much greater LD decay rate than highland lines (Supplementary Fig. S5). Indeed, the distance between two SNPs around which LD dropped below 0.1 was 1 Mbp for lowland lines and 5.5–8 Mbp in highland populations. Such LD decay distance was greater than in other self-pollinating crops such as rice (indica, ca. 65–75kb; japonica, ca. 150–200 kb), and sorghum (Sorghum bicolor, ca. 150 kb), indicating lower resolution in association mapping for identifying QTLs of quinoa.

Despite this slow LD decay rate, we performed an association mapping analysis with 3,156 SNPs and 12 phenotypic traits in our quinoa population. We identified 3,091 and 4 significant MTAs (P < 0.01 for GLM, FDR <0.05 for MLM), respectively, during association mapping using the GLM and MLM methods (Supplementary Table S11 and S12). The MLM method detected far fewer MTAs than the GLM method. This discrepancy suggests that the MLM may yield fewer false-positive results than the GLM, because the MLM used the kinship as well as the population structure. The strong population structure and well-differentiated phenotypes among populations might underlie the much smaller number of MTAs identified by MLM. Taking into account the low LD decay rate as well as the strong population structure, these observations suggest that GWAS may not be a suitable technique for identifying QTLs in quinoa.

4. Future perspectives

Quinoa has been recognized as a valued crop with the potential to contribute to global food security due to its exceptional nutritional profile and great adaptability to harsh environments. There is thus an increasing effort to expand quinoa production into non-native areas worldwide.,, To date, however, low heat tolerance, stem lodging, negative effects of long photoperiods, and high susceptibility to biotic stresses such as diseases and pests have been major limiting factors in recent quinoa expansion to areas where the pseudocereal is not adapted and for effective large-scale cultivation., Conversely, the High Andes is a major production area that provides the world with high-quality quinoa, but they are also considered to be marginal arable land (with yearly rainfall below 200−500 mm; over 3,700 m above sea level) and are in the front line of climate change. Indeed, the High Andes faces severe problems in quinoa production, mainly caused by extreme weather events induced by climate change, such as drought and frost. Thus, quinoa breeding is now necessary to yield maximization for potential production and yield optimization for production stability inside and outside its countries of origin. In this light, our current findings provide a useful basis for the molecular elucidation and genetic improvement of quinoa. One hundred and thirty-six newly developed inbred quinoa lines, which consist of three clusters grouping lines into northern highland, southern highland, and lowland sub-populations, will be powerful resources not only for quinoa breeding but also for the elucidation of the molecular mechanisms underlying its great tolerance to abiotic stresses and high nutritional value. Our genotype-informed phenotypic analysis of quinoa under high salt conditions and other important growth traits will greatly contribute to the molecular dissection of traits and breeding in quinoa. We hope that our work will contribute to fill in the gap of the past 500 years in quinoa breeding and studies spanning the 16th to the 20th century, in sharp contrast to other crops that originated from the Andes, including potato, tomato, chili pepper (Capsicum annuum L.), cotton (Gossipium hirsutum L.), and tobacco (Nicotiana tabacum L.), which have all expanded worldwide via the Spaniards since the 16th Century.

Recently, quinoa has been studied as a model plant to understand the molecular mechanisms of abiotic stress tolerance in plants. In the case of salt tolerance, the roles and functions of inorganic ions for osmotic adjustment, tonoplast channels, epidermal bladder cells,, and xylem ion loading have been reported. Recently, the role of peroxisome proliferation in heat and drought stress and heat stress responses in quinoa shoots have been also demonstrated. Apart from these, other molecular analyses concerning evolution of quinoa genome, the decision to flower, inflorescence patterning, and abscisic acid (ABA)-mediated seed dormancy and germination have also been reported.

To date, though, the quinoa research community lacks the resources of gain- and loss-of-function technologies to study gene function in situ, which constitutes a significant limiting factor for molecular analysis in quinoa. There is therefore an urgent need to develop transformation and genome editing techniques in quinoa to allow the dissection of gene function. In contrast to rice and maize (Zea mays L.), which require high fluence rates to keep plants healthy, quinoa plants grow very well in growth chambers normally used for the model plant Arabidopsis (Arabidopsis thaliana L.), thereby enabling researchers to shift their focus from the weedy model plant Arabidopsis to the new model plant quinoa with great practical value and appropriate inbred lines. Further technical development for the study and breeding of quinoa should facilitate breeding of other staple crops based on the mechanisms of tolerance to abiotic stresses described in quinoa.

Data availability

The raw Illumina reads used in this study are available from DDBJ/EMBL/NCBI under the accession numbers listed in Supplementary Table S1.

Supplementary data

Supplementary data are available at DNARES online.

Click here for additional data file.