Wetland areas of different states of North-East India (Assam, Meghalaya, and Manipur), Western India (Maharashtra, Goa), and Southern India (Kerala) were explored during the years 2017–2020 to survey and collect biological samples of the genus Nymphaea. Sixteen populations of N. micrantha and 27 populations of N. nouchali consisting of 90 and 92 individuals, respectively, were collected from different locations in Assam, Meghalaya, Manipur, Goa, Maharashtra, and Kerala. The species were identified by Prof. Santhosh Nampy and Dr. A. K. Pradeep, University of Calicut, Kerala, and the voucher specimens were deposited at the herbarium of Department of Botany, University of Delhi (DUH14710-DUH14732) and Botanical Survey of India, Eastern Regional Centre, Meghalaya (ASSAM96585-ASSAM96610). The plant material of 3–5 individuals per population was collected, dried in silica, and brought to the laboratory for DNA extraction. To avoid sampling from the same plant and capture more diversity, the samples were collected from the plants growing at a greater distance than 5 m. Individuals collected from a single wetland were assigned to one population. All the accessions studied in the present study are tabulated in Table 1.

To assess the cross-species transferability of the SSR markers, a few other Nymphaea spp. were also collected (Table 1). Three species, N. Pubescens, N. malabárica, and N. rubra, were collected from the germplasm growing in the wetlands of Malabar botanical garden, Kozhikode, Kerala. Two other species viz. N. × khurooi (earlier named as N. alba var. rubra) and N. caerulea were collected from North-eastern Hill University (NEHU), Shillong, and Yamuna Biodiversity Park, Delhi, respectively.

Genomic DNA was extracted from the silica dried leaf tissues brought from the field using CTAB (hexadecyltrimethylammonium bromide) method described by Doyle and Doyle (1987) with a few modifications. The modifications involved an extended incubation of 1 h in CTAB buffer, an extended RNase treatment of 1 h, a phenol-chloroform extraction (phenol:chloroform: isoamyl alcohol; 25:24:1), and final precipitation by adding 1/10 volume of 3 M sodium acetate and an equal volume of chilled ethanol. The quality and quantity of the extracted DNA were estimated by electrophoresis on a 0.8% agarose gel and using a NanoDrop spectrophotometer (NanoDrop, Wilmington, DE, United States). The final concentration of the DNA was adjusted to 50 ng/μl and the DNA was stored at 4°C.

The genome size of N. micrantha was estimated by flow cytometry. About 1 cm² of the young fresh leaf was used for sample preparation using the CyStain^® PI Absolute P kit (Sysmex, Germany) following manufacturer instructions. Solanum lycopersicum L. “Stupicke’ polnı’ rane”’ 1C = 0.98 pg (Doležel et al., 1992) was used as an internal standard reference. In brief, the nuclei of the standard and the sample were isolated, stained, and analyzed simultaneously. The nuclear DNA content was determined using the CyFlow^® Cube 8 flow cytometer (Sysmex, Germany) for each sample. For both, the standard reference and N. micrantha samples were prepared in replicates, and each was run three times with a minimum of 2,500 nuclei count per sample. The fluorescence values were converted to the DNA content using the following formula:

