For this study, 240 ORSC accessions from the International Rice Germplasm Collection (IRGC) previously genotyped by Kim et al. (2016), were selected for phenotypic evaluation (Supplementary Table 2). This collection included 102 ORSC accessions from Southeast Asia, 79 accessions from South Asia, 55 accessions from East Asia and four accessions from Australasia. Of these, 222 accessions were evaluated in a screenhouse at the IRRI, Los Baños, Philippines and 130 were evaluated under greenhouse conditions in the United States, with 112 accessions phenotyped in more than one location. At IRRI two plants were grown of each accession with one plant per pot and each pot was considered a single replication. Seed of each accession was planted between June 20 and July 4, 2008, and the plants were phenotyped for the 44 morphological traits related to the vegetative, reproductive and harvest growth stages (Supplementary Table 3). Twenty-two of these 44 traits were measured five times on an individual plant to improve the accuracy of the measurement (Supplementary Table 4).

In the United States, 104 accessions were phenotyped at both Cornell University (CU) in Ithaca, New York and the Dale Bumpers National Rice Research Center (DB) near Stuttgart, Arkansas. Three plants were grown of each accession with one plant per pot and each pot was considered a single replication. At CU, 104 accessions were planted between August 31 and September 8, 2008, and each plant was phenotyped for 16 morphological traits. At DB, 129 ORSC accessions were grown in the greenhouse with three plants per accession, one plant per pot, each pot was considered one replication, and plants were phenotyped for 18 morphological traits. In Arkansas, each accession was phenotyped in two different years, thus a total of six plants were characterized for most accessions. The accessions were planted August 24 to September 6, 2007 (86 accessions); August 1 to September 12, 2008 (124 accessions) and August 12 to 13, 2009 (41 accessions). Supplementary Table 2 lists the location where each accession was grown.

Supplementary Table 3 summarizes the 57 phenotypic traits characterized across the three locations and indicates the location(s) where the trait data were collected. Additional details about the collection of phenotypic data at IRRI are included in Bioversity International et al. (2007). Also listed are the corresponding Planteome acronyms and the trait ontology or crop ontology terms (Cooper et al., 2018). At IRRI, for the 22 quantitative traits with five measurements per plant, the within-plant variances were very small, thus a mean was calculated for subsequent analyses. The replicated trait data collected for each accession by location was used to calculate a best linear unbiased estimator (BLUE) for each trait and accession by location as described in the next section. The BLUEs are listed by location in Supplementary Table 4.

Using the replicated trait measurements taken at each location and a relationship matrix constructed from the GBS genotypes (Kim et al., 2016), we fit multi-trait Bayesian hierarchical models (Greenberg et al., 2011) to the data from each site separately. The model is:

Location parameters of the form y_i. are row-vectors of the corresponding parameter matrices that have data points as rows and traits as columns. The overall intercept μwas modeled with a high-variance (Iσ02 = 10⁶) Gaussian prior. The errors were modeled using a multivariate Student-t distribution with three degrees of freedom to dampen the effects of outliers (Greenberg et al., 2010; Greenberg et al., 2011). We computed the parameters using a C++ program based on an openly-available library³ that implements methods described by Greenberg et al. (2011).

The marker effects enter through the eigenvectors U of the relationship matrix, weighted by square roots of their eigenvalues. We estimated the relationship matrix from all non-singleton SNPs in each data set using the van Raden method (VanRaden, 2008). We used all the eigenvectors that correspond to non-zero eigenvalues. If this regression is performed using a Gaussian prior on the coefficients γ_j., it is the Bayesian analog of the (matrix-variate) mixed effect model (Kang et al., 2010; Hoffman, 2013). The accession means μj.a⁢c⁢care analogous to BLUEs. Modes of their posterior distributions were used in the subsequent analyses.

The covariance matrices Σ^x were modeled using Wishart distributions with weakly-informative Wishart priors with two degrees of freedom for Σ^e and Σ^s and four degrees of freedom for Σ^a (Gelman et al., 2004; Greenberg et al., 2011). This implies a uniform prior on narrow-sense (marker) heritability. Estimated model parameters can then be used to calculate marker heritability:

