A total of 102 accessions of J. curcas were collected from the following 11 countries: 72 from China, 7 from India, 7 from Vietnam, 5 from Burma, 4 from Thailand, 2 from Indonesia, 1 from Laos, 1 from Zambia, 1 from Burkina Faso, 1 from Mali, and 1 from Brazil. The accessions from China are highly representative of within-country dispersal, with 69 accessions collected from 6 provinces in southwestern and southern China, including 3 hybrid varieties. These materials were kindly provided by the Tropical Eco-Agriculture Institute of the Yunnan Academy of Agricultural Sciences and Institute of Tropical Bioscience and Biotechnology, CATAS. The ITS spacer region of 102 J. curcas accessions were sequenced in this study, and 20 additional sequences (from plants in India, Madagascar, Cape Verde, Spain, Brazil, and Mexico) were obtained from GenBank (National Center for Biotechnology Information; NCBI). The sample code, location, latitude, longitude, altitude, country, haplotype, and GenBank accession number of all accessions are listed in Table 1. Jatropha gossypifolia was used as an outgroup for phylogenetic analysis [1].

Total genomic DNA was extracted with a Plant Genomic DNA Kit (TIANGEN Biotech, Beijing, China). The nrITS spacer regions were amplified using external primers of the following sequences: ITS1 (TCCGTAGGTGAACCTGCGG) and ITS4 (TCCTCCGCTTATTGATATGC) [16]. PCR was conducted in a 50-μL mixture reaction volume containing 5.0 μL 10 × Taq Buffer, 1 μL dNTP Mix (10 mM each), 0.5 μL Taq DNA polymerase (5 U/μL), 3.0 μL 25 mM MgCl₂, 2.0 μL of each primer (10 μM), 2.0 μL of genomic template DNA (2.5 ng/μL), and additional ddH₂O to the final volume (Vazyme Biotech, Nanijing, China). The PCR amplification procedure for nrITS was 3 min at 94 °C for predenaturation followed by 35 cycles of 1 min at 94 °C for denaturation, 1 min annealing at 56 °C, and 1 min at 72 °C for primer extension; this was followed by a final primer extension of 10 min at 72 °C on a BIO-RAD S1000™ Thermal cycler. Purified samples were sequenced by BGI China. All sequences generated in this study were deposited in GenBank under accession numbers: KP190940 _ KP191040 and KP191042 (Table 1).

Multiple sequences were aligned using the Clustal W algorithm in MEGA 6.06 [17] and refined by manual adjustment with BioEdit version 7.2.5 [18]. The homogeneity of the base composition was detected for Id-test, nucleotide substitutions, and transition/transversion ratio and was calculated with MEGA version 6.06 [17]. Substitution saturation was estimated and transversions were calculated under the F84 model using the software DAMBE version 5.3.8 [19]. DnaSP 5.10 [20] was used to estimate several molecular diversity indices including haplotype diversity (Hd), nucleotide diversity (π), Watterson's diversity of segregating sites (θ_w), and the average of nucleotide differences (k). Estimates of evolutionary divergence were calculated using the Maximum Composite Likelihood Model. Selection neutrality was tested to detect historical demographic expansions by Tajima's D, Fu and Li's D, and Fu's Fs methods. The demographic history was assessed by the distribution of pairwise sequence differences (mismatch distribution) and site frequency spectra using the program DnaSP 5.10 [21].

Bayesian inference (BI) analysis was performed using MrBayes version 3.1.2 [22]. The GTR + G model was identified as the best-fit model using MrModel Test 2.3 [23]. A four-chain Markov Chain Monte Carlo (MCMC) was run for 1,000,000 generations and two simultaneous analyses were performed. Trees were sampled every 100 generations. The first 5000 trees were discarded as burn-in, and the remaining trees were used to construct a 50% majority rule consensus phylogram. Maximum parsimony (MP) analysis was implemented using PAUP version 4.0b10[24]. The following were in effect: heuristic searches with 1000 random addition sequence replicates, ten trees were saved at each step, tree-bisection-reconnection (TBR) was used to swap branches, and MULTrees option. A MP network topology between haplotypes was reconstructed using Network 4.6.1.3, followed by the median joining (MJ) algorithm to resolve intermediate nodes [25].

