Genomic insight into the origin, domestication, dispersal, diversification and human selection of Tartary buckwheat

Abstract

Background: Tartary buckwheat, Fagopyrum tataricum, is a pseudocereal crop with worldwide distribution and high nutritional value. However, the origin and domestication history of this crop remain to be elucidated.

Results: Here, by analyzing the population genomics of 567 accessions collected worldwide and reviewing historical documents, we find that Tartary buckwheat originated in the Himalayan region and then spread southwest possibly along with the migration of the Yi people, a minority in Southwestern China that has a long history of planting Tartary buckwheat. Along with the expansion of the Mongol Empire, Tartary buckwheat dispersed to Europe and ultimately to the rest of the world. The different natural growth environments resulted in adaptation, especially significant differences in salt tolerance between northern and southern Chinese Tartary buckwheat populations. By scanning for selective sweeps and using a genome-wide association study, we identify genes responsible for Tartary buckwheat domestication and differentiation, which we then experimentally validate. Comparative genomics and QTL analysis further shed light on the genetic foundation of the easily dehulled trait in a particular variety that was artificially selected by the Wa people, a minority group in Southwestern China known for cultivating Tartary buckwheat specifically for steaming as a staple food to prevent lysine deficiency.

Conclusions: This study provides both comprehensive insights into the origin and domestication of, and a foundation for molecular breeding for, Tartary buckwheat.

Keywords: Artificial selection; Buckwheat; Domestication; Genomics; Migration.

Fig. 1

Geographic distribution, population structure and genomic diversity of Tartary buckwheat accessions. Geographic distributions of 567 Tartary buckwheat accessions. The radius of each pie represents the sample size in each region and the colors indicate the proportions of HW (Himalayan wild accession), SL1 (Southwest landrace 1), SL2 (Southwest landrace 2), NLI (Northern Landrace-Within China), NLO (Northern Landrace-Outside China). XZ, Xizang Province; SC, Sichuan Province; YN, Yunnan Province; GZ, Guizhou Province; HuB/HN/JX, Hubei/Hunan/Jiangxi Province; HB/NM/LN, Hebei/Inner Mongolia/Liaoning Province; SNX/SX, Shannxi/Shanxi Province; QH/GS/NX, Qianghai/Gansu/Ningxia Province. B The maximum-likelihood phylogeny of 567 Tartary buckwheat accessions and model-based clustering analysis with different numbers of ancestry kinship (K= 2-6). Different colors indicate different groups based on the population structure. C PCA plots of 567 Tartary buckwheat accessions and outgroup. Colors represent the membership at K = 6 (Fig. 1b). D Nucleotide diversity (π; within circles) and population divergence (F_ST; between circles) for the five groups (the outgroup population was not included). E Group-specific LD decay plots

Fig. 2

Demographic history and dispersal of Tartary buckwheat. A Heatmap showing the similarity of five population through outgroup f₃ matrix. B Divergence times of the five populations. The range of predicted divergence time was shown. C Outgroup f₃ statistics biplot measuring genetic similarity. Diagonal line marks the f₃ statistics for G2/G5. Different groups representing accessions collected from different areas. G1, Himalayan region; G2, Sichuan; G3, Yunnan; G4, Guizhou; G5, Qinghai-Gansu; G6, Inner Mongolia-Hebei; G7, Hunan-Hubei-Jiangxi; G8, Poland; G9, Slovenia and G10, France. D Phylogenic tree of outgroup, Himalayan located group (G1), northern China located group (G5-G7) and outside China located group (G8-G10). E Pairwise fixation index (F_ST) of the mini-groups of Tartary buckwheat. F Gene flow between populations estimated using Treemix. Yellow and orange lines between populations indicate gene flow. G The possible spread of Tartary buckwheat from its origins in the Himalayas. Ten groups representing the population along the route are indicated. The average predicted divergent times are shown

Fig. 3

