Assessment of the elite accessions of bael [Aegle marmelos (L.) Corr.] in Sri Lanka based on morphometric, organoleptic, and elemental properties of the fruits and phylogenetic relationships

Aegle marmelos L. (Bael) is a native tree fruit species in the Indian subcontinent and Southeast Asia. Bael is a popular fruit because of its significant nutritional and medicinal properties. However, bael is an underutilized fruit species in Sri Lanka. Thus, Fruit Crop Research and Development Station of the Department of Agriculture of Sri Lanka has selected five elite bael accessions (Beheth Beli, Paragammana, Mawanella, Rambukkana, and Polonnaruwa-Supun). We assessed these five accessions for the variation of the fruit size, pulp, organoleptic preference, elemental properties, genetic diversity, and evolutionary history. The fruits at the golden-ripe stage were collected during the peak fruiting seasons in 2015, 2016, and 2017. The fruit size, pulp, shell thickness, and seed size were measured and subjected to the General Linear Model (GLM) and Principal Component (PC) Analyses. The fruit pulp was distributed among a group of 30 taste-panelists to rank for the parameters: external appearance, flesh color, aroma, texture, sweetness, and overall preference. The rank data were subjected to association and PC analyses. The elemental contents of the fruit pulp samples were measured using Inductively coupled plasma mass spectrometry and subjected to GLM and PC analyses. We observed a significant diversity in fruit size, organoleptic preference, and elemental contents among bael accessions. Rambukkana and Polonnaruwa-Supun yield the biggest and most preferred fruits. We used trnH-psbA, atpB-rbcL spacer, matk-trnT spacer, and trnL markers to construct phylogenies. Sri Lankan bael split from an Indian counterpart, approximately 8.52 MYA in the Pliocene epoch. However, broader germplasm of Indian bael must be assessed to see the presence of any independent evolution within Sri Lanka.

Introduction Aegle marmelos L. (Bael) is an indigenous tree fruit species in the Indian subcontinent and Southeast Asia. Bael is a perennial crop in India, Sri Lanka, Pakistan, Bangladesh, Myanmar, Thailand, Vietnam, the Philippines, Cambodia, Malaysia, Java, and other southeastern Asian countries [1,2]. Bael is a sacred tree in India. The gardens of many Indian Hindu temples have bael trees [3]. The ripe fruit, which contains a delicious pulp, is the most valuable part of the bael tree [4]. People mainly consume bael as fresh fruit. However, the value-added products of bael, such as drinks, traditional sweets, jam, and pudding, are available in the market [5]. People prefer bael primarily because of its rich taste and ability to cure constipation [6]. Bael fruit is an expensive commodity in supermarkets and street fruit-stalls. All parts of the bael tree possess medicinal values [7]. Thus, bael is famous as a valuable crop with immense medicinal and nutritional potentials [8]. There are many reports available on the medicinal and industrial values of bael in India [9][10][11][12][13]. However, bael is an underutilized crop species in Sri Lanka [14].
Bael can grow well under adverse conditions. Bael trees thrive well in high altitudes as high as 1,200 m and withstand without any significant growth retardation at extreme temperatures such as -7˚C and 50˚C [1]. The prolonged droughts cease fruiting in bael. However, the tree can flourish well under low moisture levels in alkaline, stony, and shallow soils. For example, the 'oolitic-limestone' soils of South Florida provide successful growth and fruit set. Fruitful bael trees dominate many vegetations in arid soils and marginal desert lands of India and Sri Lanka [1].
Bael is the solitary species in the genus Aegle of the subfamily Aurantioideae, which has a well-established phylogeny [15]. However, only the samples collected from India provided the basis for estimating the phylogenetic position, origin, and distribution from the natural range of bael [16,17]. None of the published phylogenetic studies have considered the bael germplasm in Sri Lanka. Thus, the origin of bael germplasm and the genetic diversity of bael germplasm in Sri Lanka are not known.
The Fruit Crop Research and Development Station (FCRDS) of the Sri Lankan Department of Agriculture assessed the bael germplasm in Sri Lanka based on the fruit quality and the information collected from the national network of agricultural extension workers. Thereby, FCRDS identified five elite bael accessions for mass propagation using grafting and micropropagation. However, no details exist on the variation of the fruit size, pulp yield, organoleptic preference, elemental content, and genetic diversity of the elite bael accessions. Therefore, the present study assessed the morphometric variation, organoleptic desire, and elemental composition of the fruits of elite bael accessions in Sri Lanka. We also studied the phylogenetic relationships and molecular dating of the elite bael accessions in comparison to the related germplasm worldwide.