A total of 122 nrITS sequences from J. curcas accessions and one from J. gossypifolia were aligned with a length of 614 bp. Sequence characteristics and nucleotide polymorphisms were estimated including the total number of sites (excluding sites with gaps/missing information; n), conserved sites, variable sites, parsimony-informative sites, segregating sites, identical pairs (ii), transitional pairs (si), transversional pairs (sv), and transition/transversion ratio (R). All these characteristics are listed in Table 2. The average nucleotide frequencies were 20.29% A, 35.33% C, 14.67% T, and 29.71% G. The G + C content varied from 64.3% in MX1 to 65.1% in CHH12 and CHS1. The average G + C content (65.04%) was considerably higher than the average A + T content (34.96%). Maximum likelihood estimation of substitution patterns and matrix rates were estimated under the Tamura–Nei model (Table 3). Transitional substitution rates of A/G, T/C, C/T and G/A were 5.93%, 18.28%, 7.59%, and 4.05%, respectively. A test for substitution saturation revealed a nearly linear regression, indicating that mutation of different sequences showed no saturation effects (Fig. 1).

All parameters linked to the genetic variation of ITS were estimated among the accessions. Genetic parameters were estimated by k, π, Hd, and θ_w, which exhibited values of 1.5944, 0.0027, 0.5405, and 0.0160, respectively (Table 4). The nucleotide datasets revealed moderate genetic diversity of nrITS among accessions. A total of 20 haplotypes were identified among 122 J. curcas individuals. The pairwise genetic divergence (GD) between ITS haplotypes was evaluated and ranged from 0.000 to 0.017, whereas the average of all haplotype sequences was 0.002. The GD value between H17 (MX2), H4 (CHGZ3), and H5 (CHGZ4) exhibited the highest divergence of 0.017.

Neutrality tests were performed to determine whether the ITS region was subject to selection or evolved neutrally. Both values were investigated with Tajima's D and Fu and Li's tests. The value for Tajima's D was –2.62985 (P < 0.001), while that for Fu and Li’ s D was –8.25200 (P < 0.02) and that for Fu and Li's F was –7.07538 (P < 0.02; Table 4). The significant negative values suggest that the variation deviates from selective neutrality. To clarify the cause of deviation from neutrality, Fu's Fs was also estimated (Fu's Fs = –4.768; P < 0.005). Tajima's D, Fu and Li's, and Fu's tests consistently revealed significant negative values for J. curcas lineages, which supports selection and non-neutral evolution. This suggests a bottleneck effect and recent demographic expansion of intergenic spacers of ribosomal DNA in J. curcas accessions [26,27]. The mismatch distributions consisted of a significant unimodal curve (r = 0.0815; P = 0.5535) consistent with the hypothesis of sudden population expansion related to the distribution of J. curcas (Fig. 2).

A 50% majority rule consensus tree was inferred from Bayesian analysis with posterior probabilities (PP). The parsimony analysis revealed 60 most parsimonious trees, with a consistency index (CI) of 0.9333 and a retention index (RI) of 0.9259. The MP topology (not shown) was congruent with the BI tree, which was well resolved (Fig. 3). The BI tree revealed a moderately resolved phylogeny and the J. curcas accessions were distinctly divided into three phylogeographic groups: Groups A, B, and C. Accessions of J. curcas in Group A were noticeably split into three major lineages: Clade I, Clade II, and Clade III, which were sister to Group B. Clade I included eight accessions: India (2), Madagascar (1), Cape Verde (2), and Mexico (3). Clade II consisted of 12 accessions: Brazil (1), Madagascar (1), Spain (1), Cape Verde (2), and Mexico (7). Clade III comprised 17 accessions: China (13), India (1), and Vietnam (3). Group B had a short branch and included 84 accessions: China (58), Burma (5), Indonesia (2), India (6), Laos (1), Thailand (4), Vietnam (4), Brazil (1), Burkina Faso (1), Mali (1), and Zambia (1). Group C consisted of an independent clade comprising a single Chinese accession (CHS1).

The haplotype network revealed a significant genealogical relationship among J. curcas accessions, which was generally compatible with the phylogenetic analysis (Fig. 4). H11 was in the central position of the network and was composed of seven individuals from Spain (1), Cape Verde (2), and Mexico (4), connecting H2, H3, H13, H14, H16 and one missing haplotype. These additional haplotypes were one mutational step removed from H11 and exhibited a typical star-like genealogy structure, which suggests that H11 represent an ancestral haplotype [28]. Another single dominant haplotype (H1) was present in 82 individuals and distributed in a central position. This appeared to be a founder haplotype for H3, H6, H7, H8, and H13. This star-like pattern of H1 in the network suggests a rapid demographic expansion and evolution from a small founding population [28,29]. In addition, haplotypes H17, H18, H19, and H20 from Mexico were identified and exhibited numerous mutational steps from H11, suggesting remote phylogenetic relationship and significant genetic divergence among Mexican haplotypes.

