For SSR marker analysis in this study, 182 O. rufipogon accessions from 10 different countries viz. Bangladesh (n = 8), Cambodia (n = 31), Thailand (n = 24), Myanmar (n = 16), Taiwan (n = 4), Vietnam (n = 20), Nepal (n = 20), Laos (n = 8), Papua New Guinea (n = 13), and India (n = 38) were undertaken. These germplasm accessions were originally procured either from the International Rice Research Institute (IRRI), Philippines or from National Rice Research Institute (NRRI), Cuttack and being actively maintained at Punjab Agricultural University (PAU), Ludhiana. The rice cultivars, Punjab Rice 114 (PR114), Punjab Rice 121 (PR121), Punjab Rice 122 (PR122), and Punjab Basmati 3 (PB3) were used as negative checks. The upland rice variety Vandana was selected as a positive control due to the presence of 90 kb of phosphorus uptake 1 (Pup1) locus. The list of accessions undertaken along with their country of origin is given in Supplementary Table S1. Out of 182 accessions, 67 were sequenced for complete coding sequences (CDS) of PSTOL1 (Table 1).

Genomic DNA of 182 accessions along with cultivated varieties was isolated using modified Cetyltrimethyl ammonium bromide (CTAB) method of Saghai-Maroof et al. (1984). Eleven SSR markers comprised of 7 dominant SSR markers (Pup1-K41, Pup1-K42, Pup1-K43, Pup1-K46, Pup1-K48, Pup1-K52, and Pup1-K59) in INDEL region and 4 co-dominant (Pup1-K4, Pup1-K5, Pup1-K20, and Pup1-K29) markers located in Pup1 genomic region (Chin et al., 2010, 2011) were used for SSRs genotyping (Supplementary Table S2, Supplementary Figure S1). PCR amplification was performed in a 20 ul reaction mix with the following thermal conditions: 94°C for 4 min, followed by 35 cycles of 94°C for 1 min, 55°C for 1 min and 72°C for 1 min and the final extension of 7 min at 72°C.

The Oryza sativa cv. Kasalath sequence (Accession AB458444.1) from position 275,525–276,499 bp covering 975 bp CDS region of PSTOL1 was used for designing sequencing primers (PSTOL1 forward: 5′-ATAGCAGGCATTTCTGGCTCA-3′ and PSTOL1 reverse: 5′-CCATGACAGCTGATTGCCTT-3′). The amplicons were purified using Wizard® SV 96 PCR clean up/Gel extraction kit from Promega, USA, as per manufacturer's protocol. Sequencing reaction was performed using ABI Big-dye Terminator v3.1 chemistry and sequenced using ABI Sequencer 3730XL. Hi-fidelity long-read DNA polymerase (Phusion Taq) from Promega, USA, was employed to obtain the required amplicon size. A minimum of three replications was carried out for the confirmation of single nucleotide polymorphism (SNPs).