where N_acc is the number of accessions, and the rest of the notation is as described above. We used the variance of genome-estimated breeding values (GEBV).

rather than values from the marker covariance matrix Σp,pabecause the latter reflects the additive covariance only when the prior on the principal component regression (see above) is Gaussian. The Σp,ps values estimate the portion of total genetic variance not explained by genotyped markers. They include the non-additive effects and additive effects not tagged by SNPs. The sum of the GEBV and background variance, divided by total phenotypic variance, is broad-sense heritability (Lynch and Walsh, 1998).

Using estimates of phenotypic values for each accession, we fit a Gaussian mixture model (Robert, 1996) using a variational Bayes approach (McGrory and Titterington, 2007). The model is

where the notation is the same as for the hierarchical model above and μ_m⋅ is a row vector of group means. Each group has a separate mean and covariance matrix. The prior population size α is set to a small value, shrinking groups with few members to zero. The number of groups must be set a priori.

We perform inference by running the model multiple times from different starting points and picking the best fit, as assessed using the deviance information criterion (DIC, McGrory and Titterington, 2007). Small DIC values indicate better fit. Our implementation of the model is available as the R package MuGaMix⁴. We use two methods to assess how many groups are statistically supported: the DIC and the number of groups with at least two accessions. We include the pipeline (implemented in R) that runs the mixture model, processes the results, and generates all the plots and summary tables presented in this manuscript in the Supplementary Material File 1.

For population structure analysis, the SNP datasets from Kim et al. (2016) and Huang et al. (2012) were merged into a single dataset. First, the Kim et al. (2016) SNPs from ORSC and O. sativa accessions were imputed using Beagle 5.0 (Browning et al., 2018). Next the intersection between the SNPs in the imputed Kim et al. (2016) and the (already imputed) SNPs in the Huang et al. (2012) was determined and a merged file was generated using BCFtools version 1.14 (Danecek et al., 2021). Population structure analysis was done using fastSTRUCTURE version 1.0 (Raj et al., 2014) with K values ranging from 2 to 9.

The local ancestry of the ORSC individuals was determined using RFMix (Maples et al., 2013). The training panel included ten individuals from each subpopulation (W1-W6) with the highest global ancestry as reported by Kim et al. (2016) and included individuals from each of the five O. sativa subpopulations, indica (IND) (10 individuals), aus, (AUS) (9 individuals), temperate japonica (TEJ) (9 individuals), tropical japonica (TRJ) (10 individuals), and aromatic (ARO) (6 individuals). Of the 60 ORSC individuals in the training panel, 49 were phenotyped in the current study. Prior to RFMix analysis, the genotypic data were imputed and phased using Beagle 5.0 (Browning et al., 2018).

Broad-sense and marker heritability distributions for each trait and location are presented in Supplementary Table 5. Heritabilities range from 0.98 to 0.33 and are generally moderate to high. Markers account for at least 33% of total genetic variance across all traits and locations, but it is the generally high H² that gives confidence in within-site reproducibility of our phenotypic measurements.

Fourteen traits were measured at more than one site. We tested among-site reproducibility of these measurements by estimating correlations among phenotypes evaluated on the same accessions and generally see positive correlations (Supplementary Figure 2 and Supplementary Table 6), again suggesting that phenotypic measurements are reproducible across sites despite the diverse climatic conditions. One notable exception is panicle number (PNNB). Values of this trait are positively correlated between Cornell and Dale Bumpers (temperate and sub-tropical sites, respectively), but both are negatively correlated with IRRI measurements (tropical site), suggesting a genotype by environment interaction.

Using estimates of 32 phenotypic values (binary traits were excluded) for each accession in the IRRI data set (n = 222), we fit a Gaussian mixture model (Robert, 1996) using a variational Bayes approach (McGrory and Titterington, 2007). To infer the underlying number of groups that is statistically supported and biologically meaningful, we ran analyses assuming two to 15 groups a priori and assessed model fit (Supplementary Figure 3) using the deviance information criterion and the number of non-empty groups. Both metrics lend statistical support for the division of our accessions into as many as 10 phenotypic groups. However, to make biologically relevant inferences we need a parsimonious classification that is not only statistically significant, but also robust to model uncertainty, reproducible, and interpretable given other biological information. The variational Bayes approach we used is fast but produces only point estimates of parameters. To test robustness of these estimates to initial conditions, we started the model fitting process from five independent sets of initial values. Model fit metrics were consistent across initial value sets (Supplementary Figure 3).