The high-quality genomic DNA of N. micrantha was used for genome sequencing. The sequencing library was prepared using the NEBNext Ultra DNA library prep kit for Illumina (New England BioLabs, United States) following the manufacturer’s instructions. In brief, the genomic DNA was fragmented by sonication, followed by adaptor ligation (TruSeq Index3 adaptor; GATCGGAAGAGCACACGTCTGAACTCCAGTCACTTAGG CATCTCGTATGCCGTCTTCTGCTTG). The adaptor-ligated fragments were then PCR amplified using NEBnext universal (P3; AATGATACGGCGACCACCGAGATCTACACTCTTTCCCTA CACGACGCTCTTCCGATC*T) and Index3 tagged NEBnext indexed primer (P7; CAAGCAGAAGACGGCATACGAGATG CCTAAGTGACTGGAGTTCAGACGTGTGCTCTTCCGATC) and gel purified. The concentration of the prepared library was quantified by Qubit 4.0 Fluorometer (Thermo Fisher Scientific, United States) and sequenced on Illumina HiSeq 2000 platform in paired-end sequencing mode (Illumina, United States) by Macrogen Inc., (South Korea). The FastQ files containing the raw reads were submitted to the sequence read archive (SRA) at National Centre for Biotechnology Information (NCBI) under the accession number PRJNA750447. The raw data was filtered for high-quality reads [Reads with ≥ 70% HQ bases (Q ≥ 20)] using the NGS QC toolkit (Patel and Jain, 2012) at default parameters. High-quality reads obtained were used for de novo assembly using Velvet (v1.2.10), with the k-mer length of 55. The assembled sequences were filtered to all contigs greater than or equal to 1,000 bp to remove fragmentation of the assembly.

The assembled sequences of the N. micrantha genome were mined for perfect microsatellites using MISA (MicroSAtellite identification tool; Thiel et al., 2003). The minimum numbers of repeats used to select the microsatellites were sixteen for mono-, nine for di-, six for tri-, and five for tetra-, penta- and hexanucleotide repeats. For functional annotation, the SSR containing sequences were searched against the NCBI non-redundant (nr) protein sequence database using the BLASTx program with an e-value cutoff 1e^–5 or better. Further, Gene Ontology (GO) terms were assigned by adopting the terms that had been assigned to the top hit sequences of the database and mapping them to the higher-level categories using the Functional Annotation Tool of DAVID Bioinformatics Resources 6.8 (Huang et al., 2009).

Sequences harboring SSRs with a minimum repeat length of 20 bases and optimum flanking regions (≥50 nucleotides on both flanks of SSR) were used to design primers. The web-based program, Primer3Plus¹, was used to design primer pairs from the SSR flanking regions. The parameters for designing primers were—primer length of 18–25 with an optimum of 22 bases; product size range of 100–500 with an optimum of 250 bases; annealing temperature between 5 and 65°C with an optimum value of 60°C and GC content between 40 and 80%, with an optimum of 60%. Other parameters, such as overall self and 3′ self-complementarity, penalty weight, and sequence quality scores, were kept at the default setting. A total of 217 primer pairs were designed and synthesized at Eurofins Genomics India Pvt. Ltd., India (Supplementary Table 1). Additionally, the 43 genic SSR markers were designed from the transcriptomic data (unpublished data, S. Parveen and S. Goel) (Supplementary Table 1).

The 43 genic-SSR, designated as NymTr_1 to NymT_43 and a subset of 100 non-genic SSR markers, designated as Nym_NGS1 to Nym_NGS100 were chosen for experimental validation. These markers were tested for amplification using two N. micrantha DNA samples. The PCR reaction was conducted in a total reaction volume of 15 μl containing 1x PCR buffer, 2 mM of MgCl₂, 0.2 mM of each dNTP, 0.4 mM each of forward and reverse primers, 50 ng of template genomic DNA, and 1.25 U of Taq DNA polymerase (iNtRON Biotechnology, South Korea). The DNA amplification was performed in a SimpliAmp™ thermal cycler (Applied Biosystems, United States) with the following cycling parameters: initial denaturation at 94°C for 3 min followed by 30 cycles at 94°C for 20 s, optimized primer annealing temperature (between 5 and 60°C) for 20 s, DNA extension at 72°C for 30 s, and a final extension at 72°C for 7 min. The amplified products were initially resolved on 2% agarose and visualized under UV-light after ethidium bromide staining.