Sample collection
We assessed five elite bael accessions of Sri Lanka in the present study. Table 1 presents the climatic and soil conditions, and Global Positioning System (GPS) coordinates of the original locations of the accessions and their names. We randomly hand-picked fruits from each accession within two-week-window of the peak fruiting seasons at golden-ripe stage (middle of the season: within February to March) of the years 2015, 2016, and 2017. The fruits were placed in cardboard boxes under ambient conditions and transported to the laboratory for measurements. We finished the data collection within 24 hours of harvesting.

Fruit size measurements
The fresh fruit weight (g), length (cm), width (cm) perimeter (cm), and diameter (cm) were measured from 12 fruits per accession in each year. The shell thickness (cm) was obtained from each fruit at six different random places using a Vernier caliper. The pulp, shell, and seed related parameters were measured from five fruits per accession in each year. The pulp (including seeds) was scraped out and collected using a tablespoon from each fruit separately and weighed in grams. The shell weight (g) was also measured. The percentage weight of pulp with seeds, percentage weight of shell, and pulp to shell ratio were calculated. The pulp was excavated using a spatula to pick seeds. The number of seeds per fruit was counted. The weight of seeds (g) per fruit was also measured. The net pulp weight (i.e., pulp weight without seeds), percentage weight of net pulp, and mean-number of fresh fruits at the golden-ripe stage from each accession that is required to harvest 1 kg of net pulp were also calculated.

Fruit size data analysis
The size measurements, shell thickness, and the number of seeds were subjected to the General Linear Model (GLM) Procedure in Statistical Package SAS 9.4 (SAS Institute, Cary, NC, USA, 2013). The significant mean differences were obtained using the Tukey option. The fruit weight (g), length (cm), perimeter (cm) and inner diameter (cm) were normalized to the range of 0-1 and subjected to Principal Component (PC) analysis using the statistical package Minitab 16 (Minitab Inc., USA). The first two PCs were used to draw the PC biplot to depict the variation of fruit size in bael accessions. The fruit pulp, shell, and seed measurements were subjected to the GLM Procedure and Tukey option in SAS.

Organoleptic assessment of ripe bael fruits
We conducted the organoleptic panel study on ripe bael fruits using 30 panelists. Each panelist provided written informed consent before participation in the taste panel. The panel ranked the fruits in a scale of five-levels (1: least preference; 2: low preference; 3: average preference; 4: high preference; 5: highest preference) for the parameters; external appearance, flesh color, aroma, texture, sweetness, and overall preference. The data generated from the taste panel were subjected to FREQ procedure in SAS, and relative preferences were calculated as row percentages for bael accessions. The weighted scores were assigned to each accession for each organoleptic parameter by summing up the products of row percentages and the associated ranks. The weighted scores were subjected to PC analysis in Minitab. The PC biplot was drawn between two major PCs to infer the diversity of accessions based on the organoleptic properties.

Inductively coupled plasma mass spectrometry (ICP-MS) procedure for elemental analysis
Thirty-one elements, including 13 macro and micronutrients in ripe bael pulp samples, were quantified using the ICP-MS method. The pulp extracted from the fruits collected in 2017 was assessed in triplicate for each accession. The fruit pulp was oven-dried in heat resistant crucibles at 80˚C for three-five days until the constant dry weight was reached. This step made samples into hard rock or dried gum like substance. The samples were crushed using mortar and pestle. After crushing, the semi powdered samples were sieved through a fine mesh to obtain finely powdered samples. Approximately 0.2 grams of powdered samples were subjected to digestion using 5 ml of conc. HNO 3 and 2 ml of H 2 O 2 for 30 minutes at 180˚C using Mars-6 Microwave Digester (CEM; Mathews, NC). As a control, 2 ml of deionized water was subjected to the same procedure. The resulted liquid samples after digestion were filtered through 125 μm filter papers. The colorless liquid filtrate samples obtained were diluted up to 50 ml in volumetric flasks. Five milliliters of aliquots of digested samples and the standard comparison sample were then subjected to ICP-MS using Thermo ICapQ analyzer (Thermo-Fisher Scientific Inc., Bremen, Germany). The elemental contents were calculated using the formula given in van de Wiel, (2003) [18]. The elemental data were subjected to GLM procedure in SAS, and the significant mean differences were obtained using the Tukey option. The elemental data were normalized to the range of 0-1 subjected to PC analysis in Minitab. The two first two PCs and three PCs were separately used to draw the PC biplot and PC triplot, respectively.