Moving beyond model fit estimates, we tested stability of group membership across model specifications and starting value sets. To anchor these inferences, we started by noting taxonomic classification of our accessions, as found in the IRRI-GRIN-Global database. Based on traditional classification methods, accessions in this database are listed as either O. rufipogon, O. nivara, or, if the taxonomic assignment is uncertain, as Oryza spp. or Oryza hybr. We combined the latter two categories into one, calling it Oryza spp. We tracked each accession, starting from the species designations, across phenotypic groups inferred given different a priori group numbers (N_G) using a Sankey plot (Supplementary Figure 4; Kennedy and Sankey, 1898; implemented in the ggsankey R package)⁵. A line, which is like a ribbon in such plots, represents a single accession, while boxes arranged along the y-axis are groups. Groups derived from models assuming a different number of groups are shown along the x-axis. We can thus follow a given accession as its line connects to each group across models with different N_G. If we set N_G = 3, the model infers one large group comprising most accessions. The second largest group is composed predominantly of accessions classified as Oryza spp. Increasing the a priori N_G to four, we can assign most O. rufipogon and O. nivara accessions to separate groups. If we further increase N_G, distinct groups corresponding to the traditional species designations are maintained, with an additional group that comprises most Oryza spp. accessions. However, additional small groups inferred with higher a priori N_G do not appear to maintain coherent membership across model runs (Supplementary Figure 4). Thus, while there is statistical support for finer classification, these groupings are not robust to model specification. We therefore focused our analyses on inferences assuming four a priori groups.

Given that models with N_G = 4 appear to be the most parsimonious, while still capturing biologically meaningful stratification of the ORSC, we next tested these models for robustness. We started model fitting from 20 sets of initial conditions and averaged posterior probabilities that a given accession belongs to a particular group across runs. We see a set of accessions that reproducibly belong to the P1 or P4 groups (Figure 1A), corresponding to O. rufipogon and O. nivara, respectively. The P3 group has the least reproducible membership. This is because for some initial conditions, a small number of accessions are assigned to it, while for others these lines are placed in P2 (see Supplementary Material 1 for a detailed analysis of these shifts). This may be due to the approximate nature of the variational Bayes inference. Therefore, we added the posterior probabilities of P2 and P3 assignment in subsequent analyses and called the resulting group P2/P3.

The resulting P1 phenotype group contains 106 accessions, of which 82 (78%) are classified as perennial O. rufipogon species. In contrast, P4 has 73 accessions with 43 (59%) classified as O. nivara. About half of the P2/P3 group are designated as Oryza spp. in the IRGC, with the remaining accessions about equally divided between O. nivara and O. rufipogon (Supplementary Figure 5 and Supplementary Table 2).

To further test the reproducibility of phenotypic group inference, we used data collected at DB and CU. These data sets include different, but overlapping, phenotype sets that were measured on fewer than half the accessions represented in the IRRI data (Supplementary Tables 2, 3). Despite this, we can still identify a similar set of three groups (Figure 1B). This time, setting N_G = 3 a priori is sufficient to separate O. rufipogon and O. nivara from each other and from a group that largely corresponds to the Oryza spp. designation from IRRI-GRIN-Global. Taking these observations together with the above statistical robustness analyses, we see that a core set of accessions can be reproducibly assigned to a phenotypic group, while inference for others is less certain.