Primers producing a clear unambiguous band were selected and further analyzed on 6% polyacrylamide gel for all the accessions of N. micrantha and N. nouchali. Amplicon size was determined by comparing it with a 100 bp DNA ladder (GeneDireX, Inc.). Unambiguously amplified alleles for each SSR locus were scored individually and coded as for their integer size in base pairs (bp). To avoid genotyping errors due to PCR-based amplification failures, the software Micro-Checker 2.2.3 (Van Oosterhout et al., 2004), was used to estimate the errors in allele scoring, presence of null alleles, large allele dropouts, or stuttering during PCR amplification. Null allele frequency (NAF) for each locus was calculated based on the formulas of Chakraborty et al. (1994) and Brookfield (1996). The polymorphic microsatellite markers were further assessed for their cross-species transferability in five Nymphaea species viz. N. pubescence, N. malabarica, N. rubra, N. × khurooi, and N. caerulea (see Supplementary Notes).

Various statistical genetic diversity estimates such as an average number of alleles (N_a), overall, an effective number of alleles (N_e), observed (H_o) and expected (H_e) heterozygosity, F_st = Nei’s differentiation index (F_st), and Gene flow (N_m) were calculated for each locus using GenAlEx 6.5 (Smouse and Peakall, 2012). The polymorphism information content (PI) for each of the primer pairs has been calculated using the package POLYSAT (Clark and Jasieniuk, 2011), with the following formula.

Where n is the number of alleles, and P_i and P_j are the frequencies of alleles i and j, respectively (Botstein et al., 1980). Deviations from Hardy–Weinberg equilibrium were estimated for each locus in each population and across all populations based on Fisher’s exact test using default settings in GENEPOP version 4.7.5 (Rousset et al., 2020). The significance of deviations was adjusted using sequential Bonferroni correction (Rice, 1989). Deviation from hardy-weinberg equilibrium (HWE) for each locus in each population was analyzed using a heat plot (Supplementary Figure 1) generated using the PEGAS package in R (Goss et al., 2014). Pairwise linkage disequilibrium (LD) between SSRs was computed based on log-likelihood ratio statistics (G-test) in GENEPOP.

To estimate within population diversity, the genetic diversity parameters were calculated using GenAlEx 6.5 and poppr R package. Observed mean number of alleles (A_o), the mean number of effective alleles (N_e), percentage of polymorphic loci (PPL), number of private alleles (A_p), and mean observed heterozygosity (H_o) were calculated using Genalex and the Nei’s unbiased gene diversity/expected heterozygosity (H_e; Nei, 1978), Shannon–Wiener index of diversity (H; Shannon, 2001), Simpson’s index (lambda; Simpson, 1949), and evenness (E₅) were calculated using poppr R package. Private alleles were defined as those discovered only in the population considered, discarding those that were found only once as they could reflect genotyping errors.

To evaluate genetic structure in our dataset, we used three different methods: the model-based Bayesian method implemented in the STRUCTURE version 2.3.4 (Pritchard et al., 2010) and maximum likelihood (ML) estimation method implemented in the SNAPCLUST (Beugin et al., 2018); and the model-free DAPC (Discriminant Analysis of Principal Components) (Jombart et al., 2010).

The STRUCTURE uses Markov Chain Monte Carlo (MCMC) approach to estimate every individual’s admixture proportions for a predefined K value (Pritchard et al., 2010). We ran an analysis in the STRUCTURE version 2.3.4 for 150,000 MCMC replications after 50,000 burn-in steps. About 10 replicates each were performed for K values ranging from 1 to 10. The optimum number of populations (K) was estimated using a web-based program, the STRUCTURE HARVESTER (Earl, 2012). The STRUCTURE HARVESTER uses the method outlined by Evanno et al. (2005) and searches for a mode in the ΔK distribution.

The DAPC combines the advantages of principal component analysis (PCA) and discriminant analysis (DA) (Jombart et al., 2010). It first transforms the data using PCA and then performs DA on retained principal components. The number of groups or genetic clusters was defined based on K-means clustering of principal components using find.clusters function in the ADEGENET 2.1.1 package (Jombart, 2008). The optimal K-value was identified by running k-means sequentially for K-value ranging from 1 to 10 and comparing different clustering solutions using BIC (Bayesian Information Criterion). The best-fit K was chosen based on the point at which the elbow in the curve of BIC values as a function of K was observed. The cross-validation function, xvaldapc, was used to determine the minimum number of PCs (Principal Components) to be retained for accurate assignment of every individual to different groups. Finally, the individuals were assigned to clusters based on the posterior membership probabilities of every individual as a measure of the admixture proportion originating in each cluster.