For comparative sequence analysis, the PSTOL1 sequences were trimmed to remove any poor quality region at both ends. Multiple sequence alignment was performed using Clustal W of MEGA version 7.0 (Kumar et al., 2016). Kasalath sequence was used as a reference for detection and determination of SNPs position among the PSTOL1 sequences obtained from O. rufipogon accessions. The identified SNPs were manually confirmed using chromatograms. DnaSP version 5.0 and Selecton server (http://selecton.tau.ac.il/, Stern et al., 2007) were used to calculate summary statistics for nucleotide diversity (π), the number of segregating sites, non-synonymous (k_a), and synonymous (k_s) mutations and the ratio of k_a/k_s is to estimate positive/purifying selection of a given amino acid, the number of haplotype and Tajima's D test.

Bioinformatics toolkit (http://toolkit.tuebingen.mpg.de/) was used to predict protein structures of all sequences. Homology modeling approach was employed using the Modeler to determine the structure of proteins based on the known structure of template protein. Protein domains were predicted and compared using Pfam (http://pfam.xfam.org/search) and Prosite (http://prosite.expasy.org) online tools. The protein models were checked for the quality using the Ramachandran plot developed using Procheck through PDBsum (http://www.ebi.ac.uk/thornton-srv/databases/pdbsum). The modeled protein structure was visualized and compared in UCSF Chimera (Pettersen et al., 2004). All the structures were superimposed for observing structural variations.

A phylogenetic tree was generated by MEGA7.0 software using the alignment file obtained earlier. The molecular phylogeny was inferred using the Maximum Likelihood method with 1,000 bootstrap (Tamura and Nei, 1993). All positions containing gaps and missing data were eliminated along with other default settings of the software.

For functional validation of PSTOL1 haplotypes toward phosphorus uptake efficiency, eight accessions with seven different haplotypes along with positive control Vandana and negative control PR121 were grown in replicates under low and high phosphorus conditions in the greenhouse facility following the protocol of Gamuyao et al. (2012). High and low P growth conditions were established by maintaining the NaH₂PO₄concentration in the hydroponic media as 100 μM and 10 μM, respectively. The eight accessions (CR 100013, CR 10013A-H2, IRGC 104639-H3, IRGC 104712-H4, IRGC 100588-H8, IRGC 105569-H9, IRGC 81989-H11, and IRGC 106506-H17) along with controls were grown in hydroponics for about 2 months. Due to poor germination of accessions representing remaining haplotypes, they were not included in the present study. The seeds were germinated on the wet filter paper, and four seedling replicates per accession were assayed for each phosphorus treatment. After 10 days of germination, seedlings were transferred to the Styrofoam trays suspended in Yoshida growth media (Yoshida et al., 1976). The nutrient media was changed at every third day. Data for the root length, shoot length, final dry root, and shoot weight were taken after 60 days in growth media. Phosphorus content in roots was measured using Inductively Coupled Plasma Spectrophotometer after digestion in a mixture of HNO₃, HClO₄, and H₂SO₄ (3:1:1) according to the protocol described by Neelam et al. (2011). The morphological data on root and shoot traits under study along with phosphorus content on dry root weight basis was subjected to the statistical analysis. Student's t-test was applied for testing the significance of differences among the means of O. rufipogon accessions and the controls.

The analyzed co-dominant SSR markers were found monomorphic among all the 182 O. rufipogon accessions and indica rice cultivars (PR114, PR121, PR122, PB3, and Vandana) (Supplementary Table S3). For dominant markers (Pup1-K41 to K-59), the presence of Vandana alleles was detected in the majority of the O. rufipogon accessions as well as in the modern rice cultivars except for marker K-46. Rice cultivars PR114, PR121, and PR122 did not show any amplification for K-46 marker. This indicates the specificity of K-46 marker for the assessment of phosphorus starvation tolerance.

From comparative sequence analysis, 53 nucleotide changes (52 SNPs and 1 nonsense mutation) across the exon were observed (Table 2). Both types of conversions i.e., transitions (n = 39) and transversions (n = 14) were observed, while the G/A transition was most common (28.30%). Higher transitions indicated more of the synonymous substitutions were present among genotypes and hence no conformational changes in the structure of proteins were observed. Based on the nucleotide diversity present among O. rufipogon accessions, haplotypes were identified using DnaSP software v5.0. Out of 53 identified, 10 SNPs at position 174, 260, 283, 303, 358, 410, 554, 633, 647, and 738 were found as singleton whereas 43 SNPs were parsimony informative sites with a minimum frequency of occurrence in two or more O. rufipogon accessions (Figure 1). The overall nucleotide diversity (π) of the identified PSTOL1 alleles was found 0.00758, which indicates low variance in the average number of nucleotide differences per site between two sequences. The number of mutations (n = 53) and the number of segregations sites (S = 53) were same, suggesting their positive selection. The value of Tajima's D obtained is negative (−1.09788) supporting the above-said statement. Presence of fewer haplotypes was observed than the number of segregating sites indicating the lower frequency of rare alleles present in the population. A total of 17 haplotype groups was formed, revealing genotypes divergence at PSTOL1 gene among studied O. rufipogon accessions (Table 3). The haplotype H1 carried two O. rufipogon accessions, one each from Vietnam and Thailand with only one segregating site at position 816. Major haplotype group H2, harbors 44 O. rufipogon accessions along with Vandana and Kasalath indicating the sequence similarity among them. The other 16 haplotype groups had O. rufipogon accessions ranging from one to three. The haplotype H3 and H4 shares same phylogenetic clade, but having different nucleotide segregating sites. The O. rufipogon accessions (IRGC 106150, IRGC 106156, IRGC 93200, IRGC 100588, IRGC 83819, and IRGC 93215) under the smaller phylogenetic clade B grouped into different haplotype i.e., H7, H8, H13, and H16 indicating the presence of rich allelic divergence for the PSTOL1 gene in these accessions.

A total of 28 differences in amino acid sequences with a comparison to the variety Vandana and Kasalath were identified (Table 2). The ratio of non-synonymous/synonymous site (k_a/k_s) was found 1.52, suggesting that the amino acids were under positive selection and favored by the environment. The amino acids at position 25, 31, 35, 55, 72, 85, 87, 95, 101, 115, 120, 127, 142, 146, 156, 157, 161, 185, 202, 209, 216, 219, 246, 253, 273, and 283 highlighted by yellow color were under positive selection (Supplementary Figure S2).

Protein structures for Kasalath and 67 accessions of O. rufipogon belonging to different haplotype groups were superimposed and analyzed for structural differences. The non-synonymous mutations concentrate around the ATP binding site (LEU45, ARG47, GLY48, VAL53, ALA65, GLU112, MET114, TYR113, SER118, LYS168, GLN170, and LEU173). (Supplementary Figure S3). The Ramachandran plot of Kasalath revealed more than 99.3% residues were at the core and allowed region and only two residues were present in the disallowed region. Similar results were obtained for protein models of other accessions. All the accessions had a three-dimensional structure similar to the reference Kasalath except accession IRGC 106336 (Figure 2). The protein structure of Kasalath and other accessions displayed 14 helices and two beta-pleated sheets and nine strands, while accession IR 106336 showed only three helices, one sheet, and five strands. In accession IRGC 106336, the PSTOL1 sequence revealed the presence of premature stop codon at position 137 and its further domain analysis using Pfam and Prosite revealed that it encodes partial protein kinase domain instead of protein tyrosine kinase as encoded by Kasalath. The Prosite analysis for PSTOL1 protein in Kasalath predicted the features as the kinase domain from codon 39–319, nucleotide phosphate binding site (NP_BIND) at position 45–53, ATP-binding site (BINDING) at position 67 and the proton acceptor site as active site at position 166 whereas IRGC 106336 had partial protein kinase domain from 39 to 136, NP_BIND from position 45–53, ATP-binding site (BINDING) at position 67, with absence of proton acceptor site i.e., active site (Supplementary Table S4).

Phylogenetic analysis at PSTOL1 locus revealed divergence among O. rufipogon accessions (Figure 3). The different node colors correspond to the different mutations present in O. rufipogon accessions and vice-versa. Two major groups were observed. The clade A is consisted of 59 O. rufipogon accessions which can be further divided into 8 subgroups. Out of 59 accessions, 44 were found to have similar sequences as that of the reference sequence, whereas others having either single or more substitutions as compared to the reference. Among 44 accessions 6 are of Nepal origin, 6 are of Cambodia origin, 9 are of Indian origin, 8 are of Thailand origin, 4 are of Papua New Guinea origin, 4 are of Vietnam origin, 1 is of Laos origin, and 6 are of Myanmar origin. This indicates a common evolutionary relationship of these O. rufipogon accessions with aus type variety Kasalath. The clade B includes eight accessions, two from Laos (IRGC 106150, IRGC 106156), three from Nepal (IRGC 93200, IRGC 93215, IRGC 93283) and one each from Vietnam (IRGC 83819), Taiwan (IRGC 100588), Myanmar (IRGC 81989).

Genotypic variation for root and shoot length, dry root and shoot weight and phosphorus content on dry roots basis were examined under phosphorus sufficient and deficient conditions after 2 months of the experiment (Tables 4A,B). Genotypic differences were observed among all O. rufipogon accessions under both growing conditions. All of the genotypes under phosphorus sufficient conditions had almost double root length as compared to the deficient conditions. Though, not much difference was observed in shoot length and shoot weight for all the genotypes when compared under phosphorus sufficiency and deficiency conditions. The O. rufipogon accession IRGC 106506 (H17) showed the best root and shoot length under phosphorus-deficient conditions when compared to other genotypes and control. Among H2 haplotype O. rufipogon accession CR 100013A performed better than Vandana for all the traits studied. In terms of root and shoot weight, O. rufipogon IGCC 106506 (H17) showed the best root weight followed by IRGC 81989 from H11 haplotype. Approximately, 1.5 and 2.3 times higher phosphorus content was found in the O. rufipogon accession IRGC 106506 when compared to indica cultivar PR 121 and Vandana respectively. Other than that, O. rufipogon accession IRGC 104639 from H3 haplotype also showed 1.2 times higher phosphorus content when compared to PR121 indicating their potentiality toward improving elite cultivars for phosphorus uptake.

Our results with Pup1 specific markers on O. rufipogon accessions and O. sativa cultivars revealed no allelic differences for almost all markers other than K-46 and K-05 (Supplementary Table S3). Chin et al. (2010), also observed the Kasalath specific alleles for markers K-41, K-43, and K-48 in lowland/irrigated (indica, japonica, aus, and traditional or modern) rice cultivars, representing nonusefulness of these markers for marker aided selection for phosphorus uptake efficiency. Similarly, the markers K-42, and K-29 were not found linked with PUE by Sarkar et al. (2011) while assessing indica germplasm. It should be taken into consideration that Gamuyao et al. (2012) ruled out other co-dominant and INDEL markers as indicative of PUE except for K-46. This dominant marker was found useful for MAS in the progenies involving Kasalath as Pup1 donor variety and Asian lowland rice varieties (without this gene) by Chin et al. (2010, 2011) and Pariasca-Tanaka et al. (2014), supporting our results. Mukherjee et al. (2014), assessed 108 genotypes from different states of India for phosphorus acquisition efficiency with gene specific markers and closely linked SSR marker RM1261 and reported no association between markers and PUE. The same has been observed when they studied a RIL population developed from a cross between Gobindabhog (with PSTOL1 gene) and Satabdi (PSTOL1 absent). The notion that PSTOL1 specific marker was not indicative in the case of indica germplasm (Mukherjee et al., 2014) is more likely due to the complex nature of Pup1 locus and different genetic background and environment where this gene has to express.

The germplasm survey with Pup1 specific markers of Kasalath indicated entire inserted region of 90 kb among studied O. rufipogon accessions. The probable explanation for this could be a continuous gene flow between O. sativa and O. rufipogon populations throughout the history of domestication (Vaughan et al., 2008). Also, O. rufipogon accessions from South and Southeast Asia are considered as the wild progenitor of domesticated rice (Oka, 1988; Molina et al., 2011) and hence chances of recent hybridization events needs to be accounted for the phenomena. While studying allelic diversity at PSTOL1, Vigueira et al. (2016) reported presence/absence polymorphism in 12 of the O. rufipogon accessions out of 40 studied, along with the loss of function mutation in one accession and 56 synonymous and nonsynonymous substitutions in 28. He explains this phenomenon as long- term balancing selection at PSTOL1 locus for maintaining both functional and non-functional alleles among the accessions of O. rufipogon and indica, aus, tropical japonica cultivars. Though, none of the alleles found conferring superior phenotype than Kasalath in their study, whereas in our case functional alleles were observed.

Our results on molecular diversity at PSTOL1 locus, suggests the presence of lower diversity among O. rufipogon accessions from South Asia and Southeast Asian nations. This is expected as they share common geographical boundaries. This result is in accordance with several studies conducted on the assessment of genetic diversity of Asian wild rice using RFLP, microsatellite markers, SINEs, sequence based polymorphism, ISSRs, chloroplast, and low-copy nuclear markers (Joshi et al., 2000; Cheng et al., 2002; Rakshit et al., 2007; Xu et al., 2007; Huang et al., 2012). In a study, conducted by Huang et al. (2012) on the phylogeography of Asian wild rice using 42 genome-wide sequence tagged sites demonstrated that O. rufipogon accessions were grouped into two genetically distinct clades (Ruf-I and Ruf-II). The O. rufipogon accessions from South Asia and Indochinese Peninsula (Thailand, Myanmar, Cambodia, Vietnam, and Laos), were clustered into one group (Ruf-II), supporting our results. The presence of few accessions from Nepal, Laos, Vietnam, Cambodia, and Myanmar into intermediate and second major clade may be interpreted as an admixture. Moreover, the humid tropical plain areas in the Indo Peninsula zone act as a transitional region for evolutionary studies is likely the reasons for observed admixture.

It is always worthwhile to look for better alleles of a gene for creating and maintaining natural genetic diversity. Our results demonstrated the presence of 17 different haplotypes within 975 bp sequence of PSTOL1 locus indicating rich nucleotide variation among studied O. rufipogon accessions. Two of the O. rufipogon accessions under H17 and H11 were found performing better than the positive control under phosphorus-deficient conditions. Functional allelic variants were observed and utilized for improving various agronomically important traits in cereal crops (Ellis and Setter, 1999; Bhullar et al., 2010; Ravensdale et al., 2012; Vasudevan et al., 2014; Ashkani et al., 2015). The wheat powdery mildew resistance gene Pm3 with 17 identified functional alleles is a remarkable example of natural variations present in GenBank accessions and can be efficiently utilized for conferring broad-spectrum disease resistance. For rice blast resistance gene, functional orthologous have been found in wild rice O. rufipogon accessions (Lv et al., 2013; Xu et al., 2014; Ashkani et al., 2015) defining their utility in widening the genetic base of cultivated rice varieties. A novel allele of PSTOL1 gene is identified in O. glaberrima (CG14) and being transferred to the NERICAs (New Rice for Africa) cultivars using allele-specific markers. In their study, they identified 3 novel alleles in 10 studied O. rufipogon accessions and also the presence of Kasalath alleles for INDEL markers which is consistent with our results. The successful efforts for the transfer of PSTOL1 were made by Gamuyao et al. (2012) through marker assisted backcross breeding to Asian rice cultivar IR74 with increased root growth and phosphorus uptake efficiency.

The presence of PSTOL1 in all O. rufipogon accessions raise the question regarding the functionality of different alleles under phosphorus deficiency. Validation of haplotype groups showed the significant difference for root and shoot length and biomass as compared to PR121 and Vandana under both phosphorus sufficient and deficient conditions. The correlation between root elongation, higher root and shoot biomass of genotypes under P-deficiency is considered as one of an important indicator of higher phosphorus uptake efficiency. A number of reports, including QTLs on P deficiency induced root elongation in plants were published (Steingrobe et al., 2001; He et al., 2003; Ma et al., 2003; Wissuwa, 2005; Li et al., 2007; Rose et al., 2013). Near isogenic line of “Nipponbare” with Pup1 QTL from “Kasalath” showed high P content, high tillering and high root growth under P-deficient upland conditions (Wissuwa and Ae, 2001; Wissuwa et al., 2002). The O. rufipogon accession IRGC 106506 showed the highest root growth under P-deficiency and thus is the best option for transferring this novel allele to elite cultivars for improving P starvation tolerance.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.