Having established a robust phenotype-driven grouping of ORSC accessions, we next wanted to define a minimal combination of phenotypic traits that is sufficient to reproduce these divisions. Our mixture model allows both phenotype means and covariance to vary among groups. We started by ranking traits according to how well their values correlate with our groups (Figure 2A). We used a multivariate regression with group identity probabilities as the response variable and trait values as predictors. Since the last group membership is fully determined by the other three, only P1 and the combined P2 and P3 membership probabilities were used. We estimated regression coefficients and their covariance and used them to calculate a Hotelling T² statistic (Hotelling, 1931) of each trait-group association. This is a multivariate version of the square of Student’s t statistic widely used to determine statistical significance of individual regression coefficients. However, in this case we were not interested in the significance of associations. Instead, we ranked traits by their T² values (Figure 2A). Two traits, culm length (CULT) and number of empty spikelets per panicle (UNFILLED), stand out for their association with phenotypic groups. Their Hotelling T² statistics (0.71 for CULT and 0.70 for UNFILLED) are over twice the value of the next seven highest trait-group associations for which T² ranges from 0.17 to 0.28 (Table 1 and Supplementary Table 7). Accessions belonging to the P1 group are on average taller, with P4 the shortest and P2/P3 intermediate. P2/P3 accessions typically have the most empty spikelets per panicle, while P4 accessions have the least (Supplementary Figure 6).

To find pairs of traits that change correlations across groups, we first estimated within-group associations and then used a permutation test (randomly assigning accessions to groups) to assess statistical evidence for across-group variation (Figure 2B). The most significant correlation changes occur between the number of empty spikelets (UNFILLED) and either panicle length (PNLG) or culm diameter (CUDI). These correlations are positive (0.54 and 0.48, respectively) only in P4 (Supplementary Table 7).

To test if subsets of traits can reproduce our phenotypic groups, we took 11 traits that belong to pairs that change correlations with probability at least 0.001 (Table 1). This set also includes the top two traits (CULT and UNFILLED) whose values correlate best with group assignments. To expand the trait subset, we added seven phenotypes with T² > 0.17. Hotelling values showed a marked drop below this cut-off. These two sets have four traits in common [CULT, UNFILLED, CUDI, and anther length (ANTLT)]. We also used the union of the two sets, containing 16 traits, or half the number in the full data set. Although covariance matrices include both correlations and variances, all traits with appreciable among-group variance changes are already accounted for by using correlation and mean differences (Supplementary Figures 6, 7).

We used three data subsets (nine traits associated with phenotypic groups by value only, 11 that switch correlations, and 16 that were associated either by value or correlation) to re-fit our mixture model, setting the a priori number of groups at three. We ran our analyses 20 times and averaged posterior probabilities of accessions belonging to each group. Using all three data sets, we consistently see three groups that separate O. rufipogon from O. nivara. The nine- and 16-trait data sets assign most P2/P3 accessions to P1 (Figure 2C). The 11 traits most likely to change correlations across groups appear to be sufficient to recapitulate the grouping inferred from the whole data set. The two traits, CULT and UNFILLED, most closely aligned with groups by value and are included in this collection.

Our approach identified a minimally sufficient set of traits that captures most of the phenotypic group structure in the IRRI data. This does not mean that other sets of traits could not achieve similar results but an exhaustive search for all possible combinations that could achieve the same result was not undertaken considering time and the computational expense. Instead, we pursued two approaches using the Cornell and Dale Bumpers data sets to find alternative sets of traits that could be used to re-infer the three major phenotypic groups identified using the full data set. First, we used the process outlined above to identify traits associated with the phenotypic groups that were inferred from IRRI data by value or that change correlations with other traits. Supplementary Figures 8, 9 visualize the Dale Bumpers data analyses and Supplementary Figures 10, 11 visualize the Cornell data analyses with the actual values in Supplementary Table 7. This approach did not identify additional trait subsets that could reliably identify phenotypic groups (Supplementary Figures 12A,C). Second, we identified traits evaluated at Dale Bumpers and Cornell that are the same or similar to the 11 phenotypes (Supplementary Table 3) that defined the grouping patterns in the IRRI data. Re-running our mixture model with these analogous traits reproduced the original groups, albeit at lower resolution (Supplementary Figures 12B,D). Using fewer than 11 traits thus appears to significantly degrade group assignment precision and reproducibility.

