Bulliform Phytolith Size of Rice and Its Correlation With Hydrothermal Environment: A Preliminary Morphological Study on Species in Southern China

In the last decade, our understanding of rice domestication has improved by new archaeological findings using advanced analytical techniques such as morphological and morphometric analyses on rice grains, spikelet bases and phytoliths, and ancient DNA analysis on rice remains. Previous studies have considered the size of rice bulliform phytoliths as a proxy for tracking the domestication process. These phytoliths are often abundant and well preserved in sediments, and their shape is under the control of numerous genes, which may shift toward larger sizes by genetic mutation in domestication. Therefore, it has been assumed that the bulliforms of domesticated rice are usually larger than those of wild ones; however, morphometric data supporting this assumption are lacking in the literature, thereby requiring additional evidence to test its veracity. In this study, the vertical and horizonal lengths of bulliform phytoliths were measured in four rice species (domesticated Oryza sativa and wild Oryza rufipogon, Oryza officinalis, and Oryza meyeriana) from different regions of southern China. We found that the bulliform morphometric data of wild and domesticated rice overlapped and that there was no statistically significant difference between them. Therefore, bulliform size could not be used as a diagnostic indicator to distinguish domesticated rice from wild species and is a supporting rather than conclusive proxy for determining the domesticated status of rice in archaeological research. We further found that larger rice bulliform sizes likely occurred at the locations with higher temperature, precipitation, and water levels, indicating hydrothermal environment is an alternative factor influencing the size of rice bulliform phytoliths. For further archaeological use of an increasing size trend of bulliform phytoliths to reveal the process of rice domestication, we present some suggestions for controlling the influence of hydrothermal factors. Even so, the combination of bulliform phytolith size with other established criteria is strongly suggested to provide precise identification of wild and domesticated rice in future research.


