Genetic Differentiation Revealed by Selective Loci of Drought-Responding EST-SSRs between Upland and Lowland Rice in China

Upland and lowland rice (Oryza sativa L.) represent two of the most important rice ecotypes adapted to ago-ecosystems with contrasting soil-water conditions. Upland rice, domesticated in the water-limited environment, contains valuable drought-resistant characters that can be used in water-saving breeding. Knowledge about the divergence between upland and lowland rice will provide valuable cues for the evolution of drought-resistance in rice. Genetic differentiation between upland and lowland rice was explored by 47 Simple Sequence Repeats (SSRs) located in drought responding expressed sequence tags (ESTs) among 377 rice landraces. The morphological traits of drought-resistance were evaluated in the field experiments. Different outlier loci were detected in the japonica and indica subspecies, respectively. Considerable genetic differentiation between upland and lowland rice on these outlier loci was estimated in japonica (Fst = 0.258) and indica (Fst = 0.127). Furthermore, populations of the upland and lowland ecotypes were clustered separately on these outlier loci. A significant correlation between genetic distance matrices and the dissimilarity matrices of drought-resistant traits was determined, indicating a certain relationship between the upland-lowland rice differentiation and the drought-resistance. Divergent selections occur between upland and lowland rice on the drought-resistance as the Qsts of some drought-resistant traits are significantly higher than the neutral Fst. In addition, the upland- and lowland-preferable alleles responded differently among ecotypes or allelic types under osmotic stress. This shows the evolutionary signature of drought resistance at the gene expression level. The findings of this study can strengthen our understanding of the evolution of drought-resistance in rice with significant implications in the improvement of rice drought-resistance.


Introduction
How plants have developed resistance to biotic and abiotic stresses is a fundamental question in plant biology. To address this question, plant ecotypes grown in contrasting environments are studied as they are likely under divergent selections, resulting in adaptive divergence [1]. Although ecotypes of wild species adapted to different ecosystems are served as ideal systems for studying the evolution of resistance to biotic and abiotic stresses [2][3][4][5], literatures on different ecotypes of a crop to a specific stress are still rare.
The Asian cultivated rice (Oryza sativa L.) is one of the most important cereal crops in the world by providing stable food for . 50% of the total global population (ref). However, rice production has encountered severe challenges due to frequent droughts and water shortages [6]. The utilization of natural variations of watersaving and drought-resistant characters in upland rice is an effective solution to improve the drought-resistance for rice [7][8].
Upland rice has adapted to the water-limited and rain-fed rice ecosystems, facing the higher risk of drought [9]. On the contrary, lowland rice is commonly planted in fields with irrigation facilities with relatively lower risks of drought. Upland rice may accumulate more drought-resistant genetic variances than lowland rice, leading to potential adaptive divergence of this rice ecotype on drought-resistance [7]. Therefore, studying differentiation between upland and lowland rice provides ideal systems for understanding the evolution of drought resistance and exploring drought-resistant genetic resources in crops. Previous attempts using neutral markers to explore differentiation between upland and lowland rice were not very successful [10][11], while another study using the functional-based markers has detected some ecotype-specific genetic features between upland and lowland rice [12].
The simple sequence repeats (SSRs) located in the ESTs (expressed sequence tags) are developed and widely applied in population genetics and breeding since the last decade [13]. The expression of drought-responding EST-SSRs is up-or down-regulated by drought stress. Although these genic SSRs are not always selective among ecotypes, they have a higher probability of being associated with the drought-resistance. In this study, 47 drought-responding EST-SSR loci were selected to explore the genetic differentiation between upland and lowland rice ecotypes from China. The questions to be addressed are: 1) Does considerable differentiation occur on droughtresponding EST-SSR loci between upland and lowland rice ecotypes? 2) Is the differentiation between upland and lowland rice related to drought-resistance? 3) How do the upland-and lowlanddifferential alleles of the selective loci respond to drought stress? Answering these questions will facilitate our understanding of the evolution of drought-resistance in rice and give cues in rice breeding for drought-resistance.