While our phenotypic groups largely correspond to IRRI-GRIN-Global species designations, we wanted to know how the groups related to previously characterized genotypic subpopulations. The accessions used in this study were genotyped by Kim et al. (2016) who reported six ORSC subpopulations (W1-W6). Based on RFMix (Maples et al., 2013) analysis, we identified introgressions among ORSC subpopulations and between the ORSC and five O. sativa subpopulations (IND, AUS, TRJ, TEJ, ARO). We partitioned accessions with more than half of their genome introgressed from O. sativa into a separate “admixture with O. sativa” group (ADM/OSAT). By means of a Sankey plot, we related the species designation, phenotypic group and genetic subpopulation for the ORSC accessions with complete information (Figure 3A and Supplementary Table 2). Across the ORSC subpopulations, the majority of the W1, all W3, and about half of the W6 accessions (as per Kim et al., 2016) clustered in the P1 group, taxonomically designated as O. rufipogon. Conversely, all of the W5 accessions and most of the W2 clustered in the P4 group, taxonomically designated as O. nivara. Finally, accessions with significant O. sativa introgression fell almost exclusively into the P2/P3 group and are classified as Oryza spp. The P2/P3 group identified using the reduced set of 11 traits also includes most of the W6 subpopulation (Supplementary Figure 13). Eight accessions taxonomically classified as O. rufipogon in the IRGC were genetically designated as W2 or W4 by Kim et al. (2016) and clustered with the P4 group in this study. These accessions had less than 20% of their genome from O. sativa and were re-classified as O. nivara in our analyses (including in the results reported above) and are identified in Supplementary Table 2.

Surveying the individual introgressions across the genome, there were no obvious patterns identifying particular genomic regions underlying the phenotypic group differentiation (Figure 3B). This result was not surprising, given that these groups are distinguished by multiple polygenic traits. We do see that those accessions from the P2/P3 group were more likely to harbor O. sativa introgressions (Figures 3B,C) and the median portion of the genome coming from O. sativa is also larger in P2/P3 than in P1 or P4 accessions (Supplementary Figures 14B, 15). However, P1 and to a lesser extent P4 harbor a few accessions with significant contributions from O. sativa or ORSC subpopulations typically associated with a different phenotypic group, suggesting that similar phenotypic syndromes can be achieved through divergent genotypes.

Most O. sativa introgressions come from the indica subpopulation (Figure 3 and Supplementary Figure 14B). Only accessions belonging to the P2/P3 group show appreciable contributions from O. sativa japonica. Introgressions from aus are relatively more prevalent in accessions from both P2/P3 and P4.

To compare the population structure of the collection of ORSC analyzed by Kim et al. (2016) (n = 286) with that generated by Huang et al. (2012) and subsequently analyzed by Wang et al. (2017) (n = 435), a merged SNP dataset was constructed and analyzed using fastSTRUCTURE (Raj et al., 2014). The combined dataset contained 55,213 SNPs shared between them. Structure analysis of the merged data resulted in six subpopulations that match nearly perfectly with the six subpopulations reported by Wang et al. (2017), (Figure 4A). The six subpopulations, designated W1 to W6 by Kim et al. (2016) have a clear connection to those designated Or-A to Or-F by Wang et al. (2017) as illustrated by shared merged groups. The main differences are due to the thresholds used to classify accessions designated as admixed, and the different size and composition of the two collections. Where Wang et al. (2017) used a threshold of > 80% ancestry to classify accessions as belonging to one of the wild groups, Kim et al. (2016) used a threshold of >60%. As a result, roughly half of the accessions designated as W1 by Kim et al. (2016) were identified as admixed in the combined analysis. It is also noteworthy that Kim et al. (2016) had a smaller proportion of wild accessions from China (W6), and a larger proportion from Papua, New Guinea (W3) and from Nepal (W5), which impacted the subpopulation structure that emerged in the analysis. Nonetheless, there is a clear correspondence between genetic groups in the two studies, with W1 and W3 (predominantly O. rufipogon as described by Kim et al., 2016) corresponding to Or-A (Wang et al., 2017), W6 (Chinese O. rufipogon) to Or-B, and W2 (predominantly O. nivara) to Or-C. W4 (aus-like wild ancestor as described by Kim et al., 2016) is split between Or-D (wild) and Or-E (feral), and interestingly, the accessions clustered as Or-E by Wang et al. (2017) were recognized as a separate group by Kim et al. (2016) when a higher K value (K = 8) was used in the analysis.

