Background

Rice production in Vietnam is of great value for export and providing daily food for more than 96 million people. However, agricultural production, especially rice cultivation, is inherently vulnerable to climate variability across all regions in Vietnam. Based on the records of monthly precipitation and temperature from 1975 to 2014 (Nguyen et al. 2019), the areas of highest crop production in the low lying Mekong and Red River Deltas are particularly vulnerable to the increasing threat from climate change. In 2017, the total planted area of rice in Vietnam was 7.7 million hectares. This includes 4.2 million hectares in the Mekong River Delta and 1.1 million hectares in the Red River Delta (GSO-Database 2017). These are also the areas where most of the population of the county is concentrated. In the Mekong River Delta, the damaging effects of salinisation and drought to rice production have increasingly manifested themselves in recent years (Parker et al. 2019; Son et al.  2018; Tran et al. 2019; Yen et al. 2019).

Vietnam possesses a vast diversity of native and traditional rice varieties due to its geographical situation, latitudinal range and diversity of ecosystems (Fukuoka et al. 2003). This diversity constitutes a largely untapped and highly valuable genetic resource for local and international breeding programs. Vietnamese landraces are disappearing as farmers switch to modern elite varieties. To limit this erosion of genetic resources, several rounds of collection of landraces, particularly from the northern upland areas, have been undertaken since 1987. Thousands of rice accessions have been deposited in the Vietnamese National Genebank at the Plant Resources Center (PRC, Hanoi, Vietnam), together with passport information detailing their traditional name and province of origin. One hundred and eighty-two traditional Vietnamese accessions were selected for a genotype by sequencing (GBS) study in 2014 (Phung et al. 2014). This study yielded 25,971 single nucleotide polymorphisms (SNPs) that were used to describe four Japonica and six Indica subpopulations. These subpopulations were classified by region, ecosystem and grain-type using passport information (province and ecosystem) and phenotyping. This dataset had subsequently been used for genome-wide phenotype-genotype association studies (GWAS) relating to root development (Phung et al. 2016), panicle architecture (Ta et al. 2018), drought tolerance (Hoang et al. 2019b), leaf development (Hoang et al. 2019a) and jasmonate regulation (To et al. 2019) and phosphate efficiency (Mai et al. 2020; To et al. 2020).

An international effort to re-sequence Asian rice accessions known as the “3000 Rice Genomes Project” (3K RGP) has provided the rice community with a better understanding of Asian rice diversity and evolutionary history, as well as providing valuable knowledge to enable more efficient use of these accessions for rice improvement (Wang et al. 2018; Wing et al. 2018). However, only 56 of these accessions originated from Vietnam, suggesting that the rice diversity within this country may not be fully captured within the 3K RGP. While the original 3K RGP analysis described nine subpopulations (Wang et al. 2018), subsequent reanalysis had shown that the 3K RGP could be further subdivided into fifteen subpopulations (Zhou et al. 2020).

The primary objective of this study was to gain a clear understanding of the rice population and genome structure in Vietnam, leveraging the vast diversity of rice varieties due to the reasons previously highlighted, and then to contextualize this analysis in the 3K RGP. In addition, we aim to assess the phenotypic variability between subpopulations and expand on previously published QTLs using a larger dataset. For this, we newly sequenced 616 Vietnamese rice accessions using whole-genome sequencing (WGS), most of them being native landraces. One hundred sixty-four of these rice accessions were in common with a previous study [8] based on a genotyping-by-sequencing (GBS) approach. We supplemented this dataset with all 56 Vietnamese genotypes from the 3K RGP to form a native diversity panel with 672 accessions. We analysed this diversity panel of 672 accessions to explore how breeding and environmental pressures have shaped the rice genome in Vietnamese accessions. We also carried out a comprehensive analysis of the population structure of the 3635 rice genomes obtained from joining our diversity panel and the complete 3K RGP datasets. We completed a GWAS on the diversity panel with 672 accessions (and separately for the Japonica and Indica subtypes within it) on thirteen phenotypes, which are available for around two-thirds of the samples. Our results highlight genomic differences and trait associations in traditional Vietnamese landraces, which are likely the product of adaption to multiple environmental conditions and regional culinary preferences in a very diverse country.

Results

Sequencing Rice Diversity from Vietnam