Variation of FtGULO controls disease resistance during Tartary buckwheat domestication. A-B Selective sweeps identified through comparisons between HW and SL (A) and HW and NL (B) using XP-CLR (cross-population composite likelihood-ratio test). The dashed line represents the top 5% of values therefore scores in these regions were regarded as selective sweeps. C Local Manhattan plot of GWAS signals on Chr 8 for resistance to R. solani AG4-HGI 3. The dashed line represents the threshold (-log₁₀P = 5). D Schematic diagram of FtGULO gene structure. Two SNPs in the promoter of FtGULO were marked as red letters and result in haplotypes (Hap) A and T. E Box plots show disease index in plants carrying the two haplotypes (Hap). n_Hap-A = 8, n_Hap-T = 234. P values were calculated using a two-tailed t-tests. F Expression of FtGULO in accessions harboring the two haplotypes. Error bars indicate the ± s.d., n = 6. Significance was tested using one-way ANOVA. G Transcription activity of FtGULO promoters with two haplotypes. H Disease index of accessions among HW, NL and SL groups. n_HW = 10, n_NL = 96, n_SL = 140. Significant was tested using two-tailed t-tests. *, P < 0.05. I Frequencies of the two haplotypes in the HW, NL and SL groups. J Subcellular localization of FtGULO-GFP fusion protein transient expression in N. benthamiana leave cells. Scale bars, 10 µm. (K-L) Relative expression levels of FtGULO during R. solani infection (K) and MeJA treatment (L). Histone H3 was used as the internal reference. M Disease index of Arabidopsis lines heterologously expressing FtGULO. Significant differences were identified using one-way ANOVA. n = 6. N Phenotypes of Arabidopsis WT lines and lines heterologously expressing FtGULO with and without infection with R. solani AG4-HGI 3. Scale bars, 1 cm

Fig. 4

Variation of FtPK controls salt resistance differences between north and south populations of Tartary buckwheat. A Selective sweeps identified through comparisons between SL and NL using XP-CLR (cross-population composite likelihood-ratio test). The dashed line represents the top 5% of values therefore scores in these regions were regarded as selective sweeps. B Manhattan plot of GWAS signals for salt resistance in Tartary buckwheat accessions. The dashed line represents the threshold (-log₁₀P=5). C Schematic diagram of FtPK gene structure. Two SNPs in the promoter of FtPK are marked with red letters and result in haplotypes (Hap) 1 and 2. D Box plots show salt resistance in two haplotypes (Hap). n_Hap-1 = 13, n_Hap-2 = 120. P value was calculated using two-tailed t-tests. E The expression level of FtPK in accessions with the two haplotypes. The error bars indicate the ± s. d, n = 6. The P value was calculated using one-way ANOVA. (F) Transcription activity of FtPK promoters with two haplotypes. G Differentiation salt resistance of accessions among HW, NL and SL groups. n_HW = 7, n_NL = 93, n_SL = 51. Significant differences were tested using two-tailed t-tests. *, P < 0.05. H Frequencies of the two haplotypes in the HW, NL and SL groups. I Confocal microscope image showing nuclear localization of FtPK-GFP fusion protein upon transient expression in N. benthamiana leaf cells. Scale bars, 10 µm. J Phenotypes of Arabidopsis lines heterologously expressing FtPK and subjected to salt stress. K Root length of Arabidopsis lines heterologously expressing FtPK and subjected to salt stress. Significant differences were tested using two-way ANOVA with Tukey HSD test. There was an effect of treatment (F = 11.044, df = 1, P = 0.004) and an effect of genotype (F = 4.478, df = 3, P = 0.018)

Fig. 5

Structural variation of FtXIP controls the domestication of easily-dehulled type Tartary buckwheat. A Genome features of EDT. The outermost circle represents each chromosome of the genome. The second to fifth circles indicate gene density, SNPs density, deletion density, and insertion density, respectively, using a window size of 500-kb. B Gene dot map between easily-dehulled type buckwheat (EDT) and difficult-dehulled type (DDT) Tartary buckwheat. C Diagram representing the generation of the EDT x DDT recombinant inbred lines (RILs). D Genome wide Δ(SNP index) plot of the population derived from a cross between EDT and EDT. The black lines indicates tricube-smoothed Δ(SNP index), and the gray lines indicate corresponding two-sided 99% confidence intervals. E Insertions and deletions larger than 50 bp and within 5 kb of genes in the chr 2 QTL intervals. F Expression of genes with insertions and deletions in the QTL intervals in the seed coats of EDT and DDT at the 20-day after pollination (DAP) stage. Each small square represents the differentially expression level of a gene between EDT and DDT. Square with gene ID exhibited the differentially expressed genes. The red gene ID represents FtXIP. G Schematic diagram showing the deletion of 1,140 bp in the promoter region of MqXIP gene. H Transient expression assay was conducted to compare the transcription activity of MqXIP and an empty vector. I The expression level of XIP in DDT and EDT Tartary buckwheat. The error bars indicate the ± s. d, n = 6. The P value was calculated using two-tailed t-tests. P < 0.05