The groups identified by Wang et al. (2017) also show reasonable correspondence to the original three groups proposed by Huang et al. (2012), (Figure 4B), though Huang et al. (2012) did not explicitly call out the large proportion of ORSC accessions that are highly admixed with O. sativa as both Wang et al. (2017) and Kim et al. (2016) did. Given the different compositions of the two germplasm collections, this merged genomic analysis allows us to interpret our findings about the three phenotypic groups in a larger context.

There is support for geographic subpopulation structure whereby accessions collected from the Malay Archipelago (Malaysia, Indonesia, Philippines, Papua New Guinea), which represents the southeastern extreme of the geographic distribution of the ORSC, are phenotypically classified as P1 in this study; they are also predominantly classified as O. rufipogon in the IRGC (Supplementary Figure 16). Accessions from India, Nepal, and Sri Lanka, representing the northwestern extreme of the geographic distribution, are predominantly classified as P4 in this study and O. nivara in the IRGC. Accessions collected from China were genetically classified as O. rufipogon with substantial admixture from O. sativa (Kim et al., 2016), but they were phenotypically more similar to annual O. nivara in this study. Accessions collected across the coastal regions of South and Southeast Asia are recognized phenotypically as P2/P3 and are typically the most highly admixed with O. sativa. A great diversity of phenotypic variation is found in mainland Southeast Asia where ORSC accessions manifesting both P1 and P4 phenotypic syndromes co-exist, often sympatrically or in contiguous environments. These geographical distributions of O. rufipogon and O. nivara are clearly shown in Supplementary Figure 16B where only the locations of the “core” accessions, those for which the species, phenotypic group and genotypic subpopulation agreed, are shown. The fact that cultivated forms of rice, predominantly representing the indica subpopulation, overlap with wild rice habitats in many parts of Asia creates endless possibilities for cross fertilization and the emergence of novel phenotypes. In addition to their inherent value and interest, these naturally occurring interspecific populations represent valuable, ecologically viable pools of variation for plant breeders interested in identifying novel forms of climate resilience for crop improvement.

The traits measured in this study emphasize highly heritable, easily measured morphological and physiological traits that can be readily observed in populations of crop wild relatives grown under controlled conditions. Many of the same traits are evaluated by genebanks and widely used for identification, tracking and management of ex situ germplasm collections. Of the traits measured here, four were significantly associated with phenotypic groups and with life history; these included days to heading, plant height, percent filled spikelets per panicle (seeds/panicle), and spikelet/seed size (related to grain weight). These traits provided information that was valuable for differentiating annual and perennial life habits and were similar to traits evaluated by Grillo et al. (2009) in two O. rufipogon/O. nivara mapping populations. We also measured several traits that were potentially informative about mating habit, including anther length, style length and stigma length, but only anther length discriminated the groups based on both Hotelling T² value and correlation estimates (Table 1). These traits related to life history and mating habit support our understanding of the biological significance of our phenotypic groups.

The single most informative trait for distinguishing ORSC phenotypic groups based on the Hotelling T² values across all three datasets was seed width (HULGRWD) (T² value = 1.15) measured at Dale Bumpers, which served as a proxy for spikelet width (SPKWD) measured at IRRI, and seed volume, (T² value = 1.70) calculated from seed width and length. Spikelet width and seed width were highly correlated (r = 0.62) and both width and volume were informative for distinguishing groups. Seed traits were not measured at Cornell and therefore could not be used for further confirmation.