DNA extraction, PCR and sequencing
The genomic DNA was isolated from immature bael leaves using the Qiagen DNeasy Plant Mini Kit (Cat. No.: 69104). The PCR was carried out in the thermal cycler TP600 (Takara Bio., Otsu Shiga, Japan) for four plant DNA barcoding markers (trnH-psbA, atpB-rbcL spacer, matk-trnT spacer, and trnL) (S1 Table). The forward and reverse primer details and the PCR conditions are also given in S1 Table. The PCR cocktail (15 μl) contained 5× Go Taq Green Master Mix (7.5 μl), 10 μM forward and reverse primers (0.5 μl each), 10 μM spermidine (3.5 μl) [19], and 0.5 μl of template DNA (60 ng/μl). The DNA from the apple variety, Spartan; and distilled water in place of template were used as the positive and negative controls, respectively. The PCR products were size separated using 2% agarose gel electrophoresis and observed for the presence of expected bands (S1 Fig). The PCR products were purified using the PCR Purification Kit, QIAquick (Cat. No.: 28104). The purified PCR products were subjected to 3× cycle sequencing in the Genetic Analyzer 3500 (Cat. No.: 622-0010, Applied Biosystems).

Phylogenetic analysis
The raw outputs from the sequencing were initially visualized in MEGA 7 [20] to determine the initial and end noise. Separate alignments were created employing the Clustal W algorithm [21] for different markers (matK-trnT, atpB-rbcL, trnL-trnF and trnH-psbA intergenic) in MEGA 7 [20]. Then the initial and end noise regions were trimmed, and alignment was checked manually as automated methods could add unwanted INDELs into the alignment. Then the four datasets were combined, and the single alignment was exported to PAUP 4.0a [22]. The UPGMA tree was built employing uncorrected pairwise distances, and the gaps were treated by distributing proportionally to non-ambiguous SNPs so that INDELs can be treated equally. All the substitutions were treated equally. The final dendrogram was further modified using FigTree v1.4.3 [23]. The sequences were adapted from Bayer et al., (2009) [16] for matK-trnT, atpB-rbcL and trnL-trnF intergenic spacers (S3 Table) to determine the phylogenetic position of Sri Lankan A. maemelos of the Rutaceae phylogeny. Initially, the alignments were created separately for each marker, and finally, the combined dataset was created in MEGA 7 [20]. A tree search was carried out in Maximum Likelihood (ML) formwork in RAxML [24] by carrying out a rapid bootstrap algorithm [25]. The algorithm was run for 1000 iterations, and the evolutionary model to compensate the dataset was chosen as GTRGAMMA [26]. Then the best ML tree was used to imprint the bootstrap (bs) values by drawing all the bipartitions into a single tree topology. The final majority rule consensus tree was further modified using FigTree v1.4.3 [23]. The tree was constructed in the Bayesian framework to support the ML tree. Initially, a model selection was carried out in Akaike information criteria (AIC) [27] and corrected Akaike information criteria (AICc) [28] in the J-model test [29] to infer the best nucleotide substitution model to explain our dataset. Then, the best model parameters were used to carry out the tree construction in the Bayesian framework. Using MrBays [30], four Markov chain Monte Carlo (MCMC) were implemented for 50 million iterations to sample the best trees in the tree space. The analysis was set to retain a tree in every 5000 chains, and the initial 10% of the trees were discarded as burn-in. From all the preserved trees, the 50% majority-rule consensus tree was drawn, and posterior probabilities (PP) from each branch were inferred. Finally, the congruence of Bayesian and ML trees was patterned, and the PP values were added to the nodes of the ML tree. The tree constructions and model selection were carried out in the CIPRES Science Gateway [31].