Whole-genome sequencing was carried out on 616 rice accessions. Five hundred eleven of the accessions were obtained from the PRC (Plant Resource Centre, Hanoi, Vietnam, http://csdl.prc.org.vn), together with their passport data, which shows that they were collected from all eight administrative regions of Vietnam (Table S1). The remaining samples were obtained from AGI’s collection (Agricultural Genomics Institute, Hanoi, Vietnam). Three reference accessions (Nipponbare, a temperate Japonica; Azucena, a tropical Japonica; and two accessions of IR64, an Indica) obtained from the PRC, were included in the dataset. A total of 1174 Giga base-pairs (Gbps) of data was generated for the 616 samples representing an average sequencing depth of 30x for 36 “high coverage” samples and 3x for 580 “low coverage” samples (Table S1). These 616 newly-sequenced accessions were classified into 379 Indica and 202 Japonica subtypes, with the remaining 35 (including the Aus and Basmati varieties) being classified as admixed, based on the STRUCTURE (Pritchard et al. 2000) output for K = 2 using a subset of 163,393 SNPs.

Population Structure of Rice within Vietnam

The population structure of rice within Vietnam was analysed using the diversity panel of 672 samples, comprising 616 newly sequenced accessions and 56 Vietnamese genotypes from the 3K RGP. We assigned the 672 samples to four Japonica subpopulations and five Indica subpopulations (Table S1) using (i) the population structure information obtained from the STRUCTURE analysis (Fig. 1), (ii) the previous characterisation of a panel of Vietnamese native rice varieties using GBS (Phung et al. 2014), and (iii) the assessment of the optimal number of subpopulations (Fig. S1) using the method described in Evanno et al. (2005). Subpopulations were named as in Phung et al. (2014), except that we considered the I6 subpopulation to be part of the I3 subpopulation. Although the previous study used a limited number of GBS markers, 129 of the 164 common samples were assigned to the same subpopulations in both studies. Most differences were due to samples being classified as admixed in either one of the studies. We classified 48 (11%) of the Indica (Im), and eight (4%) of the Japonica samples (Jm) as admixed. The reference varieties Nipponbare (Temperate Japonica), Azucena (Tropical Japonica), and IR64 (Indica) were classified as J4, J1 and I1, respectively.

Fig. 1
figure 1

Population structure and location of the Indica and Japonica subpopulations within Vietnam. a STRUCTURE results (mean of 10 replicates) at K = 4 for 211 Japonica subtypes. The cut off for inclusion in each subpopulation is 0.6. The number of samples in each subpopulation is shown above, a further 8 samples were classified as admixed. b STRUCTURE results (mean of 10 replicates) at K = 5 for 426 Indica subtypes. Each colour represents one subpopulation. Each accession is represented by a vertical bar and the length of each coloured segment in each bar represents the proportion contributed by each subpopulation. The cut off for inclusion in each subpopulation is 0.6. The number of samples in each subpopulation is shown above, a further 48 samples were classified as admixed. c STRUCTURE results for the I5 subpopulation expanded to show individual samples. d The proportion of each population originating from each of the 8 regions in Vietnam (based on a subset of 377 samples, 54% of Indica samples and 85% of Japonica samples)

Full size image

Each Indica subpopulation contained shared ancestry (admixed components) with other Indica subpopulations (Fig. 1a). The admixed components are shown in detail for the 43 samples in the I5 subpopulation (Fig. 1c) namely 38 samples from our dataset and the following five samples from the 3 K RGP; IRIS 313–11384 (IRGC 127275), B184 (IRGC 135862), IRIS 313–11383 (IRGC 127274), IRIS 313–10751 (IRGC 127577) and IRIS 313–11893 (IRGC 127519). The Japonica subtropical J1 subpopulation shared ancestry (between 0 and 25% of the genome) with the Japonica tropical J3 subpopulation, whereas the two temperate subpopulations, J2 and J4 shared ancestry dominantly with each other. The tropical J3 subpopulation contained four samples with around 20% of the haplotypes in common with the temperate J4 subpopulation. Using the passport information available from the PRC, the proportion of each subpopulation originating from each of the “administrative regions” of Vietnam is shown in Fig. 1d. Only the I1 and I2 Indica subpopulations were collected from the Mekong River Delta regions, I2 being almost exclusively grown there whereas I1 was more widespread than I2. The I4 and J4 subpopulations were mainly collected from the Red River Delta areas. The J1 and J3 subpopulations were closely related; the J1 subpopulation was predominantly from the North of Vietnam whereas the J3 subpopulation was concentrated around the South-Central Coast region. Small variations in the percentage of reads mapping were observed for each of the subpopulations (Fig. S2).

A Principal Component Analysis (Fig. 2a and b) showed the relationship between these nine Vietnamese subpopulations. Concerning the Vietnamese genotypes from the 3K RGP dataset included in the diversity panel, the Indica I1 subpopulation included two XI-1B modern varieties and eight admixed (XI-adm) accessions. I2 included fourteen XI-3B1 genotypes, which comprises Southeast Asian accessions, and similarly, I3 and I4 included one and ten XI-3B2 genotypes, respectively. Finally, I5 included five XI-adm accessions and clustered distinctly away from all the other subpopulations (Fig. 2a). On the other hand, J1 included the two subtropical (GJ-sbtrp) accessions from the Vietnamese 3K RGP genotypes, and J3 included one tropical (GJ-trp1) accession from the Vietnamese 3K RGP genotypes (Fig. 2b). These results correlate well with the latitudinal distinction between these subpopulations. J2 and J4 included two and one temperate (GJ-tmp) accessions, respectively; and split into two clear subpopulations in Vietnam compared with the East Asian temperate subpopulation described by the 3K RGP.

Fig. 2
figure 2

PCA analysis of Indica and Japonica Vietnamese subpopulations. a PCA analysis of 426 accessions from Vietnam using the top two components to separate the five Indica subpopulations. The ellipses show the 95% confidence interval. b PCA analysis of 211 accessions from Vietnam using the top two components to separate the four Japonica subpopulations. The ellipses show the 95% confidence interval

Full size image

Population Structure of the Combined 3635 Asian Cultivated Rice Genomes

Six hundred twelve of the 616 newly sequenced accessions from this study and the 3023 accessions from the 3K RGP were combined and classified into 9 and 15 subpopulations (Table S2), and compared with the subpopulations from the 3 K RGP analysis (Wang et al. 2018; Zhou et al. 2020). For clarity, we used the prefix Jap- and Ind- to label these subpopulations from our analysis.

When the combined dataset of 3635 samples was classified into nine subpopulations (Fig. S3a), we found that 95% of the 3 K RGP accessions (2882 out of 3023) were assigned into the same subpopulations. The remaining 5% lines were either (i) previously classified as admixture and our analysis placed into a subpopulation, or (ii) were previously classified in a subpopulation and were now classified as admixture. The 612 newly sequenced Vietnamese accessions were placed in three Indica clusters (187 accessions), three Japonica clusters (176 accessions), the Basmati and Sadri aromatic cB group (11 accessions), or the Aus cA subpopulation (one accession). In more detail, the three Indica clusters included three Im accessions in the East Asian cluster (Ind-1A), seventy-six I1 accessions in the cluster of modern varieties of diverse origins (Ind-1B), and 108 accessions (I2, I3 and Im) in the Southeast Asian cluster (Ind-3). Whereas, the three Japonica clusters included 54 accessions (J2, J4 and Jm) in the primarily East Asian temperate cluster (Jap-tmp), 119 accessions (J1, J3 and Jm) in the Southeast Asian subtropical cluster subpopulation (Jap-sbtrp) and three J3 accessions in the Southeast Asian Tropical subpopulation (Jap-trp). Any remaining accession with admixture components over 65% either Indica or Japonica were classified as Ind-adm (191 accessions) or Jap-adm (27 accessions), respectively. Finally, the remaining accessions were considered as Admix (19 accessions). Notably, all thirty-seven I5 accessions were placed in Ind-adm, and ten of the sixteen J3 accessions were placed in Jap-adm.

When the combined dataset of 3635 samples was reclassified into 15 subpopulations (K15_new, Fig. S3b), we noticed the following differences in the distribution of subpopulation compared to the 3K RGP analysis for the same number of 15 subpopulations (K15_3KRGP); we did not observe the division of the Aus samples into cA-1 and cA-2, and we subdivided the Indica subtypes and Japonica subtypes into eight and five subpopulations, respectively. A Principle Coordinate (PCO) analysis of the Indica and Japonica subpopulations is shown in Fig. 3, highlighting our new eight Indica and five Japonica subpopulations (In addition the Vietnamese and 3K RGP subpopulations are shown in Figs. S5 and S6).

Fig. 3
figure 3

PCO analysis of Indica and Japonica Vietnamese subpopulations. a PCO analysis of 1605 Indica samples (omitting the samples classified as XI-adm and Ind-adm outside Vietnam for clarity). The ellipses show the 95% confidence interval for the K15_new subpopulations (the K15_3KRGP and five Vietnamese Indica subpopulations are shown in Fig. S5). X = PC1, Y=PC4, Z = PC5. b PCO analysis of 982 Japonica samples (omitting the samples classified as GJ-adm and Jap-adm outside Vietnam for clarity) showing the K15_new subpopulations (the K15_3KRGP and four Vietnamese Japonica subpopulations are shown in Fig. S6) X = PC3, Y=PC4, Z = PC5

Full size image

The relation between the subpopulations in our comprehensive analysis (3635 accessions) and the 3 K RGP (3023 accessions) was as follows: (i) The Ind-1A, Ind-1B.1 and Ind-1B.2 were equivalent to XI-1A, XI-1B1 and XI-1B2, respectively. Forty-three of the Vietnamese I1 accessions were in the Ind-1B.1 subpopulation, and the remaining 102 I1 accessions were classified as admixed. (ii) The Ind-2 was equivalent to XI-2A and XI-2B, and as expected, this geographically distant South Asian subpopulation was not present in Vietnam. (iii) The previously observed split of the Indica-3 subpopulation into 3A and 3B was also observed in our analysis, where Ind-3.1 was equivalent to XI-3A and did not contain any Vietnamese accessions. (iv) The remaining Ind-3.2, Ind-3.3 and Ind-3.4 were a rearrangement of the XI-3B1 and XI-3B2 subpopulations. (v) The 89 Vietnamese I2 accessions belonged to Ind-3.2, which was a subset of XI-3B1. (vi) Ind-3.3 contained 16 of the 37 Vietnamese I3 accessions. (vii) 72% of the accessions in Ind-3.4 were from Vietnam, which contained 13 of the 37 I3 accessions, 61 of the 62 I4 accessions, and all I5 accessions. Within Ind-3.4, the admixture components of I3, I4 and I5 subpopulations (Fig. S7) showed that I3 accessions were highly admixed, some I4 and I5 accessions were completely within Ind-3.4, while other I4 and I5 accessions showed admixture with Ind-3.3 (I5) or Ind.2, Ind-3.2, and Ind-3.3 (I4). To clarify these relations, a principle component analysis (PCA) with a reduced number of accessions was carried out using the 723 sample dataset (672 Vietnamese accessions and 51 genotypes from neighbouring Southeast Asian Countries; Fig. S8), this supported the close relationships of I2 with XI-3B1, I4 with XI-3B2, I5 with XI-adm, J1 with GJ-sbtrp, and that both J2 and J4 were within GJ-tmp.

Phenotypic and Genetic Diversity Analysis of the Vietnamese Indica and Japonica Subpopulations

Phenotypic measurements for 19 traits were scored in field conditions in the Hanoi area by breeders from the Agricultural Genomics Centre (AGI) for approximately two-thirds of the samples in our study. For five of these traits, additional scores were also included from trials by the Vietnamese Plant Resource Centre. In addition, phenotypic data were available for eleven of the traits in 38 of the 56 samples sourced from the 3K-RGP dataset (Table S3, S4). Finally, the grain length to grain width ratio (GL/GW) was calculated to give a total of 20 traits (Table S5). Scores were available for between 328 and 503 of the 672 samples (Indica subpanel, 170–297 samples and Japonica subpanel, 134–178 samples).

Using a t-test, significant differences in measurements were found between the Indica and Japonica subtypes for ten of the traits; these are detailed in Table S5 and histograms are shown in Fig. 4 for selected phenotypes. The Indica subtypes had significantly (p-value < 0.0001) higher values for grain length to width ratio, leaf pubescence, culm number, culm length, and floret pubescence. In contrast, the Japonica subtypes had significantly higher values for grain width, leaf width, flag leaf angle, panicle length, and floret colour. The Indica I1 subpopulation (mostly elite varieties) was the most phenotypically distinct when compared to the rest of the Indica samples (mostly native landraces). I1 samples had longer grains (p-value = 2.2e-16), earlier heading date (p-value = 9.9e-12), higher culm strength (p-value = 2.2e-16), shorter leaf length (p-value = 2.7e-14) and shorter culm length (p-value < 2.2e-16). Similar values were obtained when comparing I1 to just the I5 subpopulation (Fig. 4). The I5 subpopulation was not phenotypically distinct (p-value < 0.001) from the other landrace subpopulations I2, I3 and I4, except for a significantly lower measurement of leaf pubescence (p-value = 0.0007). The Japonica J2 subpopulation had a significantly lower grain length to width ratio than J1 (p-value = 1.8e-13) and J3 (p-value = 5.7e-07). A correlation analysis carried out between the 20 phenotypes (Fig. S9) showed that the highest correlation (r = 0.6) was between leaf length and culm length (excluding the correlation between grain length to width ratio and grain length and grain width). Histogram and correlation plots are available for the 13 traits used for the GWAS analysis in Fig. S10 comparing the Indica and Japonica subtypes and in Fig. S11 comparing subpopulations I1 and I5. Further boxplots showing the phenotypic distribution according to subpopulation for culm length, grain length, grain width and heading date are available in Fig. S12.

Fig. 4
figure 4

Histograms comparing the Indica and Japonica subtypes and the I1 and I5 subpopulations. Histogram are shown for 8 of the 13 traits used in the GWAS analysis. The Japonica and Indica subtypes are shown in green and purple respectively and underneath a histogram is shown for a subset of the Indica values comparing subpopulations I1 and I5. The mean value is shown by a dotted line and the p value (T-test) is shown at the top of each plot. A ggpairs histogram and correlation plot is available for all 13 traits in Figs. S7 and S8

Full size image

The Japonica subtypes had a lower nucleotide diversity (π = 0.000912) than the Indica subtypes (π = 0.00167). Looking at the individual subpopulations (Table S6), the elite I1 subpopulation is the most diverse (π = 0.00144), and the I5 subpopulation is the least diverse (π = 0.00103). Regions of the genome with low diversity in all Indica subpopulations, and regions with low diversity in specific subpopulations, were observed when plotting diversity along each chromosome (Fig. S13). The J3 subpopulation is the most diverse of the four Japonica subpopulations. (π = 0.000697). Large genomic regions with very low diversity were observed in chromosomes 2, 3, 4 and 5 in all Japonica subpopulations (Fig. S14).

Genome-Wide Genotype-Phenotype Association Analysis

Three independent GWAS were conducted using the full panel (672 samples, 361,191 SNPs), the Indica subpanel (426 samples, 334,935 SNPs) and the Japonica subpanel (211 samples, 122,881 SNPs). Thirteen (13) of the 20 traits were suitable for GWAS based on the variance (Coefficient of Variation < 56% for the full panel). The full list of phenotypic measurements is available in Table S3. We found 643 significant phenotype-genotype associations. These associations were organised into 21 QTLs (Table 1, Table S7). The GWAS Manhattan and Quantile-Quantile plots are available in Fig. S17 and Fig. S18.

Table 1 21 QTLs identified for plant description traits in the full panel, and Indica and Japonica subpanels. Detailing for the QTL analysis; significance threshold -log10 (p value) ≥ 8.0; panel in which significant associations were detected, highest level of significance for all panels and overlaps with publish QTLs for Vietnamese rice populations or for the 3 K RGP
Full size table

We used both the STRUCTURE and kinship matrix to control for population structure, however the QQ-plots suggest that p-values are still inflated, especially for the full panel, regarding several traits, such as GL GW ratio, heading date, diameter internode, and floret pubescence. We used a strict cut off of -log10(p) > 8.0 and required at least two significant SNPs per QTL, as detailed in Table 1.

The 21 QTLs contained 1730 genes and covered a total of 11 Mbp over ten chromosomes, and contained 453 SNPs with a significant association to a trait in at least one diversity panel (Fig. 5). The list of genes within each QTL is available in Table S8.

Fig. 5
figure 5

The distribution of 21 QTL. 21 significant associations for 8 of the 13 traits (−log10 (p value) ≥ 8.0). The 33 individual associations for the full panel and the Japonica and Indica subpanels were merged to form the 21 final QTLs. The QTLs for grain length, grain width and grain length/width ratio were merged into QTLs for grain size, these are labelled in brown. The remaining QTLs are labelled in black; Leaf width (LW), Panicle Length (PL), Heading Date (HD), Floret Pubescence (FP), Diameter Internode (DI). Regions smaller than 100 kb are extended to 50 kb either side of SNP with maximum p value. Centromeric regions are shown as 100 kb regions in dark grey

Full size image

Seventeen QTLs were identified in the full diversity panel significantly associated with eight traits: grain length, grain width, grain length-to-width ratio, leaf width, panicle length, floret pubescence, heading date and internode diameter. A further 4 QTLs associated with grain length and grain width were observed only in the Japonica subpanel. Three of the QTLs, which were found in the full panel, were also observed in the Indica subpanel.

The set of 3.8 M SNPs (see methods), representing one SNP every 99 bases, was annotated based on the potential effect of each SNP on protein function using SnpEff (Table S11). 526,138 (4.79%) of the SNPs were in genes. There were 21,639 (0.197%) SNPs in 11,125 genes classified as having a putative “High impact” effect (E.g. Exon changes, frameshifts, gene fusions or rearrangements, protein structural changes, etc.). Following additional minimal allele frequency (MAF) filtering, in the Indica dataset (MAF 5%, 2,027,294 SNPs), there were 11,906 “High impact” SNPs in 7396 genes and in the Japonica dataset (MAF 5%, 1,125,716 SNPs), there were 6240 “High impact” SNPs in 4439 genes of which 2818 were present in both Indica and Japonica.

None of the 453 SNPs with a significant association was annotated as resulting in protein changes (“High impact” SNPs). However, “High impact” effects were identified in other SNPs within the QTL. Among the total 1730 genes in the 21 QTLs, we annotated 309 genes with “High impact” SNPs in the Indica subpanel, 248 genes with “High impact” SNPs in the Japonica subpanel, including 137 “High impact” SNPs common between the two sets. In addition, we looked for overlaps with the QTL in five published Vietnamese studies (Hoang et al. 2019a; Hoang et al. 2019b; Phung et al. 2016; Ta et al. 2018; To et al. 2019), which used 25,971 SNPs in 182 samples (164 in common). We found that 2_GL and 6_GS overlapped with QTL for panicle morphological traits (Ta et al. 2018); 2_GL overlapped with QTL9 for secondary branch number, and spikelet number (SBN and SpN), and 2_GS overlapped with QTL12 for secondary branch average length (SBL). 4_GW_jap overlapped with “q1” for longest leaf length (LLGHT) (Phung et al. 2016).

Discussion

Indica and Japonica Rice Subpopulations within Vietnam

Whole-genome sequencing of 616 Vietnamese rice accessions, predominantly landraces, plus 56 Vietnamese genotypes previously sequenced by the 3K RGP, provides us with a diversity panel to clarify the structure of rice subpopulations in Vietnam. Here, we describe five Indica subpopulations and four Japonica subpopulations using phenotypic measurements from this study, passport information available from the Vietnamese National Genebank (PRC), and the agronomic and geographical annotations from Phung et al. (Phung et al. 2014). In general terms, our population structure within Vietnam agreed with the previous study, which used a smaller number of markers and 182 samples and is approximately a third of our diversity panel (Phung et al. 2014). Subpopulation I1 is the most phenotypically distinct of the Indica subpopulations and shows typical phenotypes of ‘elite’ varieties, such as short height, strong culm strength, long slender grains and a short growth-duration (less than 120 days from sowing to harvest). I1 accessions are grown throughout Vietnam in irrigated ecosystems but predominantly in the Mekong River Delta in the south of the country. Subpopulation I2 is mainly composed of long growth-duration (over 140 days), tall varieties grown in the rainfed lowland and irrigated ecosystems of the Mekong River Delta with a broad diversity of grain shapes. The remaining three Indica subpopulations are intermediate between I1 and I2 for growth-duration, height and culm strength, have a broad diversity of grain shapes, and are not grown in the Mekong River Delta. Subpopulation I3 has the highest proportion of upland varieties but also includes some lowland varieties from the “South Central Coast” region many of which were classified as an independent subpopulation (I6) by Phung et al. (Phung et al. 2014). Subpopulation I4 is mainly grown in the rainfed lowland and irrigated ecosystems of the Red River Delta. Subpopulation I5 is grown in a range of ecosystems but concentrated around the North Central Coast and Red River Delta regions, but excluding the Northwest region suggesting that it is the main lowland subpopulation. The J1 and J3 subpopulations are closely related upland varieties and the J2 and J4 subpopulations are closely related lowland varieties. Subpopulation J1 is mostly composed of medium growth-duration upland varieties from the mountainous regions in the North of Vietnam, with long large grains typical of upland varieties. Subpopulation J2 is grown throughout Vietnam in a range of ecosystems but has consistently short grains. Subpopulation J3 is mainly grown in the “South Central Coast” region and has long large grains. Subpopulation J4 is primarily grown in the Red River Delta region in lowland and mangrove ecosystems and has short grains.

Root traits measured by Phung et al. (2016) can be used to gain some information about the drought tolerance of the subpopulations. The J1 and J3 upland subpopulations have deeper and thicker roots than the thinner shallower roots in the J2 and J4 subpopulations, which are grown in irrigated and mangrove ecosystems (Phung et al. 2016). This suggests that the J1 and J3 subpopulations, which are grown mainly in rainfed upland regions, would be more drought tolerant than the others. Similarly, the I3 subpopulation has the deepest and thickest roots. It would, therefore, be more drought tolerant than the I1 and to a lesser extent the I5 subpopulation, which has the thinnest, shallowest root systems.

A Comprehensive Analysis of the Available 3635 Asian Cultivated Rice Genomes

The comprehensive analysis of the combined 3635 Asian cultivated rice genomes obtained by joining our diversity panel with the full 3K RGP dataset resulted in a similar assignation to the previous 3K RGP analysis in 84% of the cases. The largest differences were that the 3K RGP split the cA and XI-2 subpopulations, while our analysis split the GJ-tmp and rearranged the two XI-3B subpopulations into Ind-3.2, Ind-3.3 and Ind-3.4. The single temperate subpopulation (GJ-tmp) from the 3K RGP is further split in our study between the Jap-tmp.1 and Jap-tmp.2 subpopulations, with 88% of the samples in Jap-tmp.2 coming from Vietnam and forming the J2 subpopulation. These differences are likely due to changes in the distribution of genetic variants in subpopulations expanded within Vietnam.

Vietnamese Rice Subpopulations in the Context of the 3K RGP Asian Cultivated Rice Subpopulations

The Indica I1 subpopulation, which contains a high proportion of elite varieties, clustered with the X1-1B1 subpopulation of modern varieties. The Southeast Asian native subpopulations (XI-3B1 and XI-3B2) clustered with the I2 and I4 subpopulations, respectively. I3 appeared to include both XI-3B1 and XI-3B2 accessions. The subpopulations from East and South Asia (XI-1A, XI-2A, XI-2B, XI-3A) had no representatives from Vietnam and fell outside of the Vietnamese subpopulation clusters, as expected. Our four Vietnamese Japonica subpopulations relate to the tropical (J1), subtropical (J3) and temperate (J2 and J4) Japonica subpopulations from the 3K RGP according to their latitudinal origin from South to North Vietnam, respectively.

The most exciting subpopulation is I5. When all 3635 samples were considered, the subpopulation XI-3.4 included half of the I3, all but one of I4 and all I5 Vietnamese accessions, as well as half of the Southeast Asian native XI-3B2 genotypes from the 3K RGP. The remaining XI-3B2 were classified as Indica admix (Ind-adm). However, when only the Vietnamese samples were considered in the analysis, I5 clustered distinctly away from I3 and I4 subpopulations (Fig. 2a) and included five accessions from the 3K RGP, which had very low shared ancestry (admixture components) with other 3K RGP samples. Notably, Vietnamese landrace IRIS 313–11384 (IRGC 127275) had little shared ancestry with any other Vietnamese 3K RGP genotypes. Remarkably, a recent study on genomic signals of admixture and alien introgression in a core collection of 948 accessions representative of the earlier Asian Rice Landraces (Santos et al. 2019) included IRIS 313–10751 (IRGC 127577) and IRIS_313–11383 (IRGC 127274) from the I5 subpopulation.

Genome-Wide Association Analysis in Vietnamese Rice Landraces Highlighted 21 QTL

We have also extended upon seven published GWAS (Hoang et al. 2019a; Hoang et al. 2019b; Mai et al. 2020; Phung et al. 2016; Ta et al. 2018; To et al. 2019; To et al. 2020), which focussed on specific traits but used a smaller number of markers and a third of the samples from the Vietnamese dataset. We took a similar approach of carrying out the analysis on both the full panel and the Indica and Japonica subpanels. There are some interesting overlaps between the QTLs  from the various studies, notably, the overlap of QTL for panicle morphology with our QTL for grain size (2_GL and 6_GS). These previous studies found QTL in the full panel and in the Indica subpanel, but not in the Japonica subpanel. However, we found QTL for grain size that were only present in the Japonica subpanel, and all the QTL found in the Indica subpanel were also found in the full panel. These differences probably reflect our larger dataset. Comparing our results with the GWAS results from the 3K RGP (https://snp-seek.irri.org/) (Mansueto et al. 2017; Mansueto et al. 2016), the QTL 5_GS on chromosome 3 is in the same region as a marker associated with grain length, and the QTL 10_GS on chromosome 5 is in the same region as a marker associated with both grain width and grain length. Underlying these two QTL, there are genes that have a putative role in the control of grain size in rice (Li et al. 2018), namely GS3 (Os03g0407400) in 5_GS and GSE5 (LOC_Os05g09520, Os05g0187500) in 10_GS. Functional nucleotide polymorphism (FNP) can be caused by either a SNP or Indel, as has been shown for a number of yield related traits such as grain size (Kim et al. 2016). However, we were only able to detect SNPs in our dataset. We also looked for genes with “High impact” SNPs in QTL, relevant candidates include bip130 (Zhou et al. 2019) (LOC_Os05g02260, Os05g0113500) with a stop gain mutation underlying the QTL 9_PL for panicle length and OsSPX-MFS3 (LOC_Os06g03860, Os06g0129400) (Wang et al. 2015) with a splice acceptor variant at the end of an intron underlying the QTL 11_GL for grain length.

Subpopulation I5 Constitutes an Untapped Resource of Cultivated Rice Diversity

The analysis restricted to Vietnamese accessions allowed us to observe differences among the accessions within the country. Although 38 accessions (including two genotypes from the same accession in our study) are deposited in the PRC in Hanoi, and the remaining five accessions are available from the 3 K RGP, there is limited information from the passport and phenotypic data to be able to understand the distinctiveness of this subpopulation fully. Further analysis of this subpopulation should encompass ‘Indica specific genes’ which may have been overlooked in our study as we used a Japonica reference.

Phung et al. (Phung et al. 2014) described subpopulation I5 as “medium growth-duration accessions from various ecosystems of the North and South Central Coast regions, with rather small and non-glutinous grains”. Our I5 accessions are predominantly from the Red River Delta and contiguous coastal departments, the “North Central Coast” and “Northwest” administrative regions, but remarkably excluding the higher altitude Northwest region in the North, the more upper “Central Highlands”, as well as the whole Mekong River Delta in the south. This suggests that I5 accessions are common traditional low yielding lowland varieties with specific environmental or culinary values.

Comparing the Vietnamese subpopulations to the fifteen Asian rice subpopulations identified from the 3K RGP highlighted the I5 subpopulation as a potential source of novel variation as it forms a well-separated cluster. Subpopulation I5 originates from lowland areas such as the Red River Delta and adjacent regions. For the range of phenotypes measured in this study, the I5 subpopulation did not differ phenotypically from the other landraces, which have undergone breeding selection within Vietnam. However, compared to the ‘elite’ I1 subpopulation, I5 accessions have shorter grains, take longer to flower, having lower culm strength, longer culms and leaves.

Conclusions

In this study, we generated a large genome-variation dataset for rice by sequencing 616 accessions from Vietnam and supplemented these with data obtained from the 3K RGP. Using this resource, we incorporated the Vietnamese rice diversity within the population structure of the Asian cultivated rice. Firstly, we incorporated ~ 50 representative samples into our dataset to define the five Indica and four Japonica subpopulations found in Vietnam using STRUCTURE population analysis. We then added a further ~ 50 samples from outside Vietnam and carried out a PCA analysis. Subsequently, we merged both our Vietnamese and the 3K RGP datasets and used admixture for the global population analysis. These approaches showed comparable population structure.

A GWAS analysis yielded associations for grain characteristics, panicle length, heading date and leaf width. These traits are likely under selection and associated with yield. Together with previously published QTLs using a subset of our Vietnamese rice dataset (Hoang et al. 2019a; Hoang et al. 2019b; Phung et al. 2016; Ta et al. 2018; To et al. 2020; To et al. 2019), we obtained a comprehensive set of QTLs that can be used to further understand the breeding potential of varieties in Vietnam. The locally adapted varieties and subpopulations provide a source of novel alleles that can be exploited in rice breeding to develop a new generation of sustainable ‘Green Super Rice’ with lower input needs, enhanced nutritional content and suitability for growing on marginal lands (Wing et al. 2018). Also to develop climate-smart rice varieties, which are of particular relevance to Vietnam’s high production areas in the low lying deltas, which are currently severely impacted by climate change.

Materials and Methods

Sequencing of 616 Accessions from Vietnam

We sequenced a total of 616 rice accessions, 612 accessions from Vietnam and three reference accessions, Nipponbare, a temperate Japonica; Azucena, a tropical Japonica; and IR64, an Indica (2 samples). Five hundred eleven accessions are available from the Vietnamese National Genebank (PRC) at http://csdl.prc.org.vn (Table S1). All Vietnamese native rice landraces were grown at Dai Dong Experimental Farm (Dai Dong commune, Thach That district, Hanoi, Vietnam) in 2015. The healthy seeds generated from one mature spikelet of the individual plant in each landrace were harvested and dried separately. After that, the selected seeds (35–40 seeds/landrace) were incubated and sown for 2 weeks to collect leaf samples (30 g/sample) for genomic DNA extraction. Total genomic DNA extraction of each rice landrace was made from young leaf tissue using the Qiagen DNeasy kit (Qiagen, Germany). DNA concentration and purity of the samples were measured by the UV-VIS NanoDrop ND-2000 spectrophotometer (Thermo Fisher Scientific) at OD 260/280 nm and OD 260/230 nm wavelengths.

Sequencing was performed by Genomic Services at the Earlham Institute (Norwich, UK). Around 1 μg of genomic DNA from each sample was used to construct a sequencing library. For the 36 high coverage samples (prefix: SAM) the Illumina TruSeq DNA protocol was followed, and the samples were sequenced on the HiSeq 2000 for 100 cycles. For the low coverage samples (prefix: LIB), genomic DNA was sheared to 500 bp using the Covaris S2 Sonicator (Covaris and Life technologies), and samples were processed using the KAPA high throughout Library Prep Kit (Kapa Biosystems, MA, USA). The ends of the DNA were repaired for the ligation of barcoded adapters. The resulting libraries were quality checked, pooled, and quantified by qPCR. The libraries were sequenced on a HiSeq 2500 instrument following the manufacturer’s instructions.

Phenotyping

Phenotyping experiments were conducted at the Thach That Experimental Farm of AGI in 2014 and 2015 (Dai Dong commune, Thach That district, Hanoi, Vietnam). The seeds of each rice landrace were incubated in an oven at 45 °C for 5 days to break the seed dormancy. All rice seeds were soaked in tap water for 2 days and incubated at 35-40 °C for 4 days for germinating. The fully germinated seeds of each rice landrace were directly sown in the paddy field plot (1.5m2 in the area). After 15 days of sowing, 24 seedlings of each landrace were carefully transplanted by hand in field plots (2x4m2). The fertiliser and pesticide applications were performed following the conventional methods of rice cultivation in Vietnam. The phenotypic and agronomic characteristics were carried out following IRRI’s standard evaluation system (IRRI 2002).

In addition, phenotypic data were available for eleven of the traits in 38 of the 56 genotypes sourced from the 3K RGP dataset. These eleven traits were included in our analysis because we did not observe a significant difference (p-value > 0.07) between our dataset and the 3K RGP dataset for the I2 subpopulation (Table S5).

Merging the SNP Called in the Sequenced Materials and the Complete 3 K RGP Dataset

Raw sequencing reads were mapped to the Nipponbare reference genome Os-Nipponbare-Reference-IRGSP-1.0 (IRGSP-1.0), using BWA-MEM with default parameters except for “-M -t 8”. Alignments were compressed, sorted and merged using samtools. Picard tools were then used to mark optical and PCR duplicates and add read group information. We used freebayes v1.1.0 for variant calling using default parameters. A total of 21.2 M variants were identified of which 16.4 M were SNPs, and 4.8 M were indels. The resulting VCF file was then filtered for biallelic SNPs with a minimum SNP quality of 30, resulting in 16.0 M variants. PLINK v1.9 was used to convert the VCF into a PLINK BED format. These variants were then combined with the 3K RGP 29 M biallelic SNPs dataset v1.0 by downloading the PLINK BED files from the “SNP-seek” database (https://snp-seek.irri.org) excluding variants on scaffolds and 26,553 SNPs that were flagged as triallelic upon merging, resulting in 36.9 M SNPs. The SNPs present in both datasets were then extracted and filtered using an identical approach to Wang et al. (Wang et al. 2018), resulting in 5.9 M SNPs. For that, PLINK v1.9 “--hardy” (Purcell et al. 2007) was used to obtain observed and expected heterozygosity for 100,000 SNPs. We removed SNPs in which heterozygosity exceeds Hardy–Weinberg expectation for a partially inbred species, with inbreeding coefficient (F) estimated as the median value of “1 − Hobs/Hexp”, in which Hobs and Hexp are the observed and expected heterozygosity for SNPs where “Hobs/Hexp < 1” and the minor allele frequency is > 5% and using the cut-off value of 0.479508 for the entire 3622 samples dataset. A further filtered set of 3.4 M SNPs was obtained by removing SNPs with > 20% missing calls and MAF < 1%. Finally, a core set of 361,279 SNPs was obtained with PLINK by LD pruning SNPs with a window size of 10 SNPs, window step of one SNP and r2 threshold of 0.8, followed by another round of LD pruning with a window size of 50 SNPs, window step of one SNP and r2 threshold of 0.8. Samples with more than 50% missing data in this core set were then removed, resulting in dropping seven newly sequenced samples and one genotype from the 3K RGP dataset.

Population Structure of the Combined 3635 Samples

The population structure was analysed using the ADMIXTURE software (Alexander and Lange 2011) on the SNP set obtained in the previous section. First, ADMIXTURE was run from K = 5 to K = 15 in order to compare it with the analysis from IRRI (Wang et al. 2018; Zhou et al. 2020). For each K, ADMIXTURE was then run 50 times with varying random seeds. Each matrix was then annotated using the subpopulation assignment from the 3K RGP nine subpopulations. Then, up to 10 Q-matrices belonging to the largest cluster were aligned using CLUMPP software (Jakobsson and Rosenberg 2007), these were averaged to produce the final matrix of admixture proportions. Finally, the group membership for each sample was defined by applying a threshold of ≥0.65 to this matrix. Samples with admixture components < 0.65 were classified as follows. If the sum of components for subpopulations within the major groups (Ind and Jap) was ≥0.65, the samples were classified as Ind-adm or Jap-adm, respectively, and the remaining samples were deemed admixed (admix).

Multi-dimensional scaling analysis was performed using the ‘cmdscale’ function in R (R Core Team, 2020), using a distance matrix obtained in R using the ‘Dist’ function from the amap package (Lucas 2018). The resulting file was then passed to Curlywhirly (https://ics.hutton.ac.uk/curlywhirly/) and rgl v0.100.19 (https://r-forge.r-project.org/projects/rgl/) for visualisation.

Recalling the Diversity Panel with 723 Samples

The 616 rice samples were mapped to the Japonica Nipponbare (IRGSP-1.0) reference with BWA-MEM using default parameters, duplicate reads were removed with Picard tools (v1.128) and the bam files were merged using SAMtools v1.5 (Li et al. 2009). Variant calling was completed again on the merged bam file with FreeBayes v1.0.2 (Garrison and Marth, 2012) separately for each of the 12 chromosomes, but using the option “--min-coverage 10”. Over 6.3 M bi-allelic SNPs with a minimum allele count of ≥3 and quality value above 30 and missing in < 50% of samples were obtained with VCFtools v0.1.13 (Danecek et al. 2011). BAM alignment files to the Nipponbare IRGSP 1.0 reference genome were downloaded from http://snp-seek.irri.org/ (Mansueto et al. 2017; Mansueto et al. 2016) for 107 selected samples. Alignment statistics are included in Table S11. These BAM files were merged and variant calling was similarly completed using FreeBayes v1.0.2 (Garrison and Marth, 2012) separately for each of the 12 chromosomes using the option --min-coverage 10, and filtered with VCFtools v0.1.13 as before to obtain 6.8 M bi-allelic SNPs with a minimum allele count of ≥3 and quality value above 30 and missing in < 50% of samples. The two sets of 6.3 M and 6.8 M SNPs were merged using BCFtools v1.3.1 isec to obtain 4.4 M SNPs which were present in both sets and in at least 70% of samples. These 4.4 M SNPs were then filtered to remove positions which fell outside the expected level of heterozygosity for this dataset, as previously indicated. The resulting estimate of F for the 723 samples was 0.882, so a SNP whose heterozygosity is >5x higher than the most likely value for a given frequency and the dataset’s inbreeding rate will be deemed as having an excessive number of heterozygotes. The cut-off value was 0.591, which resulted in 3.8 M SNPs passing this filter, a scatter plot indicating the SNPs which were kept and removed is shown in Fig. S15. Missing data was imputed in this latest dataset using Beagle v4.1 with default parameters (Browning and Browning 2016). A comparison using PCA, between the imputed and non-imputed SNP sets showed that imputation did not change the clustering of these 723 samples (Fig. S16). The 3.8 M SNPs were subsequently filtered for minimum allele frequency (MAF), linkage disequilibrium (LD pruning or filtering), and distance between polymorphisms (thinning) in different subsets of samples to obtain fourteen sets of SNPs that ranged from 59 K to 3.8 M SNPs, which were appropriate for the various downstream analysis described below (Table S11).

Population Structure and Diversity Analysis for the Panel of 672 Vietnamese Samples

SNP sets were filtered for MAF 5%, followed by LD filtering using PLINK --indep-pairwise 50 10 0.2, with further thinning if required. We ran STRUCTURE (Pritchard et al. 2000) v2.3.5 using the default admixture model parameters; each run consisted of 10,000 burn-in iterations followed by 50,000 data collection iterations. STRUCTURE was run using K = 2 for the 616 samples using SNP set 1 (163,393 SNPs). Samples with admixture components < 0.75 were classified as admixed, and the remaining samples were classified as Indica or Japonica. STRUCTURE was run varying the assumed number of genetic groups (K) from 3 to 10 with three runs per K value for the 672 Vietnamese samples (SNP set 9–80,000 SNPs); from 1 to 8 with ten runs per K value for the 426 Indica subtypes from Vietnam (SNP set 10–108,420 SNPs) and the 211 Japonica subtypes from Vietnam (SNP set 11–59,815 SNPs). The output files were visualised using the R package POPHELPER v.2.2.7 (Francis 2017) including the calculation of the number of clusters (K) using the Evanno method (Evanno et al. 2005; Zheng et al. 2012). Using the combined-merged clumpp output from POPHELPER, Indica (K = 5) and Japonica (K = 4) samples were classified into Indica I1 to I5 and Japonica J1 to J4 subpopulations using a threshold of > = 0.6, with the remaining samples being classified as mixed (Im and Jm). The principal component analysis (PCA) was performed using the R package SNPRelate v1.16.0 (Zheng et al. 2012) using method = ‘biallelic’. Nucleotide Diversity (π) was measured for each of the subpopulations with VCFtools v0.1.13 using 100-kbp windows and a step size of 10 kbp.

Determining the Effect of SNPs

The effects of all bi-allelic SNPs (low, medium and high effects) on the genome were determined based on the pre-built release 7.0 annotation from the Rice Genome Annotation Project (http://rice.plantbiology. msu.edu/) using SnpEff (Cingolani et al. 2012) release 4.3, with default parameters. The complete set of 3,750,621 SNPs (SNP set 2) which contained on average one variant every 99 bases was annotated. Using sequence ontology terms, the effect of each SNP was classified as described by SnpEff. A summary of the SNP effect analysis is available in Table S10.

Genome-Wide Association Analysis

Three independent analyses were conducted using the full panel (672 samples, 361,191 SNPs), the Indica subpanel (426 samples, 334,935 SNPs) and the Japonica subpanel (211 samples, 122,881 SNPs), SNP sets 12, 13 and 14 respectively (Table S11). The GWAS analysis was performed by employing the R package Genome Association and Prediction Integrated Tool (GAPIT) version 3.0 (Lipka et al. 2012; Tang et al. 2016). The covariate matrix was generated in STRUCTURE. We used the combined-merged output from POPHELPER for the full panel (K = 8), the Indica subpanel (K = 5) and the Japonica subpanel (K = 4). The covariate matrix and the kinship calculated in GAPIT were included in the GWAS model to control for false positives. The SUPER (Settlement of MLM Under Progressively Exclusive Relationship (Wang et al. 2014) method integrated into GAPIT, designed to increase the statistical power, was used to perform the association mapping analysis. The SUPER method was implemented in GAPIT by setting the parameter of “sangwich.top” and “sangwich.bottom” to CMLM and SUPER, respectively. A quantile-quantile (Q–Q) plot was used to check if the model was correctly accounting for both confounding variables. Associations held by peaks with -log10 (p-value) ≥ 8.0 (equivalent to a FDR < 0.01) were used to declare the significant associations.