Comparing this analysis to the principle component analysis (PCA) of phenotypic data collected at IRRI on a different set of 116 ORSC accessions, Banaticla-Hilario et al. (2013a) also reported SPKWD as one of the most important characters separating O. nivara and O. rufipogon accessions. The other important traits included anther length (ANTLT), days to heading (DTHD) and culm length (CULT). Of note, spikelet fertility was not included in the PCA because it was highly correlated with anther length and spikelet width, leading to the conclusion that anther length and spikelet width serve as a proxy for fertility when plants are grown under controlled screenhouse conditions. Our data on DTHD collected at IRRI differentiated P1 from the other groups. The trend is the same in the Cornell and the Dale Bumpers data, but with more overlap between groups. This may be a consequence of evaluating DTHD during the short days of fall at higher latitudes. Across our three environments, unfilled grain number (UNFILGRNB) in the IRRI and Dale Bumpers data and seed weight (Seed_S2) collected from bagged and unbagged panicles at Cornell provided useful proxies for spikelet fertility. Lastly, culm length (CULT) which was measured from the soil surface to the panicle base at IRRI and plant height (PTHT) measured from the soil surface to the tip (end) of the panicle at Dale Bumpers and Cornell were prominent differentiators across all three studies. This comparison confirms that many traits distinguishing the phenotypic groups can be identified across different environments and panels of ORSC accessions. These observations are confirmed by earlier studies summarized in Supplementary Table 1 in which O. nivara was characterized by having wider seeds, higher seed production (fertility), shorter culm length (plant height), and earlier flowering (due to photoperiod insensitivity) compared to O. rufipogon which was characterized by relatively narrow seeds, low seed production, longer culm length (taller plant) and later flowering (due to photoperiod sensitivity) (Morishima et al., 1961; Oka and Morishima, 1967; Barbier, 1989a).

Considering the high information content of SPKWD (spikelet/seed width) as a distinguishing feature of the phenotypic groups identified in this study and its utility for differentiating the O. rufipogon and O. nivara species, we recommend including seed width and length in the suite of traits collected on ORSC accessions in genebanks around the world. With the availability of scanners and appropriate software, this would be an inexpensive investment that would help improve the classification of ORSC accessions, making collections of these crop wild relatives more valuable to users.

Of the ORSC accessions phenotyped in this study, approximately 41% (∼98 accessions, including the 49 extracted for use as training sets) had less than 5% introgression from O. sativa and can be considered non-admixed and wild. Most of these (n = 91) clustered in either P1 (n = 48) where most were classified as O. rufipogon or in P4 (n = 43) where most were classified as O. nivara, while a few (n = 7) clustered in P2/P3.

Based on analyses of both phenotypic and genotypic data using different model parameters, we identified 42 “core” accessions that are consistently classified as perennial O. rufipogon and 25 “core” accessions consistently classified as annual O. nivara. These “core” accessions all have <20% admixture from O. sativa and intersect with the group carrying <5% O. sativa introgression. This collection of 67 accessions provides a useful subset of ORSC materials for studying the genetic basis of the phenotypic divergence that distinguishes the two ecotypes of wild rice in Asia.

ORSC accessions that are highly admixed with O. sativa are also of interest. In this study, 25% of accessions were classified as admixed because they had >40% of the genome introgressed from O. sativa. In the previous study by Wang et al. (2017) using a different panel of ORSC accessions and a different threshold for determining admixture, 42% of wild rice samples were reported to be substantially admixed with >20% of the genome showing introgression from O. sativa. In both studies, indica and aus introgressions accounted for the largest proportion, and japonica introgressions were rare. Admixed accessions found in P1 and P4 in this study had an average of 19.8% and 19.0% of their genomes comprised of introgressions from O. sativa, respectively, while levels of admixture in P2/P3 averaged 52.6%; indeed, two P2/P3 accessions carried ∼85% O. sativa DNA (Figure 3C and Supplementary Figure 15). Nonetheless, the genome-wide level of admixture alone, is not a good predictor of phenotype, given that 13 accessions classified phenotypically into P1 carry > 40% O. sativa DNA, and 9 accessions classified phenotypically as P4 carry > 40% O. sativa DNA. We therefore infer that it is the particular distribution of O. sativa introgressions across the genome and the subpopulation ancestry of each introgression, in combination with the ORSC genetic background that determines the phenotypic outcome. Our evaluation of this ORSC collection provides the foundation for future research to examine the genotype-phenotype relationships for specific traits of interest.