Time calibration
The software package BEAST 2.0 [32] was used for the calibration of the divergence times. Since the best nucleotide substitution model was inferred that fits the sequence datasets of the current study, the rate and shape parameters were implemented to carry out the tree search. A relaxed lognormal clock [33] was used as it can be used to define the substitution rates for each branch. The clock rates were kept being estimated since the divergence dating was done by calibrating the internal nodes of the phylogeny. A mechanistic species level process was used as the branching model because the models are more applicable in divergence dating of sample trees with multiple taxa. Thus, the fossilized birth-death (FBD) model was used as it encounters with lineage diversification utilizing both extinction and speciation [34]. Two lognormal priors were used to calculate birth and death rates. Since the present dataset mostly associated subfamily Aurantioideae, the older nodes (Family: Rutaceae) were not calibrated due to the lack of samples on the species tree.
It is essential to calibrate the nodes that include the whole spectrum of species that can be aggregated into the particular clade. Since there is a considerable amount of data in the citrus clade as well as many studies that have been carried out inferring the divergence time of citrus plants, the species tree of A. marmelos was calibrated mostly using different clades of the citrus crown. Bayer et al., (2009) [16] divided the citrus clade into two groups, i.e., mostly southern clade that includes many plants from the Australian biogeographic zone and Northern clade. Many have carried out fossil calibrations of the citrus plants using different approaches [35][36][37]. Thus, in the present study, we used consensus calibration points according to those studies. The Most Recent Common Ancestor (MRCA) priors were used to calibrate the species tree of A. marmelos. The citrus crown was calibrated using fossil calibration MRCA before an upper boundary of 10.3 million years ago (MYA) and a lower boundary of 4.3 MYA. The Northern and Southern clades were calibrated using fossil calibration MRCA priors with a lower boundary of 3.2 MYA and an upper boundary of 11.6 MYA, and 3.5 MYA lower boundary and 8.4 MYA upper boundary, respectively. In time calibrations of species trees, it is essential to calibrate one basal lineage that connects most of the taxa. Thus, the subfamily Aurantioideae crown was calibrated using MRCA fossil calibration prior to 12.1 MYA lower bound and 28.2 MYA upper boundary. Four MCMC trees were implemented to run a heuristic tree search for 50 million iterations with an initial 10% burn-in. TreeAnotator [32] was used to draw the final Maximum Clade Credibility (MCC) tree. The reliability of the tree search was checked by assessing the Effective sample sizes (ESS) and Autocorrelation time using TRACER v1. 4 [38]. The final MCC tree was further modified using FigTree v1.4.3 [23].

Variation of the fruit size parameters
The weight, length, width, perimeter, inner diameter, and shell thickness of the fruits were not significantly affected by year (P>0.05). The mean fruit weight was significantly highest in PS (followed by RA and MA) in all three years (P<0.05) ( Table 2). Although the mean fruit weight of PS was significantly higher than that of RA, the mean fruit length, width, perimeter, and inner diameter were not significantly different among them (P>0.05). The mean fruit size parameters were markedly lower in BB and PA ( Table 2). The mean shell thickness of BB was the highest (0.4 cm) compared to that of all other four accessions (0.3 cm) (P<0.05) ( Table 2). Furthermore, all the fruit size parameters followed the same trend throughout the collection periods (2015-2018). The highest net pulp percentage was observed in PS, followed by RA, indicating the superiority of these two accessions compared to the other three (Table 3).  Table. The PCA for the fruit size parameters, weight, length, width, and diameter yielded three principal components (PC) to explain the 100% variance (S4 Table). The PC loading status in terms of eigenvalues is shown as a scree plot given in S2A Fig. The weight and length were separately associated with the overall variance of fruit size. In contrast, circumference, width, and diameter got the same effect on fruit size variance because of the near-spherical shape of the bael fruits (S2B Fig). In the PC biplot drawn among PC1 and PC2, two distinct clusters could be seen for the fruit size. One cluster contains PS and RA fruits; however, the relatively larger size of PS fruits is apparent due to the clear separation. The other cluster holds MA, PA, and BB, which have relatively smaller fruits. However, it is evident from the separation that MA fruits are bigger than the BB and PA fruits (Fig 2).

Organoleptic preference on bael fruits/accessions
The organoleptic parameters; external appearance, flesh color, aroma, texture (an indicator of flesh to fiber ratio), sweetness and overall preference (visual appeal and the taste) were significantly associated with the type of accession (χ2 is ranging from 60.0-368.7; P<0.05). The highest preferred external appearance was observed for PS, followed by RA. The PA and BB fruits were not highly preferred (Fig 3A). The most sorted flesh color was reported for RA, and all accessions got ranked above average for the flesh color ( Fig 3B). The most appealing aroma