Plant materials
A total of 377 rice landraces from Yunnan, Guizhou, Guangxi, Jiangsu, and Hebei province in China were used to study the differentiation between upland and lowland rice ecotypes (Table 1). Based on the information provided by the National Core Germplasm Database (http://crop.agridata.cn/A010110. asp), the landraces of upland rice from the five provinces were accounted for ,80% of the total upland rice germplasm in China. The predefined subspecies and ecotypes of rice landraces were provided by the institutes that collected them. Thus, crop materials in this experiment can be divided into four groups: jap-upland, jap-lowland, ind-upland, and ind-lowland. Additionally, 22 common wild rice strains were also used as reference while studying the differentiation among these two rice ecotypes. Twenty-two drought-resistant traits of 56 japonica and 49 indica materials were further evaluated in the field experiments rice landraces.

Rice DNA extraction and SSR genotyping
Rice genomic DNA was extracted from the 10-day-old seedlings following the common cetyltrimethyl ammonium bromide (CTAB) protocol. Three individual seedlings of a landrace were mixed together to include the genetic variations within a material. To focus on the drought-related genetic features, 47 polymorphic drought-responding EST-SSRs (SSRs locating in drought-responding ESTs) across 12 rice chromosomes were selected for genotyping (Table S1). The three primers system was applied in the polymerase chain reaction (PCR), containing a pair of EST-SSR primers and a rice-universal M13 sequence. The M13 sequence (59 to 39: CACGACGTTGTAAAACGAC) was labeled with fluorescent dyes at the 59 end. The forward primer of the EST-SSR was also joint with the rice common M13 sequence at the 59 end. The 20 ul reaction system contains 16buffer (Mg2+), 2 mM each of dNTP, 50 ng of genomic DNA, 1 U of Taq  polymerase, 16 mM of fluorescent dye labeled M13 sequence,  4 mM of the SSR primer liked with M13, and 20 mM of the other  SSR primer. The PCR program started with a denaturation  period of 4 min at 94uC, followed by 34 cycles of 40 s at 94uC,  30 s at their respective Tm (Table S1), and 40 s at 72uC, and ended after 10 min final extension at 72uC. The PCR products were then analyzed on ABI 3130XL (Applied Biosystems, USA) using ROX500 as the internal standard. The resulting chromatograms were analyzed and scored by Peakscanner ver. 1.0.

Field drought-resistance evaluation
The field evaluation of the drought-resistant traits was conducted at the Zhuanghang experimental station of Shanghai Academy of Agriculture Science in 2012. Two treatments were designed as control (CK) and drought (DT). Plants in CK were grown in a paddy field following the normal irrigated rice cultivation, while plants in DT were grown in a similar field nearby. The surface level of the DT field was approximately 2 meters higher than the CK field to simulate the upland condition and drain away the water more quickly. The DT field had a drought-resistance screen facility, which equipped a moveable platform capable of stretching out to keep the rainfall away when necessary. The germinated seeds were directly seeded both in the CK and DT on May 30 th . However, water in the DT field was drained away on July 15 th and it was never irrigated until August 25 th . Individuals of each rice landrace were planted in a plot with 364 grids of 20 cm intervals between two seedlings. Consequently, there were 105 plots in one experimental block. Three replicates were included and arranged following the randomized complete block design. The number of stomas per leaf area, the excised leaf water (EWR) loss rate in 4 hours, and the root-shoot ratio were measured from leaves of main tillers in CK 30 days after the beginning of drought treatment. Relative water content (RWC) of leaf samples was measured 16 days after drought both in CK and DT. The content of malonaldehyde (MDA) was measured from fresh leaf tissues 40 days after drought following the protocol of the test kit (Nanjing Jiancheng Bioengineer Institute). The flag leaf length and width were measured at the heading stage both in CK and DT. All morphological traits were measured from three individuals from each replicate. The yield related traits were measured from four individuals of each replicate after harvest. The general drought resistance of a rice landrace was quantified by the droughtresistant index (DRI) [14]. It was calculated as P d /P C *(P d /P ad ), where P d is the yield production in drought treatment, P C is the yield production in control treatment, and P ad is the average yield production of total materials in drought treatment.

Missing data and null alleles in data scoring
For genotyping, any null alleles were genotyped twice to avoid the manipulation errors. 337 putative null alleles were detected out of the total 18753 individual-loci combinations (1.8%), which was much lower than that expected by chance at the 5% level. The putative null alleles were treated as missing data in the further analyses.

Population structure and genetic diversity analysis
The model-based program STRUCTURE ver. 2.3.3 was used to infer population structure using the admixture model with a burn-in length of 50,000 and a run length of 50,000 for Markov Chain Monte Carlo iterations. Ten simulations were run for each K from K = 1 to K = 8 and the results were generated together using software CLUMPP_Windows.1.1.2. The Evanno's K was applied to determine the inferred K value [15]. The Principal Coordinate Analysis (PCoA) implemented in GenAlex ver. 6.43 was conducted to investigate the genetic distances between rice landraces and wild rice using the total 47 SSR loci. Genetic diversity was quantified by four estimators as gene diversity, number of alleles, heterozygosity, and polymorphism information

Neutrality tests for EST-SSR loci
Most of the genic markers used in this study were neutral among ecotypes as very low Fsts were recorded among rice ecotypes, similar with that using putatively neutral genomic SSRs [11]. Three Fst-based outlier tests were included to detect signs of selection among these genic SSRs respectively in the indica and japonica subspecies as indicated by the population structure. First, LOSITAN was used to detect loci under selection at the 95% confidence level [16]. An initial run with 100,000 simulations was conducted and followed by computing the distribution of neutral Fst using the putatively neutral loci derived from the initial run with 500,000 simulations. Second, the hierarchical-Bayesian method developed by Foll and Gaggiotti (2008) was applied to detect outlier loci using the software BAYESCAN (http://cmpg. unibe.ch/software/bayescan/) on the 95% posterior probabilities [17]. The parameters were set to 10 pilot runs of 5,000 iterations and additional burn-in of 50,000 iterations. Third, the software DETSEL 1.0 was applied to do pairwise comparisons to identify the outliers. For DETSEL, the coalescent simulation was conducted with the following parameters: mutation rate (infinite allele model) 0.005, 0.001, and 0.0001; ancestral population size Ne = 500, 1,000, and 10,000; population size before the split N 0 = 100; time since an assumed bottleneck T 0 = 50, 100, and 1,000 generations; and time since divergence t = 100 generations. Loci falling outside the specified ''probability region (95% levels)'' were considered as outliers [18]. Based on these outlier tests, two levels of loci sets could be defined from the total loci set as: (1) the neutral loci, which were not detected in any of the three tests, and (2) the decisive selective loci, which were detected at least twice by the three outlier tests.

Genetic differentiation among upland and lowland rice
The general level of genetic differentiation among different groups was quantified by pairewise Fst calculated in Arlequin ver. 3.1 using total loci set with 1,000 permutations. To estimate the level of differentiation on selective loci, Fst was calculated from the decisive selective loci set among upland and lowland ecotypes in japonica and indica, respectively. Japonica landraces were from 5 provinces and indica landrace were from 4 provinces. Thus, there were 10 japonica and 8 indica populations included in this study. To explore the influence of the genetic differentiation on population structures, genetic distance (GD) matrices were calculated by GenAlex ver. 6.43 with 999 permutations using the corresponding neutral or decisive selective loci. The genetic distance matrix obtained form the calculation was used for the cluster analysis according to the un-weighted pairgroup method with arithmetic averages (UPGMA) via the software NTSYS ver. 2.10e.
Testing whether the upland-lowland rice differentiation was related to drought-resistance According to previous studies, the comparison of quantitative genetic divergence (Qst) and neutral genetic divergence (Fst) can be used to detect adaptive evolution. If the Qst is significantly higher than the Fst calculated by the neutral loci, it indicates that the directional selection drives phenotypic divergence and results in ecological adaptation [19]. In this study, the Qst of each trait was calculated as: Qst = V pop /(V pop +2V ind ), in which V pop was the variance among-population and V ind was the variance withinpopulation. Although the SSRs located in the ESTs, the Fst generated among upland and lowland rice was as low as genomic SSRs [11]. This suggests most of these markers were neutral among the two ecotypes. The 95% confidence interval of Qst was calculated using R program with 1000 bootstraps. The Fst derived from neutral loci and its standard errors were calculated by Arlequin ver. 3.1 with 1000 permutations. Any significant differences between the Qst and Fst at p = 0.05 level was determined when |Qst2Fst|.2SQRT (SE Qst 2 +SE Fst 2 ). To test whether the genetic differentiation between upland and lowland rice was resulting from divergent selections on droughtresistance, Mantel tests were conducted between the individual based GD matrix using the neutral or decisive selective loci and matrix constructed from differentiated morphological traits. The differentiated traits were defined as the traits having significant differences between upland and lowland rice by independent t-test using SPSS ver. 15.0. The individual based GD matrix was constructed by GenAlex ver. 6.43. The dissimilarity matrices of differentiated traits were constructed by Hierarchical Cluster Analysis on SPSS ver. 15.0.

Test the expression of ecotype-specific alleles of the selective loci
The decisive selective loci may play a role during upland rice adapting to water-limited upland conditions. To test this hypothesis, the expression levels of upland-and lowland-preferable alleles of the decisive selective loci (E647 and E1899 in japonica, E647 and E1177 in indica) were measured. Here, upland-or lowland-preferable allele was defined as a predominant allele in upland or lowland rice whose frequency was 20% higher than that in lowland or upland ecotype. The allele frequency was calculated as the number of counts of an allele/the total number of counts of all alleles in a given group via the software GenAlex ver. 6.43. Four landraces of each ecotype conferring the ecotype-preferable alleles were included. Three-week-old seedlings in the growth chamber were treated with 20% polyethylene glycol (PEG) 6000 to cause osmotic stress (DT). Some materials were kept in nutrition solutions as controls (CK). When most of the seedlings showed signs of leaf rolling 5 hours after treatment, leaves of three seedlings in DT and CK were collected and kept in liquid nitrogen until RNA extraction. The RNA extraction was followed by the kits of TRNzol A + (TianGen Biotech Co. Ltd.). The level of gene expressions were quantified by RealTime-PCR using SYBR R premix Ex Taq (Takara Bio. Inc.). Actin was used for reference. The expression change value was calculated as: the expression level in DT over that in CK. In addition, to test whether the ecotype-preferable alleles of the selective loci had any impacts on rice drought-resistance, the drought-resistant traits were further compared among upland and lowland rice containing differently expressed ecotype-preferable alleles (E647 and E1899 in japonica, and E1177 in indica).

Population structure and general genetic diversity
A total of 477 alleles were detected in the 47 drought responding EST-SSR loci with number of alleles per locus ranging from 4-25. When running the STRUCTURE simulation using the total materials, there was a sharp peak in Evanno's DK at  Figure S1a). In japonica, the sharp peak occurs when K = 2 (DK = 7.6), which divided landraces into two groups as upland and lowland rice occupied the majority in each ( Figure  S1b). In indica, the sharp peak was also occurs when K = 2 (DK = 2760.1), which separated rice landraces of Guangxi from other regions ( Figure S1c). However, there was a lower peak when K = 4 (DK = 78.0), in which upland rice (white) was distinguishable from lowland rice (thin gray) in Guangxi ( Figure S1d). Similar to the results of STRUCTURE, PCoA indicated that japonica and indica landraces were separated along the first coordinate (x-axis) while common wild rice located at the center (Figure 1), suggesting that upland-lowland rice differentiation should be separately analyzed in japonica and indica subspecies. Differentiation between upland and lowland ecotypes was more apparent in japonica than that in indica (Figure 1). The level of genetic diversity was similar between upland and lowland rice ecotypes, while subspecies (japonica vs. indica) and species (O. sativa vs. O. rufipogon) showed considerable differences ( Table 2).

Selective loci detected by neutrality tests
Among the 47 drought-induced EST-SSR loci, 5, 9, and 5 loci were detected to be selective by Lositan, BayeScan, and Detsel, respectively, in japonica subspecies. Five loci were detected to be the decisive selective loci (Table 3). In the indica subspecies, 2, 2 and 7 loci were detected to be selective by Lositan, BayeScan, and Detsel, respectively. Three loci were considered as the decisive selective loci in indica (Table 3).

Genetic differentiation among upland and lowland rice
The values of Fst were as low as 0.0714 in japonica or 0.0302 in indica types of upland and lowland rice. The value was much lower than that between japonica and indica varieties (Table 4). This result was similar to a previous study in which Fst among different ecotype was recorded as 0.068 in japonica using putative neutral SSR markers [11]. However, the values of Fst between upland and lowland rice became much higher in japonica (0.285) and indica (0.127) when using their respective decisive selective loci. In the UPGMA clusters, upland and lowland rice populations were grouped separately when using the decisive selective loci, while populations of landraces from the same regions were clustered together when the neutral loci were used (Figure 2). These results suggest differentiation occurred between upland and lowland rice generally on the selective loci.

Divergent selection on drought-resistant traits
Drought resistant traits were evaluated in the field experiment. Significant differences were detected on 11 traits between upland and lowland ecotypes in japonica, while significant differences were detected on only 2 traits in indica ( Table 5). The comparison between Qst with neutral Fst was used to test any divergent selection on the drought-resistant traits. The values of the Qst ranged from 0.009,0.171 in japonica and from 0.009,0.072 in Table 5. 22 drought-resistant traits measured in control treatment (CK) and drought treatment (DT) and their Qst calculated in Jap-upland/lowland and Ind-upland/lowland pairs (mean 6 standard error). indica. The Qsts of the most differentiated traits were significantly higher than the neutral Fsts in japonica (0.089860.0024) and indica (0.036360.0017) materials (Table 5). This result suggested that some of the drought-resistant traits were under divergent selection during domestication.
The genetic differentiation associated with morphological differentiation on drought-resistant traits between upland and lowland rice Given that some drought-resistant traits were differentially selected in upland and lowland rice, their genetic differentiation may be associated with the evolution of drought-resistance. To test this, Mantel test was conducted between the individual-based genetic distance (DG) matrix and dissimilarity matrix constructed from morphological traits. In japonica, the dissimilarity matrix constructed from the differentiated traits was significantly correlated with the GD matrix using decisive selective loci (Figure 3b), while it was not correlated with the neutral GD matrix (Figure 3a). These results suggested a strong association between the differentiations on the drought-resistance and on the selective loci in japonica. However, the morphological dissimilarity matrix was neither significantly correlated with GD matrix from the selective loci nor with that from the neutral loci in indica (Figure 3c, 3d).

Expression of ecotype-preferable alleles and their impacts on drought-resistant traits
In japonica subspecies, the allele-6 of locus E647 was uplandpreferable (taking account for 87.9% in upland rice but only 50.6% in lowland rice), while allele-2 was lowland-preferable (taking account for 40.4% in lowland rice but 0 in upland rice).
These two alleles expressed similarly in normal conditions (CK). However, the expression of allele-6 down-regulated largely (. 50%) while the expression of allele-2 remained at the same level in DT. Thus, the expression change values of the allele-2 (1.14660.232) were marginally higher than the allele-6 (0.62560.158, p,0.10) (Figure 4a). The allele-4 of E1899 was upland-preferable (taking account for 71.0% in upland rice but only 16.7% in lowland rice) and the allele-3 was lowlandpreferable (taking account for 79.5% in lowland rice but 23.8% in upland rice). This EST was up-regulated in lowland rice ecotype while its expression was down-regulated in upland rice under DT. Thus, the expression change value of the allele-3 in lowland rice (1.83260.583) was marginally higher than the allele-4 in upland rice (0.41360.111, p,0.10) (Figure 4b).
In indica subspecies, the allele-6 of E647 was upland-preferable (taking account for 71.4% in upland rice but 50.7% in lowland rice), while the allele-7 was lowland-preferable (taking account for 49.3% in lowland rice but only 23.7% in upland rice). The allele-8 of E1177 was upland-preferable (taking account for 44.1% in upland rice but only 17.3% in lowland), while the allele-9 lowlandpreferable (accounted for 74.7% in lowland rice but only 44.5% in upland). The expression change values of E647 were similar among ecotypes or allelic types (Figure 4c), while that of E1177 was much higher in lowland ecotype or lowland-preferable alleles than in upland rice or upland-preferable ( Figure 4d).
As these ESTs differently expressed between allelic types or ecotypes, the drought-resistant traits were then compared among rice ecotypes conferring their ecotype-preferable alleles. As expected, many significant differences not previously detected among total upland and lowland rice were now detected (Table  S2). For example, the upland rice conferring the allele-8 and lowland rice conferring allele-9 of E1177 exhibited significant difference on the root-shoot ratio in indica. This was accordant to the annotated function of gene Os06g0633300 (Table 6). Other genes containing the selective EST-SSRs were also considered to be related with resistance to abiotic stress given their annotations ( Table 6), suggesting that these selective loci might play roles in rice drought-resistance.

Considerable level of upland-lowland rice differentiation on selective loci
As in previously studies, obvious differentiation was reported on drought-resistant morphological traits between upland and lowland rice [9]. However, the general level of genetic differentiation between upland and lowland rice was very low in our study,  similar to the results in previous studies [10][11]. However, Fst calculated based on the outlier loci was considerably high between upland and lowland rice ecotypes, matching the level of Fst between O. sativa and O. rufipogon (0.252) in our study. These results suggested that adaptive divergence among upland and lowland rice occurred on these selective loci. Strong selections on the adaptive loci always affected population structures [20][21]. For example, wheat grown in different ecosystems was clustered separately based on the drought-responding gene TaSnRK2.7 [21]. Similar results were found in this study as revealed by the UPGMA cluster, providing evidences that these selective EST-SSRs received uniform selections among regions. It is noteworthy that morphological differentiation between upland and lowland ecotypes was also greater in japonica than in indica, consistent with the genetic data in these two subspecies.
Upland-lowland rice differentiation was driven by divergent selection on drought-resistance Theoretically, drought-resistance alleles are more strongly selected in upland rice than that in lowland rice, due to its adaption to the water-limited environment. In this study, the general drought-resistance gained from field experiments was much higher in upland rice than that in lowland rice. The higher Qst than neutral Fst suggests some drought-resistant traits are under the divergent selections, leading genetic differentiation on the outlier loci between upland and lowland rice. The genetic differentiation between upland and lowland rice is likely associated with drought-resistance based on the Mantel tests between the genetic distance and morphological dissimilarity, especially in japonica. The formation of upland-or lowland-preferable alleles always result from divergent selections [22]. Different expression of the ecotype-preferable alleles under their favorable conditions could be considered as the signature of natural selection [23]. In this study, we found different expression patterns among uplandand lowland-preferable alleles encountering osmotic stress, adding further evidence of the adaptive divergence between upland and lowland rice. Based on these results, we conclude that genetic differentiation between upland and lowland rice is likely driven by divergent selection on the drought-resistance.

Candidate genes for drought-resistance gene from the selective loci
Recently studies disclosed that it was possible to indentify potential functional genes via studying adaptive divergence in natural populations [24], crops [5], and even in rice [25]. In this study, we found potential candidate genes that may play roles in rice drought-resistance. For example, Os06g0633300 (E1177) encodes the rice phytosulfokine 1 precursor (OsPSK1). Its overexpression can promote rice cell division [26] and root development [27]. Interestingly, upland and lowland rice conferring different ecotype-referable alleles of E1177 exhibited significant differences on root-shoot ratio, a key character of droughtavoidance. Besides, genes with similar annotations as Os01g0607400 (E647) and Os12g0563600 (E1899) were also reported to be associated with plant stress responses by previous studies in other plant species [28][29]. These results provide strong evidences that these candidate genes should have played some roles in rice drought-resistance. However, their functions in rice need further investigation. In a word, we can explore more drought-resistant genes by studying the genetic divergence among the upland and lowland rice.

Conclusions
Crops adapted to different agro-ecosystems always promote the variation of agricultural important genes. In this study, several outlier loci were detected between upland and low rice ecotypes. A considerable degree of genetic differentiation between upland and lowland rice was detected both at the DNA and gene expression level on these outlier loci. Results from this study reveals that the genetic differentiation among the two rice ecotypes is most likely driven by divergent selection on drought-resistant traits. The findings of this study not only help us to understand the underlying molecular basis of adaptive divergence, but also provide valuable implications for rice domestication and breeding, especially on the drought-resistance in rice.   Table S2 22 drought-resistant traits measured in upland and lowland rice materials containing their ecotype-preferable alleles (indicated in parentheses) in control treatment (CK) and drought treatment (DT) (mean ± standard error). The values in bold and with ''*'' indicated significant differences between upland and lowland ecotypes. (DOC)