In rice, SINE codes have been associated with life history (annual, intermediate or perennial growth habit) and found useful for distinguishing O. nivara from O. rufipogon (Motohashi et al., 1997; Cheng et al., 2003; Chen et al., 2004; Xu et al., 2007). Kim (2016) developed a SINE code and used it to classify most of the ORSC accessions used in this study (Supplementary Table 2). Using the SINE code, 74.5% of the accessions classified as P1 were identified as perennial (compared to 72.5% of accessions classified as O. rufipogon), and 62.5% of accessions classified as P4 were identified as annual or intermediate (compared to 61.7% classified as O. nivara). In P2/P3, 47.4% of accessions were classified as perennial and 52.6% as annual or intermediate using the SINE code (compared to 23.3% classified as O. rufipogon, 25.6% as O. nivara in P2/P3 with the remainder identified as Oryza spp.). Thus, the SINE-based classifications approximated the species designations in these groups, though they were not entirely concordant.

Kim et al. (2016) further analyzed 286 ORSC accessions for variation at 25 polymorphic sites in the chloroplast genome and identified unique chloroplast haplotypes associated with three wild subpopulation groups located at the geographical extremes of the species range. These may be of interest for investigating refuge populations of ORSC and elucidating phylogenetic relationships among wild populations of Oryza but chloroplast haplotype analysis alone was not robust enough to be useful as a predictor of genetic subpopulation, life habit or species in the ORSC. Nevertheless, the chloroplast data provided clear evidence of gene flow between annual and perennial populations, as well as between ORSC and O. sativa populations, and confirmed that it occurs through both pollen dissemination and seed dispersal (likely facilitated by human migration). Taken together, these data suggest that the population structure of the ORSC is the result of complex evolutionary pathways that intersect and loop back upon each other due to highly permeable ‘species’ boundaries and, as documented in this study, highly plastic phenotypic variation. Information about SINE-code and chloroplast haplotype variation (as determined by Kim et al., 2016) corresponding to each accession is summarized in Supplementary Table 2. These features, along with geographical and ecological information about where each accession was collected, disease and insect resistance, grain quality and additional use-data provided by the IRGC database add value to these genetic resources and enhance our understanding of where they come from and how they might be used in the future.

The P2/P3 group is characterized by 43 accessions with trait combinations that form a coherent group, despite the fact that they do not match either of the a priori taxonomic species descriptions. This group is comprised of 29 (67%) admixed ORSC accessions with > 40% of the genome derived from O. sativa but it also contains accessions classified as annual O. nivara or perennial O. rufipogon. The fact that the P2/P3 accessions were collected in diverse ecological and geographic regions, especially from coastal regions of China, Southeast and South Asia (Supplementary Figure 16), and yet share a common suite of emergent phenotypes suggests that particular combinations of traits may evolve repeatedly and independently when annual, perennial and intermediate forms of common wild rice come together as sympatric swarms in unsupervised settings. The patterns of admixture observed in this group suggest that diverse ORSC populations have hybridized among themselves and with local populations of O. sativa, giving rise to forms of trait variation not observed in the parental populations. In rice, several lines of evidence suggest that annual and intermediate forms have emerged repeatedly from genetically and geographically diverse populations of perennial ancestors in response to changing patterns of temperature, rainfall and CO₂ (Grillo et al., 2009; Banaticla-Hilario et al., 2013a,b; Liu et al., 2015). The ORSC populations examined here provide opportunities to deepen our understanding of incipient speciation and to identify mechanisms by which annual life forms may evolve from perennial ancestors in response to changes in the environment. They also provide material for examining what kinds of selection pressure(s) are associated with a shift toward the annual, early-flowering, seed-bearing habit or conversely, serve to maintain the perennial, stoloniferous, late-flowering and vegetatively vigorous habit. During times of climate change, a strategy of balancing selection to maintain components of annual, intermediate and perennial life forms may be favored by evolution to allow for the emergence of new resilience mechanisms. We have a lot to learn from the existence of these dynamic ORSC populations. They have survived many waves of climate change in the past and are likely to hold secrets about survival strategies for the future.