SNAPCLUST allies the advantages of both model-based approaches as used in STRUCTURE and geometric approaches like those used by DAPC and provides a fast maximum-likelihood solution to the specific genetic clustering problem (Beugin et al., 2018). We used Snapclust.choose.k to choose the more adequate number of clusters (K). The most suitable K-value was selected through Bayesian Information Criterion (BIC) by choosing the point at which the BIC value was lowest. After choosing K, we proceeded with the clustering analysis using the function snapclust, implemented in R package ADEGENET.

To further assess the population genetic structure and evaluate the genetic interrelationships among the individuals collected from different locations, unrooted neighbor joining (NJ), and principal coordinate analysis (PCoA) were performed based on pairwise Nei’s genetic distance matrix. The PCoA was generated using dudi. pco function implemented in ade4 R-package (Dray and Dufour, 2007), and the NJ tree was constructed using NJ function from ape R package (Paradis and Schliep, 2019).

The analysis of molecular variance (AMOVA) analysis was carried out to test the structure of genetic variation within and among genetic clusters inferred based on population structure analysis.

We performed hierarchical AMOVA analysis using the poppr R package (Kamvar et al., 2014) to quantify the partitioning of variance within the individuals, among individuals within the genetic clusters, and among genetic clusters by using pairwise F_ST as the distance measure. Within individual variance was calculated to quantify variation due to heterozygosity within genotype. The F-statistics was calculated to summarize the degree of differentiation at each hierarchical level. The significance of variance components and F-statistics was tested using 9,999 random permutations of data matrices.

To analyze the dynamics of genetic diversity within each genetic cluster, the measures of genetic diversity were estimated for each cluster obtained during population structure analyses. The genetic diversity indices: observed mean number of alleles (A_o), the mean number of effective alleles (N_e), percentage of polymorphic loci (PPL), number of private alleles (A_p), and mean observed heterozygosity (H_o) were estimated using GenAlEx 6.5 (Smouse and Peakall, 2012), and Nei’s unbiased gene diversity/expected heterozygosity (He), Shannon–Wiener index of diversity (H), Simpson’s index (lambda), and evenness (E5) were calculated using poppr R package. Weir and Cockerham’s (1984) estimator of pairwise F_st between accession groups/clusters were estimated using Hierfstat (Goudet, 2005).

To determine if the genetic variability is structured in geographic space, we performed an isolation by distance (IBD) test (Mantel, 1967) between geographic populations and within each genetic cluster using procedures implemented in the GenAlEx version 6.5 (Peakall and Smouse, 2006). The Mantel test was performed using the Slatkin’s (Slatkin, 1995) pairwise linearized F_ST {F_ST [F_ST/(1 − F_ST)]} and natural log-transformed geographic distance [ln (1 + geographic distance)] matrix as outlined by Rousset (1997). The statistical significance of the test was conducted by random shuffling (10,000 times) of all individuals among the geographic locations.

The spatial autocorrelation analysis was computed using the multi-locus genetic distance between individuals calculated via the method of Peakall et al. (1995) and was further explained in Smouse and Peakall (1999), and the geographic distance calculated as Euclidean distance using the coordinates of the samples’ collection sites (Peakall and Smouse, 2006). The analysis was performed with an even distance class of 6 km and 10,000 permutations. For testing the statistical significance, a two-tailed 95% confidence interval was constructed around the null hypothesis of no genetic structure (r = 0). A bootstrap of 10,000 trials was used to estimate a 95% confidence interval around the observed r values for each class.