PLOS ONE
was reported for RA, and like flesh color, all bael accessions were ranked above average for aroma ( Fig 3C). The texture of BB was the least preferred as 40%, and 60% of the panelists ranked BB texture as the least and low, respectively ( Fig 3D). The MA and PA fruits were reported for above-average sweetness, whereas all tasters ranked PS for the average sweetness. The best-preferred sweetness was reported for RA by all tasters (Fig 3E). When considering the overall preference, RA was the highest preferred accession, followed by PS. The accessions MA and PA received the average preference, whereas BB received below average overall preference ( Fig 3F). The ratings received by each accession for all the organoleptic parameters were converted to weighted scores and subjected to PCA which defined four PCs to explain the entire variability of the organoleptic parameters. The PC loading status and Eigenvalues for the weighted scores calculated based on organoleptic parameters are given in S6  4). The PC biplot for the weighted scores of organoleptic parameters provided a clear separation of the bael accessions to depict their strengths for consumer preference. All the organoleptic parameters exhibited their individual effects on the overall variance of organoleptic properties. However, the aroma and sweetness got positioned adjacently in the loading plot indicating their interrelatedness (S3B Fig). The preferred flesh color got placed distantly from the other organoleptic parameters suggesting that it has a peculiar effect on the preference for bael fruits.

Variation of the elemental contents in bael fruit pulp
The elemental profiles revealed by ICP-MS for five bael accessions are given in Table 4. The ICP-MS analysis has not detected Na in bael. The significantly highest Mg content was detected in ripe bael fruits of PA (8881.38 mg/kg), and the least amount was detected in RA (485.83 mg/kg). The significantly highest content of K was detected in MA (15607.70 mg/kg). When the quantitative data of 36 elements were subjected to PCA, a total of 13 PCs resulted. The first and second PCs explained 53%, and the first five PC explained 85% of the total variance (S6 Table; S4A Fig). The PC loading plot for PCA of the elemental contents shows the unique efforts and potential close associations in the elemental contents of the ripe fruit pulp (S4B and S4C Fig). Except for Se isotopes, the effects of all the isotope bearing elements assessed (Pb and Fe), positioned in the same area indicating that their effects on overall diversity are similar or the contents of isotopes of these elements are interrelated. We examined PC biplot and PC triplot to infer the relationships among the bael accessions based on the ICP-MS elemental contents of ripe fruits (Fig 5 and S4C Fig)

Phylogenetic analysis
The UPGMA tree drawn for combined datasets was separated into two clusters at 0.48% separating BB and MA with PA, PS, and RA (Fig 6). BB and MA had 0.094% of uncorrected pairwise genetic divergence, whereas PA diverged out from PS and RA at 0.042% of uncorrected pairwise genetic divergence. PS and RA claded sister to each other, having shallow divergence (with 0.79% uncorrected pairwise genetic divergence). The model selection resulted the GTR +I+G model (AC = 0.850, AG = 1.266, AT = 0.348, CG = 0.688, CT = 0,995, GT = 1.000; proportion of invariable sites = 0.34; Gamma shape parameter = 1.04). The majority-rule  consensus tree built in the ML framework had an almost similar topology with a 50% majority rule consensus tree constructed in the Bayesian framework. The higher bs and PP support values reinforced most of the branches. The nodes that had lower bs were supported with higher PP or vice versa (Fig 7). A highly supported clade was formed with A. marmelos sequences of the present study and the Indian A. marmelos sequence. However, the Sri Lankan bael accessions clustered sister to the Indian A. marmelos. The clade containing the members of A. marmelos clustered sister to Afraegle paniculata, which is found in Africa (bs = 100; PP = 95).

Divergence dating
The ESS for all the priors was above 200, indicating that the probing of the tree-space was done independently by the MCMC chains. In the tree search, the initial 500,000 trees were discarded as burn-in. The final tree had the mean log-likelihood value of -10779. 38 (Fig 8).

Discussion
The FCRDS of Sri Lanka has selected five elite accessions of bael for distribution as cultivars within the country. The accessions BB, PA, MA, and RA were chosen from the Wet Zone, and