Genetic diversity, estimated with k, Hd, π, and θ_w revealed a moderate amount of genetic variation and diversity. Variation in nrITS among J. curcas accessions revealed higher haplotype diversity (Hd = 0.5405) and lower nucleotide diversity (π = 0.0027), which may be associated with particular ecological conditions, a low level of recombination, and few base variations in the nuclear genome [30]. It is likely that the species encountered genetic bottlenecks and founder effects during its geographic spread [31,32]. The highest level of genetic divergence was evaluated as 0.017 between haplotypes H17 (MX2), H4 (CHGZ3), and H5 (CHGZ4), and the genetic distances between Mexican haplotypes were evidently higher than those of other regions (Table 5). There was significant genetic divergence between Mexican and Asian, African, and South American haplotypes, suggesting restricted gene flow between adjacent populations caused by geographical isolation, abundant vegetative propagation, and asexual reproduction [11–13].

The combined results of the phylogenetic and network analyses were moderately resolved, and the accessions of J. curcas were consistently divided into three genetically distinct population lineages (Groups A, B, and C). Group A was separated into Clades I, II, and III and included all the primary Mexican accessions. Group B was composed of the major lineage accessions from Asia, especially from southern and southwestern China. Combined with the data from the network analysis, our results suggest that the American and Asian accessions of J. curcas were identified as two distinct nuclear ribosomal lineages. Phylogenetic relationships between American and Asian lineages demonstrated remote geographical isolation and restricted gene flow, generating significant genetic divergence between American accessions and accessions from other regions [11,13,33].

BI phylogenetic trees showed that the primary American lineages (Mexico and Brazil) were clustered together with accessions from Cape Verde, Spain, Madagascar, and India, forming a distinct sister subgroup with the Asian lineage (China, India, and Vietnam). This suggests a relatively close phylogenetic relationship and evolutionary history and indicates that the Asian lineage may have originated from the American lineage [11,13,14]. The phylogenetic trees showed that accessions from Asia, Africa, and South America clustered together with extremely short internal branches in Group B, and shared a single dominant haplotype (H1). This indicates high genetic similarity a common origin, and a slow rate of genetic divergence among population lineages, suggesting a rapid radiation and recent demographic expansion caused by a founder effect [31,32]. The pattern of haplotype diversity, the significant negative values for the neutrality tests, and the mismatch distribution were consistent with this conclusion. Both the phylogenetic trees and the network analysis identified the presence of two separate lineages among the Chinese accessions, indicating that China's complicated geography and diverse ecosystems led to restricted dispersal and gene flow, and promoted genetic diversification in southwestern China [34].

J. curcas is believed to be native to Central America or Mexico, where it occurs naturally in forests in coastal regions [35]. It was then perhaps widely disseminated from its center of origin by the Portuguese, who transported it to the Old World [3]. The Ceara State in Brazil was also identified as a possible center of origin; however, this proposal was unfounded [36]. Despite this, the center of origin of J. curcas has remained controversial for many years. Our data confirms that Mexico or Central America is the center of biodiversity and the possible center of origin, due to the dominant haplotype frequency and the high genetic variability of accessions from this area [7,13].

Regarding the migration route, texts have mentioned the use of J. curcas prior to 1810 [3]. Previous research assumed that the Portuguese distributed J. curcas from the Caribbean via Cape Verde to Africa and Asia [3,14]. However, the migration and dispersal route remained uncertain. Our BI tree showed initial clustering of CV3 and MX1 (in Clade II of Group A), followed by the formation of a sister subclade, which suggests that the Cape Verdean accessions had prior phylogenetic relatedness with Mexican accessions, and then with the Spanish accessions. The cluster of CV2 (or CV1) and MA2 in Clade I indicates the close relationship between Cape Verdean and Madagascan accessions. Meanwhile, the genetic divergences between H1 and H14 and between H1 and H15 were 0.003 and 0.002, respectively, suggesting that accessions from Madagascar are closely related to Asian accessions, especially to the Chinese accessions. Furthermore, one Brazilian and three African accessions clustered with Asian accessions in Group B and shared one common haplotype in H1, suggesting their close relatedness and migratory connectivity. Combining the above data analyses, we deduced that J. curcas was distributed from Mexico through the Caribbean and/or Brazil by the Portuguese, via Cape Verde where there was a presumed bifurcation point, with one route going to Spain and then to China, further spreading to other regions of Southeast Asia. Another route was dispersed to Africa, via Madagascar and migrated to China, and then further spread to other regions. Our research tentatively supports previous research and provides further details of the migration route (Fig. 5).