Besides making a record of pheno-events covering onset/duration of flowering, fruiting, and production of vegetative propagules through regular field visits (2017–2020), the plantlets/propagules were grown in experimental ponds at the Department of Botany, the University of Delhi to monitor the events. The pollinated flowers/fruits (n = 20) of both species were collected from the natural populations and dissected to observe the seed set.

As N. micrantha exhibited a complete absence of genetic diversity and compromised reproductive output in nature, some of the key aspects of floral biology influencing seed/fruit-set were investigated in the plants growing in experimental ponds at the Department of Botany.

The genome size estimation of N. micrantha, based on the linear relationship of 2C peaks of the sample and the internal standard, indicated an approximate genome size of 1.11 pg per 2C content corresponding to 544 Mbp/1C (Supplementary Figure 2). Illumina paired-end sequencing resulted in 56,645,378 raw reads providing a genome coverage of ∼10 X fold based on the estimated genome size. The raw reads were subjected to quality check to remove sequence artifacts, such as low-quality reads and adapter contamination, using the NGS QC toolkit (Patel and Jain, 2012). A total of 53,511,392 (94.57%) high-quality reads showing a Q value of ≥ 20 were obtained after quality filtration. An assembly of high-quality filtered reads, with Velvet (v1.2.10) at a k-mer value of 55, resulted in 518,464 contigs with an average contig length of 596 bp (Supplementary Table 3). To reduce fragmentation of the assembly, the contigs were further filtered to retain contigs ≥ 1,000 bp resulting in 77,175 contigs with an average contig length of 1,536 bp (Supplementary Figure 3). These contigs were used for the identification of microsatellites. A total of 22,268 non-redundant microsatellite loci were obtained using the Perl script MISA. There were 16,937 sequences containing SSRs of which 4,123 sequences contained more than one SSR. About 2,464 SSRs were present in the compound formation, which were removed from further analysis (see Supplementary Notes and Supplementary Table 4 for details).

To identify the functional significance of 12,556 non-redundant SSR containing sequences, a BLASTx analysis was performed against the NCBI non-redundant (nr) protein database. A total of 3,809 (30.3%) sequences were found to have similarities to the sequences in the NCBI nr protein database, and 1,623 (12.9%) were found to have at least one functional annotation. The remaining 8,747 (69.67%) did not show significant similarity to the known sequences and were not annotated.

Gene Ontology (GO) enrichment analysis, using DAVID Bioinformatics Resources 6.8 (Huang et al., 2009), classified microsatellite sequences into three categories: cellular components (CC), molecular functions (MFs), and biological process (BP). Supplementary Table 5 provides detailed information regarding annotated SSR sequences. Fourteen hundred and thirty-six contigs were mapped to GO terms with 1,009, 1,172, and 892 assignments distributed under the cellular component, molecular function, and biological process ontology, respectively.

The 143 SSR primer pairs were analyzed for polymorphism among 90 N. micrantha individuals. Surprisingly, all the markers tested in N. micrantha produced the same DNA profile and were monomorphic.

The flowering and fruiting in N. nouchali were observed from March to October; while the N. micrantha flowers throughout the year, with peak flowering residing between April and September. Vegetative propagules were observed only in the case of N. micrantha arising from the petiolar node (Figures 4b,c). They appeared throughout the year, with copious production during the winters (November–January). Fruits of N. nouchali matured within 25 DAP (Figure 4g,h) and dehisced underwater, with multiple longitudinal splits in the capsule, producing an enormous number of seeds. On the contrary, fruit formation in the naturally growing plants of N. micrantha was not apparent (Figure 4).

The flowers exhibited protogyny, as the stigma became receptive on the day of anthesis, while anthers dehisced a day later. There was a gradual decline in the receptivity after the first day, indicating overlapping dichogamy in the species. Pollen tubes reached the ovules 6 HAP (Hours after pollination). However, only 4.33 ± 1.91% of the ovules showed receptivity. The fertility of pollen grains was significantly greater (87.89 ± 5.18%) than the viability (18.02 ± 4.51%; t = 44.36, df = 38, p < 0.001). Controlled pollinations (both selfed and outcrossed) resulted in the initiation of fruits formation but failed to mature into fruits. The pollinated flowers turned brown and started decomposing after 12–15 days of the anthesis phase (Figure 4d). Such abortive fruits were dissected to ascertain seed formation, showed a 2.8 ± 0.84 % seed set (Figure 4e). No difference in seed set was observed between selfed and cross-pollinated flowers. The results of seed set in the naturally growing populations were also comparable with manual pollination treatment.