INTRODUCTION
Asian rice (Oryza sativa L.) is one of the most important crops and forms a staple food for more than half of the global population (Nayar, 2014). Understanding its origins and domestication from wild rice (Oryza rufipogon Griff.) is thus an important aim for researchers. There are two major subspecies of domesticated O. sativa, Oryza sativa indica, which is thought to have originated in the Himalayan region, and Oryza sativa japonica, which is thought to have originated in China (Londo et al., 2006). Oryza rufipogon is generally recognized as the ancestor of Oryza sativa , with two distinct domestication events leading to the two subspecies of domesticated rice (Londo et al., 2006). However, this hypothesis is still being vigorously debated as other evidence supports a single origin of Asian rice (Huang et al., 2012). There are more than 20 wild rice species recognized in the genus Oryza (Nayar, 2014;Stein et al., 2018), which belongs to the family: Poaceae and tribe: Oryzeae. Three of these species are found in China, Oryza rufipogon, Oryza officinalis, and Oryza meyeriana (Fan et al., 2000).
In the last decade, scientific understanding of rice domestication has greatly improved by new archaeological and genetic evidence using advanced analytical techniques such as flotation combined with morphometric analysis on rice grains and spikelet bases (Liu et al., 2007;Gross and Zhao, 2014;Zheng et al., 2016), phytolith analysis (Wu et al., 2014;Huan et al., 2015;Zuo et al., 2017), pan-genome analysis (Wang et al., 2018a), and genome-wide association studies (Huang et al., 2012;Civáň and Brown, 2017;Choi and Purugganan, 2018). According to recent findings, a general consensus among scholars has been reached that rice was first domesticated in the middle and lower Yangtze River regions of southern China (Molina et al., 2011;Gross and Zhao, 2014;Silva et al., 2015;Choi et al., 2017;Zuo et al., 2017), although some consider the Pearl River region to also be a part of the original area where rice domestication occurred (Huang et al., 2012;Wei et al., 2012).
The estimated time of origin of rice domestication is before 13,000 BP, based on molecular clock analysis (Molina et al., 2011;Choi et al., 2017;Choi and Purugganan, 2018), which is much older than the earliest archaeological date of domestication (<10,000 BP) Larson et al., 2014). This disagreement may result from genetic studies that mainly focused on identifying the origins of alleles associated with domestication (e.g., sh4, rc, laba1, prog1), which likely emerged in wild rice prior to domestication (Choi et al., 2017;Civáň and Brown, 2017). On the other hand, archaeological studies attempted to detect the first appearance of morphological traits associated with domestication in archaeobotanical remains Jones and Liu, 2009). The molecular and archaeological chronologies may date two main phases in the macroevolutionary process: the emergence of a trait and the success of that trait (the trait becomes quantitatively significant within a population) (Katz, 2019). However, the most notable chronological dispute over the rice domestication is between two archaeological opinions: one suggests that the process of rice domestication may have begun around 10,000-9,000 BP (Liu et al., 2007;Zheng et al., 2007;Wu et al., 2014;Zheng et al., 2016;Zuo et al., 2017), while the other suggests that domestication of rice did not occur until around 8,000-6,000 BP Fuller et al., 2008;Fuller et al., 2010;Larson et al., 2014). This debate is largely attributable to the differences in the methods employed, and the criteria used by various authors to identify domestication in rice remain from early sites in the Yangtze River region such as Shangshan, Kuahuqiao, and Xiaohuangshan. Establishing accurate and feasible criteria for distinguishing between domesticated and wild rice is thus of prime importance.
Three important lines of archaeological evidence have frequently been used in China, including grain size and morphological characteristics, spikelet bases, and phytoliths (Liu and Chen, 2012;Fuller and Castillo, 2014;Fuller, 2018). Grain morphometrics are considered to be semi-domestication traits and may not be diagnostic indicators of early domestication (Fuller, 2007;, mostly due to the considerable variation and overlap in length between domesticated and wild populations which leads to some proportion of false assignments in ancient rice grains (Fuller, 2007). Many studies have focused on the form of spikelet base, which is thought to be the most diagnostic trait in rice remains in terms of identifying domestication status (Fuller, 2018). However, these are insufficient, because the distinctive characteristics of immaturity, shattering, and nonshattering states and/or wild, japonica, and indica rice based on spikelet bases are divergent among the criteria provided by different researchers (Zheng et al., 2007;Pan, 2008;Fuller et al., 2010;Pan, 2011;Gross and Zhao, 2014;Zheng et al., 2016), and their diagnostic power for domestication has the potential to be more reliable. More importantly, these macrobotanical remains do not preserve well in early sediments with acid soil; therefore, very few have been recovered from sites dated earlier than 9,000 BP (Zhao, 2011;Qin, 2012;Zhao and Jiang, 2016).
Phytoliths have played an important role in the identification of rice remains recovered from early archaeological sites, due to their high resistance to decomposition (Piperno, 2006;Ball et al., 2016). Double-peaked phytoliths from husks and bulliform phytoliths from rice leaves are both certainly diagnostic indicators of Oryza and show variation within and between species. A number of identification criteria based on these phytoliths have been suggested and widely used in the last 20 years (Zhao et al., 1998;Lu et al., 2002;Gu et al., 2013;Wu et al., 2014;Huan et al., 2015;Hilbert et al., 2017;Zuo et al., 2017;Wang et al., 2018b). Although the utility of these methods for distinguishing domesticated from wild rice is still under debate (Fuller et al., 2010;Fuller, 2018), they are recognized as key alternative methodologies besides using morphological domestication data of rice macroremains. Three-dimensional measurements and discriminant function analysis of double-peaked phytoliths are useful in determining the wild/domesticated nature of rice remains; however, double-peaked phytoliths usually present their side and top view under the microscope, which does not meet the requirements of morphometric analysis (present in front view) (Zhao et al., 1998;Gu, 2009), making the work arduous in most cases.
Rice bulliform phytoliths are abundant in rice leaves and are often well preserved and represented in archaeological sediments (Fujiwara, 1993;Wang and Lu, 1993). Generally, bulliform phytoliths in Oryza have a distinctive fan shape with numerous scale-like decorations on the half round side (lateral side) (Lu et al., 2002;Wang and Lu, 2012;Gu et al., 2013). Morphological measurements and number of scale-like decorations along the scalloped edge have been employed to distinguish wild Oryza species from domesticated ones. Studies on modern rice plants and paddy surface soils have suggested that bulliform phytoliths with ≥9 scale-like decorations were likely domesticated, while those with <9 were generally wild (Lu et al., 2002;Huan et al., 2015). Whether this feature is a useful domestication indicator remains inconclusive and requires further validation; in addition, genetic explanatory mechanisms of bulliform scale-like decoration variation between species remain unclear and deserve further study.
Bulliform shape of rice appears to be under the control of 16 genes (QTLs) (Zheng et al., 2003a) and phytoliths may shift toward larger sizes as a result of genetic mutation during the domestication process (Zheng et al., 2003b;Piperno, 2006;Luo et al., 2016). Some researchers have, therefore, assumed that the bulliforms of domesticated rice are usually larger than wild ones and that the trend of increasing size in rice bulliform phytoliths could reflect domestication of rice (Zheng et al., 2003b, Zheng et al., 2004. In recent years, vertical and horizontal lengths (i.e., sizes) of rice bulliforms have increasingly been used as a proxy for tracking the domestication process or determining the degree of domestication at different sites, such as the Tanghu , Zhuzhai (Wang et al., 2018b), Shunshanji (Luo et al., 2016), Shangshan, Hehuanshan, and Huxi sites (Zuo et al., 2017;Qiu et al., 2019). However, to date, morphometric data from modern rice plants supporting this method and their assumptions are missing (Pearsall et al., 1995;Zhang and Wang, 1998;Ma and Fang, 2007;Gu et al., 2013), and thus additional evidence is required to test its veracity.
Moreover, some studies have argued that changes in phytolith size may not only be triggered by domestication but also influenced by environmental factors, such as CO 2 concentrations (Ge et al., 2010), evapotranspiration rates (Issaharou-Matchi et al., 2016), and water levels in the growing habitat (Fuller, 2018). Fuller (2018) indicated that the 16 genes suggested by Zheng et al. (2003a) only explained between 37 and 54% of bulliform variation, suggesting that the environment or growing conditions also play an essential role. These authors further speculated that if shifting bulliform morphology was merely a phenotypic response to environmental conditions, it would be a less useful indicator of domestication. Thus, without data supporting the exclusion of environmental factors, bulliform phytolith measurements alone may not be an accurate identification tool for distinguishing between domesticated and wild Oryza species.
In the present study, we tested whether the size of bulliform phytoliths is an effective statistical indicator for distinguishing between wild rice and domesticated rice, based on the comparative analysis of morphometric data from different rice species. In addition, we attempted to examine how growing conditions, especially climate and water levels, influence bulliform size.
Finally, we discussed how rice bulliform phytolith morphometry can be used in archaeological research.

MATERIALS AND METHODS
In the present study, a total of 24 specimens of Oryza were sampled. Field collection of rice plants was assisted by a team of investigators from Wuhan Botanical Garden and Nanjing Institute of Geology and Palaeontology, CAS, with the permission of the owner or regulatory body for each location. The criteria for categorizing and naming the species collected were accepted from the classification scheme of the genus Oryza in the Flora of China (Liu and Phillips, 2006) (available at http://flora.huh.harvard.edu/china/PDF/PDF22/Oryza.pdf). The specimens of O. sativa and O. rufipogon in Wuhan were both mature when we collected them between September 21 and 23, 2011, but the ripening rate of O. rufipogon was very low due to low temperatures. In Hainan Province, the specimens of O. rufipogon were in anthesis and immature while the specimens of O. officinalis and O. meyeriana were mature and had begun shattering when we collected them between December 1 and 7, 2012. The specimens of O. rufipogon in Chaling, Hunan Province, were immature when we collected them between September 18 and 28, 2010. The specimens of O. rufipogon in Jiangxi and Yunnan Province were already mature when we collected them during October and November 2010, respectively. All plant samples were preserved at the Institute of Geology and Geophysics, CAS, Beijing.
There are different hydrological environments among our sampling sites. In Wuhan, specimens of O. sativa and O. rufipogon were simultaneously cultivated at paddy fields with shallow water (Figures 2a, b), which sometimes needs draining to achieve moderate draught during the pustulation and fruiting   Figures 2h, k); O. meyeriana in Sanya City grows in an understory bush on a hill and prefers shady and dry habitats (Figures 2i, l).
For each rice specimen, we selected all leaf blades from the bottom to the top of a single plant, making sure to keep the leaf blade intact for phytolith extraction. This is because there is significant difference in bulliform phytolith size among different leaf blades of the same plant and different parts of the same leaf blade, with the smaller bulliform phytoliths from the lower leaves . Bulliform size tends to decrease from leaf base to leaf apex of the same leaf blade (Wang et al., 1997). Therefore, variation in bulliform phytoliths from a few randomly selected rice leaves does not reflect the overall data, and only the selection of all intact leaves can guarantee the representativeness and reliability of the data.
Leaf blades were cleaned with distilled water in an ultrasonic cleaner, and then prepared for wet oxidation: 1) all samples were cut into 1-3-cm pieces and placed in 20 ml of 65% saturated nitric acid for over 12 h, then heated in a water bath for 20 min to oxidize organic materials completely. 2) The solutions were centrifuged at 3,000 rpm for 6 min, decanted and rinsed three times with distilled water, and then rinsed with 95% ethanol until the supernatants were clear. 3) The extracted phytoliths were mounted onto microscopic slides in neutral resin, and the residual samples were transferred to storage vials.
A Leica DM750 light microscope at 600× magnification was used for photomicrography and phytolith counting. One hundred or more bulliform phytoliths, including asymmetric types, were counted in each sample, except for the sample of O. meyeriana which produced only a few bulliform phytoliths. Two morphometric parameters, vertical and horizonal lengths (VL and HL), were measured to describe the size of bulliform phytoliths. The measurements were taken from images using the ImageJ software (version 1.48r.). Descriptive statistics of morphometric data were performed using Excel software, and the mean ± SD of each sample was plotted using the Grapher software to perform a comparative analysis. Discriminant function analysis in SPSS 24.0 software was then used to statistically determine the differences in bulliform sizes in different species.
In order to understand whether bulliform phytolith size correlated with environmental conditions and to what degree, a Pearson correlation analysis was performed. Eight environmental variables were chosen: altitude, mean annual precipitation (MAP), July precipitation (MP7), January precipitation (MP1), mean annual temperature (MAT), July temperature (MT7), January temperature (MT1), and relative humidity (HHH). Modern climatic data for the 11 sampling sites (Figure 1) were obtained from the nearest meteorological station to each site, since the spatial variation of climatic parameters exhibits a clear gradient across these locations. These data can be collected from the databases (1981-2010) of the National Meteorological Information Center, China (http://data.cma.cn/). Origin 8.5.1 software was used to conduct the correlation analysis of bulliform morphometrics and environmental variables, which were plotted into scatter plots. A linear regression was inserted into these scatter plots, and then the Pearson correlation coefficients (r) and significances (P) were taken to statistically evaluate the correlation.

Morphological Contrast of Bulliform Phytoliths in the Four Oryza Species
Overall, the bulliform phytoliths in O. sativa and O. rufipogon have similar shapes with an intact circular part of the fan, round arc of the scalloped edge, and ridge-like tubercle on the lateral side (Figure 3). Significant intraspecific morphological variation, however, was also found in both O. sativa (Figures 3a-h) and O. rufipogon (Figures 3i-p). Moreover, we noted that the bulliform phytoliths from the Xishuangbanna O. rufipogon had a very specific shape with a small circular part of the fan, angular scalloped edge, large deep decorations, and without the ridgelike tubercle on the lateral side (Figures 3q-t). This shape was not only different from that of O. sativa but also distinct among other O. rufipogon specimens.
Bulliform phytoliths of O. officinalis are mostly large and full in shape with a rounded arc of the scalloped edge, irregular large deep decorations, longer handles, and a shorter intact circular part of the fan but without the ridge-like tubercle on the lateral side ( Figure S1). From an overall perspective, this shape is similar to that of O. sativa and O. rufipogon.
Bulliform phytoliths of O. meyeriana are generally long and thin and not full in shape like the other varieties. These phytoliths are very small and appear similar to a teardrop or nail, with an angular scalloped edge, small irregular, but deep decorations, a longer handle, and shorter circular part of the fan not intact, without the ridge-like tubercle on the lateral side ( Figure S2). This shape is different from that of the other three Oryza species.

Morphometric Analysis of Bulliform Phytoliths in the Four Oryza Species
Overall, our morphometric data demonstrated that the bulliforms of domesticated rice were not always larger than wild ones, and there was no significant difference in size. Table 1 shows the mean values of vertical length (VL) and horizontal length (HL) of bulliform phytoliths from the studied samples. For original measured data, see the Datasheet S1.
In the only specimen of O. officinalis, the mean values of VL and HL of bulliform phytoliths were 51.13 ± 8.77 and 40.43 ± 7.34 μm, respectively. In all 100 bulliform phytoliths, the maximum VL was 80.00 μm, while the minimum was 32.94 μm; the maximum HL was 62.76 μm, while the minimum was 25.17 μm.
In the only specimen of O. meyeriana, the mean values of VL and HL of bulliform phytoliths were 38.02 ± 5.35 and 29.37 ± 6.29 μm, respectively. In all 34 bulliform phytoliths, the maximum VL was 50.59 μm, while the minimum was 30.49 μm; the maximum HL was 44.26 μm, while the minimum was 18.18 μm. Figure 4 presents a comparison of the sizes of bulliform phytoliths from the studied rice species. The results revealed that the values of VL and HL from O. rufipogon were scattered and widely overlapped with the other three Oryza species. Although the bulliform sizes of O. sativa were larger Parameters VL and HL were important and used in the discriminant function analysis ( Table S1). The following two canonical discriminate functions were used in the analysis: function 1 explained 89.5% of the variance, and function 2 explained 10.5% of the variance. Parameter VL had the largest absolute correlation with function 1, indicating that it contributed most to function 1; parameters HL had the largest absolute correlation with function 2, indicating that it contributed more to function 2 (Table S1). These two functions were used to plot the data (Figure 5) The 24 specimens of rice were then divided into two groups: mature and immature (Figure 6). As can be seen, the sizes of bulliform phytoliths from mature rice were scattered and partly overlapped with the immature rice. In contrast, bulliform phytoliths of immature rice were slightly larger than most mature species. There was no significant difference in bulliform size between mature and immature rice species.
We further compared the bulliform size of rice species in terms of their growing region (Figure 7). In addition, the effect of habitat wetness on bulliform phytolith size was investigated. We found that the bulliform phytolith size of specimens of O. rufipogon native to the warmer and wetter sites of Wenchang, Hainan (field no. Z23), and Dongxiang, Jiangxi (field no. Z1), and cultivated in the paddy fields in Wuhan, was  smaller than that of native species (Figure 8). Similarly, in the Wuhan paddy field, the species of O. rufipogon native to the warmer and wetter sites of Guangdong (field nos. Z13, Z15), Guangxi (field nos. Z12, Z14, Z27), and India (field no. Z41) also have smaller bulliforms (Figure 8). In the tropical Hainan, under parallel climate conditions, O. officinalis in the aquatic environment has a significantly larger bulliform size than those of O. meyeriana in a dry habitat (Figure 8). Therefore, the sizes of bulliform phytoliths from wild rice with preferable water habitats were mostly larger than those of wild rice under relatively dry conditions. Finally, Figure 9 shows the comparison of bulliform phytolith sizes between O. sativa and O. rufipogon growing in the adjacent test paddy field in Wuhan. For O. sativa, the ranges of mean VL and HL of bulliforms fell into 37-45 and 32-48 μm, respectively. For O. rufipogon, the mean VL and HL ranges were 29-41 and 24-34 μm, respectively. Thus, our data also indicated that the bulliform phytoliths of O. rufipogon may be generally smaller than those of O. sativa if they were artificially grown in the same environment.

Correlation Analysis of Bulliform Morphometrics and Environmental Variables
Summary statistics for the eight environmental variables of different sites are given in Table 3. The VL and HL values of bulliform phytoliths plotted against the different environmental variables and the results are shown in Figure 10 and Figure 11. MAP, MAT, MT1, and HHH had a positive correlation (r = 0.438-0.610) with VL and HL, and the linear regression analysis revealed that these correlations were significant (p < 0.01 or 0.05) (Figure 10), demonstrating that the bulliform phytolith sizes are affected by changes in these climatic parameters. In contrast, other environmental characteristics, including MP7, MT7, MP1, and altitude, were not significantly correlated with changes in bulliform sizes (Figure 11). FIGURE 6 | Contrast of Oryza bulliform phytolith sizes from the studied specimens depending on whether they were mature or immature. The sample codes in the figure refer to the field numbers: specimens WN-7, TS-1, and HL-3 were O. rufipogon from Wanning and Wenchang, Hainan Province, and specimen P4 was O. rufipogon from Chaling, Hunan Province. August 2019 | Volume 10 | Article 1037 Frontiers in Plant Science | www.frontiersin.org

Causes of Variations in Rice Bulliform Phytolith Morphometry
The results of the present study indicate that morphometric measurements of bulliform phytoliths from wild and domesticated rice widely overlap (Figure 4), exhibiting little diagnostic potential for taxonomic identification at the species level. These results thus support the conclusion that the morphometry of bulliform phytoliths is not as informative for distinguishing between domesticated rice and wild rice, as suggested by previous studies (e.g., Pearsall et al., 1995;Wang and Lu, 2012;Gu et al., 2013). However, the above morphometric data have often been overlooked, and an increasing number of studies have recently used bulliform phytolith size as a proxy to track the rice domestication process Luo et al., 2016;Zuo et al., 2017;Qiu et al., 2019). The bulliform phytoliths from domesticated rice really were larger than those from wild ones in the same test paddy field in Wuhan city (Figure 9), possibly indicating that domestication may result in an increase in bulliform phytolith size. The genetic and phylogenetic signal for this bulliform size variation has not been well revealed to date. Furthermore, this increase did not necessarily result from domestication and may be caused by other factors.
Previous studies have suggested that two factors, plant maturity and environmental conditions, may affect rice bulliform phytolith size. It is suggested that the bulliform phytoliths from mature rice leaves are usually larger than those from immature leaves (Zheng et al., 2003a;Qin et al., 2006;. This hypothesis may be real for the plants growing in the same location, but is not supported when comparing bulliform size of rice species from different sites (Figure 6). Because environmental conditions are seemingly much more important, and the effect of degree of maturity could be ignored when multiple environmental factors are considered. According to the results of the present study (Figures 10 and 11), from a statistical point of view, we can conclude that the larger rice bulliform phytolith sizes, as defined by the higher VL and HL values, likely occurred at the locations with higher temperature and precipitation. Therefore, the increasing trend in rice bulliform phytolith size in some archaeological records (Zheng et al., 2003b;Zheng et al., 2004;Luo et al., 2016;Zuo et al., 2017;Qiu et al., 2019) may also be caused by climatic changes during the early and middle Holocene when temperature and precipitation were gradually rising.
The present study also revealed that the growing microenvironment, such as water environment, can also influence the size of rice bulliform phytoliths. Rice growing under wetter conditions usually produced larger bulliform phytoliths than those growing under drier conditions (Figures 7 and 8). Therefore, except climate regimes, changes in wet/dry habitat for rice should be considered for the use of bulliform size to track the process of rice domestication.
It should be pointed out that the present study just revealed the hydrothermal condition as one of the environmental factors influencing rice bulliform size. Some factors such as plant genotypes, soil fertility, light period length on photosynthesis, and other abiotic factors may also cause these variations, which were  not controlled or excluded for this study. Conditional plantation experiment under controlling environments and genotypes in test paddy field is needed to further test if and how water levels and temperature can affect the bulliform phytolith size.

Implications for Archaeology of Rice Domestication
Changes in bulliform phytolith size of rice are regulated not only by domestication, which possibly represents genetic changes, but also by environmental factors. Given that environmental factors influence bulliform phytolith size of rice and that the role of genetic background has not yet been firmly established similar to established domestication traits such as non-shattering and increased seed size, bulliform measurement was considered as a semi-domestication trait ). Therefore, the use of rice bulliform phytolith size as an index for determining domesticated plants from their wild ancestors should be conditional. In other words, if the increasing size trend of bulliform phytoliths is used to reveal the process of rice domestication, the influence of hydrothermal conditions should be excluded first.
For further archaeological use of this index, we suggest that: 1) if the time series of rice bulliform phytolith size from a region is long, then the climatic changes (fluctuations in temperature and precipitation) through time should be considered, and the results from quantitative reconstructions of paleoclimate could be used as an independent variable to explain bulliform size variation; 2) the spatial scale of studied regions should be small and without a clear climatic gradient; 3) parallel comparison FIGURE 9 | Contrast of bulliform phytolith sizes from domesticated and wild rice growing at the adjacent test paddy field in Wuhan. Demonstrating that the bulliform phytoliths of these domesticated rice specimens were generally larger compared with wild specimens (Oryza rufipogon) in the same environment. WH, Wuhan. of rice domestication processes in different regions using bulliform size should consider climatic differences between the regions; and 4) the changes in rice arable systems (wet/dry growing conditions) in any studied archaeological sites should be first revealed using the promising sensitive/fixed phytolith morphotype model defined by Weisskopf et al. (2015). Even though the influence of environmental factors has been controlled or excluded, and rice bulliform phytoliths shifting toward larger sizes are interpreted as reflecting the domestication process, it is still not possible to provide a determinate range of bulliform size for identifying domesticated rice, due to the wide overlap observed in the bulliform morphometric data between modern wild and domesticated rice (Figure 4). Thus, rice bulliform phytolith size is a supporting rather than conclusive proxy for determining the domesticated status of rice in archaeological research. Combination of bulliform phytolith size with other established criteria can provide precise identification of wild and domesticated rice.
Finally, notably, frequent gene exchange occurs between domesticated and wild rice, and there is a co-evolutionary relationship between them (Song et al., 2006;Zhao et al., 2010;Ge and Sang, 2011;Choi et al., 2017). Recent large-scale genomic analysis showed that O. rufipogon populations are widely affected by the gene flow of domesticated rice, and that the existing O. rufipogon are actually a hybrid swarm (Wang et al., 2017). This indicates that it is difficult to rule out the interference of the domesticated rice gene flow when using the existing O. rufipogon species as a reference for phytolith morphological analysis. Wild rice populations in the region with higher rice farming intensity are more affected by the introgression of the domesticated rice gene, and the genetic relationship with domesticated rice is closer (Song et al., 2003;Song et al., 2006), leading to the possibility of bias in the morphometric measurements of their bulliform phytoliths. Further research utilizing archaeological rice remains combined with ancient DNA analysis (e.g., Tanaka et al., 2010;Castillo et al., 2016), may reduce the interference of domesticated rice gene flow, and thus generate credible results to establish suitable criteria for distinguishing between domesticated rice and wild rice.

DATA AVAILABILITY
The datasets generated during the current study are available from the corresponding author on reasonable request.

AUTHOR CONTRIBUTIONS
HL and CW designed research. JZ, CW, and LM collected samples. CW and JZ performed experiments. CW and YG analyzed data. CW and HL wrote the paper. All authors read and approved the final manuscript.

ACKNOWLEDGMENTS
We are grateful to Dr. Wen Zhou and Yumei Li for their assistance in the field work on rice sampling; Dr. Yajie Dong and Deke Xu for kind help with data processing. We appreciate the constructive comments from the editor and three reviewers. We also would like to thank Editage (www.editage.cn) for English language editing.