PLOS ONE
PS from the Dry Zone (Table 1). The fruit size and pulp parameters did not exhibit significant variation within accession or among the fruiting seasons of each accession (Tables 2 and 3). However, fruit size and pulp parameters were significantly different among accessions indicating higher genetic control on fruit morphometric variation than the effect of environmental factors. The higher genetic regulation of fruit size is a common phenomenon reported for various plant species such as cherry, olive, and tomato [39][40][41]. The variation of the fruit size and pulp parameters indicate that RA and PS are the most industrially desirable cultivars. The probable reason for selecting BB (Beheth Beli means medicinal bael) and PA is the preference shown by the rural people and indigenous medical practitioners towards them. The smallfruited accessions BB and PA are used in local medicinal applications regardless of the smaller fruit size and less preferred taste. The PC biplot fruit size parameters indicated that there are two size-based bael groups present within the elite accessions. However, within the two groups separately, PS and MA possess bigger fruits (Fig 2). The organoleptic assessment also revealed that RA and PS are the two best-preferred accessions with the highest preference ratings (Fig  3). The PC loading plot for weighted organoleptic scores revealed that aroma, sweetness, texture, and external appearance are more related to the overall preference than flesh color. It is

PLOS ONE
also interesting to note that as organoleptic parameters, aroma and sweetness are highly correlated to each other (S3B Fig).
The ICP-MS elemental analysis revealed that bael does not contain sodium. The contents of Mg and K detected in bael samples were comparable with the elemental contents reported in previous studies on bael [42][43][44]. It is apparent from the PC loading plot for ICP-MS-based elemental contents that the elements absorbed/present as groups. The contents of isotopes, Fe, and Pb, are vert tightly correlated to each other (S4B Fig). The absorption of elements to the plants as clusters is a known mechanism in plant species [45,46]. Further studies are needed to ascertain how the different element groups get absorbed into the plants under common regulatory mechanisms. Such studies are vital as the elemental contents within plant tissues are important as nutrients and, in some cases, as toxic substances (exceeding Lethal Dose, 50%) that can cause diseases. The PC biplots drawn for fruit size, taste, and elemental contents revealed a distinct separation of accessions in approximately similar passions (Figs 2, 4 and 5). Thus, the simultaneous selection for the bigger fruit size, better taste, and higher medicinal/ nutritional values decided by the elemental contents could have been the driving force of the domestication of bael.
According to our analysis, the divergence between two genera Afraaegle and Aegle had occurred in the early Miocene (Fig 8). The current understanding of plate tectonics where Indian plate collided with the Eurasian plate which opened the land connectivity from Africa to India in the early Miocene. This land connectivity might have opened the introduction of new plant species through animal movements. One possible hypothesis is the common ancestor of Afraaegle and Aegle could have been introduced to India and Sri Lanka in the early Miocene and evolved to be two separate lineages through isolation. The clustering of Sri Lankan and Indian bael as sister clades (Fig 7) implies that the two sources of germplasm are related to each but not quite the same. Our time analysis suggested that the movement of bael germplasm from India to Sri Lanka occurred in the Pliocene epoch. The land connectivity between Sri Lanka and India in Pliocene ice ages may have aided the germplasm movement. Soon after the Pliocene ice ages, rise of the sea levels might have acted as a robust geographic barrier that restricted possible gene flow between Sri Lankan ad Indian germplasms. As we assessed a selected set of bael accessions in Sri Lanka, it is liable to argue that the naturalization of the bael in Sri Lanka has occurred during the Quaternary period (Fig 8). Where, the slow adoption may have led to the sizeable variability in fruit morphology, taste, and elemental contents.

Conclusions
The elite bael accessions selected for the mass propagation in Sri Lanka exhibits higher diversity in fruit size, taste, and elemental contents. The largest fruits and most preferred taste were reported for the accessions, RA, and PS. Also, on average, 1.3 of PS fruits and 2.1 RA fruits are enough to extract one kilogram of net-pulp, indicating the superior industrial value in making beverage, jam, and sweets. The bael fruit pulp does not contain sodium; however, it is apparent that the elements are absorbed/present as groups. The phylogenetic and divergence dating analyses revealed that bael got introduced to Sri Lanka in the Pliocene epoch. Our analysis also showed that A. marmelos is evolutionary more closely related to genus Afraaegle which is native to African region, showing the ancestry of Aegle could be from African region.