The 57 polymorphic SSR markers produced for N. nouchali were also assessed for their cross-species transferability in five other Nymphaea spp. viz. N. malabarica, N. caerulea, N. rubra, N. pubescens, and N. × khooroi (see Supplementary Table 6). The markers with a high transferability rate would be of importance and can be efficiently utilized in conducting basic genetic studies in Nymphaea.

Nymphaea micrantha is non-native to India and, thus, is expected to possess low genetic variation. Surprisingly, however, no genetic variation was found within or between populations investigated as a single multilocus genotype (MLG) was detected in the samples studied.

Clonal reproduction is often considered an evolutionary dead-end (Stebbins, 1957). Lambertini et al. (2010) detected low genetic diversity in three clonal aquatic species, Elodea canadensis, Egeria densa, and Lagarosiphon major. Ren et al. (2005) reported very low genetic differentiation in Eichhornia crassipes (water hyacinth) populations from different localities in China. Ellstrand and Roose (1987) summarized the population genetic structure of 21 clonal plant species and revealed that two species, Taraxacum obliquum, and Gaura triangulata, did not display any apparent intraspecific variability. Pollux et al. (2007) reported that the shifting balance between sexual and asexual reproduction affected the genotypic diversity within populations Sparganium emersum. These studies suggest that the plants reproducing asexually can preserve well-adapted genomes in stable environments, which permits the fixation of well-adapted phenotypes over generations (Niklas and Cobb, 2017).

Various species of Nymphaea are known to be propagated by both seeds and bulbils (Wiersema, 1988; Ansari and Jeeja, 2009). The prevalent mode of reproduction (asexual vs. sexual) in these plants may likely influence the extent of genetic variability of the species. Our results showed that sexual reproduction in N. micrantha is significantly impaired in all the investigated populations. Manually pollinated flowers were aborted at different stages due to low pollen viability (male fitness), and poor ovule receptivity (female fitness). These results were comparable to the flowers that were marked for natural fruit-set (open-pollination). As seen in the present work as well, the species is known to propagate predominantly through leaf propagules (Wiersema, 1988). The extensive clonality accompanied by impaired sexual reproduction could be the major cause for the observed lack of genetic diversity in N. micrantha. Moreover, N. micrantha is not indigenous to India (Ansari and Jeeja, 2009; POWO, 2021), although, the time and source of the introduction are not known. It is predicted that populations in their introduced range are less variable relative to their source or native range, owing to founder effects and genetic drift (Husband and Barrett, 1991; Amsellem et al., 2000; Hagenblad et al., 2015). The successful expansion of its distribution range in India demonstrates its strong adaptability to the local ecological conditions. The non-existence of genetic variation within and between N. micrantha populations, therefore, may be attributed to (1) predominant vegetative propagation, broad tolerance range, and plasticity of genotypes; (2) introduction of the same genotype once or multiple times in the North-eastern parts of India; (3) fixation of well-adapted genomes over generations.

Our results suggest that a single genotype was introduced to India and expanded its distribution range through massive vegetative reproduction. The ability of N. micrantha to conquer water bodies in the North-East regions of India as a silent invader could also be attributed to human interventions. The ornamental and religious values of the species aided by human involvement might have played a role in its spread. Several previous studies have also demonstrated the presence of low genetic diversity in invasive and clonally reproducing plant species supporting our results. Xu et al. (2003) reported extremely low genetic diversity in the invasive alligator weed in China. Li et al. (2006) revealed a lack of genetic variation in the clonally reproducing invasive plant Eichhornia crassipes.