Genipa L. is a widespread, common Neotropical genus in the coffee family Rubiaceae, tribe Gardenieae DC, consisting of small to large trees, 8–20 m in height, rarely up to 30 m. It is found in a variety of tropical and subtropical lowland (0–900 m) habitats (Pittier and Mell 1931; Steyermark 1972; Burger and Taylor 1993; Zappi et al. 1995). The genus is well known due to its economic and cultural significance.

The most well-known species in the genus is Genipa americana L. It has many uses, for example the fruit is eaten or made into beverages and it is used as a natural blue food colorant. It is important to several indigenous groups who extract an ink from the unripe fruit, which is used as body paint (Steyermark 1972). This practice has been commercialised and the ink is marketed as a henna alternative: so called jagua tattoos. It is also important for its medicinal uses and its timber. The tree is found in gardens, it is cultivated in and around Amazonian villages (Milliken et al. 1992), and it has been proposed as a potential shade tolerant tree crop by the United Nations Conference on Trade and Development (ProFound 2005). Despite its seeming ubiquity, many uses and importance, the taxonomy of this genus is not well resolved.

Previous classifications are based on morphological data and existing phylogenetic studies have been restricted to a few loci and only for one species, G. americana (Persson 2000, 2003; Andersson and Antonelli 2005; Rakotonasolo and Davis 2006; Kainulainen et al. 2013; Kainulainen and Bremer 2014; Mouly et al. 2014; Borges et al. 2021). Genipa has been through several taxonomic expansions and contractions over the years. According to the International Plant Names Index (IPNI 2022) 76 specific names exist in the genus plus a further five infraspecific names. Early circumscriptions of Genipa, for example by Baillon (1880) or Drake del Castillo (1897), encompassed species also occurring in Africa and Asia. Molecular studies demonstrated that this broad circumscription resulted in paraphyly and its size has since been gradually reduced to a solely Neotropical genus (Persson 2000, 2003; Rakotonasolo and Davis 2006). Previous Genipa species have been found to be congeneric with a diversity of Rubiaceae genera including: Agouticarpa C.H.Perss., Aidia Lour., Alibertia A.Rich. ex DC., Benkara Adans., Bertiera Aubl., Burchellia R. Br., Casasia A.Rich., Catunaregam Wolf, Ceriscoides (Hook.f.) Tirveng., Duroia L.f., Gardenia J.Ellis, Glossostipula Lorence, Hyperacanthus E.Mey. ex Bridson, Randia Houst. ex L., Rosenbergiodendron Fagerl., Rothmannia Thunb., Sphinctanthus Benth., and Tocoyena Aubl.

The most detailed taxonomic studies of the genus have been undertaken by Steyermark in the 1970s (Steyermark 1972, 1974) and more recently by Zappi et al. (1995). Few treatments of the genus exist outside of the floras of Central and South American countries (Dwyer 1980; Burger and Taylor 1993; Mendoza et al. 2004; Steyermark and Persson 2004; Delprete and Cortes 2005; Bernal et al. 2019; Gomes 2020). These treatments recognise a different number of species and infraspecific taxa without consensus, summarised in Table 1. Indumentum in particular is considered as an important diagnostic character in Genipa. The different taxonomic hypotheses center on the treatment of G. americana as a single highly phenotypically variable species (Burger and Taylor 1993; Zappi et al. 1995) versus treating its phenotypic variation to be of taxonomic merit (Steyermark and Persson 2004). Steyermark and Persson (2004) recognise two varieties in G. americana and provide a key detailing the morphological differences between the two varieties. G. americana var. americana is distinguished by its glabrous (or nearly so) lower leaf blade and upper leaf surface not rugulose, in contrast with G. americana var. caruto which has a densely soft-pubescent lower leaf blade surface and upper leaf surface sometimes rugulose.

Genipa spruceana Steyerm. was first described by Steyermark (1972) in The Botany of Guayana Highlands. Genipa spruceana is recognized in a number of treatments (Mendoza et al. 2004; Steyermark and Persson 2004; Bernal et al. 2019) but it is not universally recognized, for example Zappi et al. (1995) states that indumentum, the main distinguishing character, are quite variable and perhaps linked to environmental conditions. Genipa caruto Kunth, commonly known as the hairy genip, was first described by Kunth (1820), and since the publication of Flora Brasiliensis (1889) it is often demoted to G. americana var. caruto (Kunth) K. Schum. or not recognized at all (Zappi et al. 1995). The most recently described species in the genus is G. infundibuliformis Zappi and Semir. This species has a more restricted distribution than other members of Genipa, having only been recorded from the Atlantic Forest of Brazil. It is easily distinguished by its distinct flower and leaf morphology.

Table 1.

Summary of Genipa taxa recognized in different works. N/A denotes that it is outside the known distribution of the species.

The global botanical databases and taxonomic data aggregators reflect various taxonomic hypotheses. Kew's The World Checklist of Vascular Plants (2023) lists three accepted species: G. americana, G. infundibuliformis, and G. spruceana. The Missouri Botanical Garden database, Tropicos (2023) lists three species: G. americana, G. chapelieri (A. Rich.) Drake and G. infundibuliformis. In the Tropicos database G. spruceana is treated as a synonym of G. americana. The Leipzig Catalogue of Vascular Plants (Freiberg et al. 2020) and World Flora Online (WFO 2023) list four species in Genipa: G. americana, G. infundibuliformis, G. spruceana, and G. chapelieri. Genipa chapelieri is a Madagascan species (Bridson and Robbrecht 1985) synonymous with G. talangninia (DC.) Drake, recently moved to Hyperacanthus talangninia (DC.) Rakotonas. and A.P. Davis in the Aidia clade (sensu Mouly et al. 2014) and therefore excluded from this study. The Flora e Funga do Brasil (2020) follows the classification of Zappi et al. (1995) and treats G. spruceana as conspecific with G. americana. The entry in the Checklist of the Plants of the Guiana Shield (Funk et al. 2007) is “G. spruceana = G. americana?”, indicating that it is a species of unknown certainty.

The tendency of adopting broad taxon concepts or lumping extends to the infraspecific taxa in Genipa. Five infraspecific names are listed in IPNI (2022) G. americana var. caruto, G. americana f. grandifolia Chodat and Hassl., G. americana f. jorgensenii Steyerm., G. americana f. parvifolia Chodat and Hassl., and G. americana var. riobranquensis Kuhlm. Most of the botanical works listed in Table 1 do not recognize these infraspecific taxa (Dwyer 1980; Burger and Taylor 1993; Zappi et al. 1995; Gomes 2020).

Genipa is widely distributed from Mexico and the Caribbean to Argentina (Fig. 1). The distribution shown in Fig. 1 is based on records from GBIF (GBIF 2020), cleaned using the package CoordinateCleaner (Zizka et al. 2019) in R (R Core Team 2020). Given the variation in the taxonomic treatments of Genipa (Table 1), the distribution of species (Fig. 1) reflects an approximate distribution based on GBIF data (GBIF 2020), as we do not know how the determination of each record was reached. Given the known differences in taxonomic classification in the genus it is likely that G. spruceana, G. americana var. caruto, and G. infundibuliformis are under-recorded and have been recorded as G. americana. Considering the economic and cultural significance of the genus, the distribution shown in Fig. 1 may result from human cultivation; this remains to be tested.

Many habitats where Genipa grows are undergoing drastically increased rates of deforestation and land conversion to agriculture (Hansen et al. 2013). This is especially critical in the Atlantic Forest, where G. infundibuliformis is distributed, where around 85% of the original area has been deforested (Ribeiro et al. 2009). While the IUCN threat status has not been assessed for Genipa in this work, some species are considered to be endangered (G. americana) and vulnerable (G. spruceana) (Ter Steege et al. 2015).

Here, we infer the phylogeny of Genipa using phylogenomic data. Target sequence capture is a genome reduction approach which allows researchers to select and specifically amplify a set of target loci across the genome, using Illumina sequencing (Andermann et al. 2020). This approach has been readily adopted for evolutionary studies as it balances cost, data scale, and computational requirements (Jones and Good 2016; Hale et al. 2020). It is suitable for DNA of limited quality that is more fragmented, such as herbarium specimens or degraded silica dried plant material (Brewer et al. 2019). A major benefit of target sequence capture is the existence of pre-designed bait kits that target known regions of the genome. One such kit is the Angiosperms 353 bait kit which targets 353 single-copy protein-coding genes and works across all angiosperms (Johnson et al. 2019).

Fig. 1.

Putative distribution of Genipa based on cleaned GBIF occurrence records.

We use two phylogenetic inference methods to identify independently evolving lineages from genomic data i) a heuristic two-step approach where gene trees are created first independently and then combined to infer a species tree, and ii) using Bayesian inference where gene trees and the species tree are co-estimated. Both methods implement the multispecies coalescent model (Rannala and Yang 2003; Degnan and Rosenberg 2009; Liu et al. 2009) for phylogeny construction. The model applies probabilistic theory to explain the evolution of alleles and accounts for the incongruence between gene trees and species trees because of incomplete lineage sorting.

Materials and Methods

Taxon Sampling—Twenty-eight Genipa samples (Appendix 1) representing all four putative taxa in the genus were sampled. The taxonomic concept of Genipa in this study follows Steyermark and Persson (2004): G. americana is recognised with two varieties; G. spruceana is treated as a separate species and we follow (Zappi et al. 1995) in the recognition of the species G. infundibuliformis. The samples represent interspecific and infraspecific variation in the genus and comprise of: 11 individuals of G. americana var. americana, eight of G. americana var. caruto, two of G. infundibuliformis, and seven of G. spruceana. Tocoyena pittieri (Standl.) Standl, also in Rubiaceae, was included as the outgroup. All specimens in this study were collected legally and the permits can be presented on request.

Methodology—Total genomic DNA was extracted from silica dried plant material using the NucleoSpin Plant II Kit (Macherey-Nagel, Düren, Germany) or DNeasy Plant Mini Kit (Qiagen, Hilden, Germany). The protocol followed manufacturer's instructions apart from the cell lysis time, which was increased to overnight to maximise DNA yield. DNA quality was assessed using a NanoDrop 2000 spectrophotometer and quantified using the Qubit 2.0. The NanoDrop 2000 and Qubit 2.0 results were used to determine samples that needed concentration by vacuum centrifugation. Gel electrophoresis was also carried out to assess DNA fragment size. Multiple extraction rounds were pooled as necessary when initial DNA quantity was low, in order to meet the minimum concentration requirements of Rapid Genomics, Florida, USA who performed target capture library preparation and sequencing. The minimum DNA sample concentration was 8.52 ng/uL. The DNA was mechanically sheared to a size of 200–500 base pairs (bp). Illumina libraries were constructed and barcode adapters for the Illumina Sequencing platform were ligated to the libraries then PCR-amplified using standard cycling protocols. Samples were pooled into 16 barcoded libraries with equimolar amounts to a total of 500 ng for hybridization. Target enrichment was performed using the Angiosperms 353 bait set (Johnson et al. 2019) targeting 353 putatively orthologous genes. After enrichment, samples were re-amplified for an additional 6–12 PCR cycles and sequenced using an Illumina NovaSeq 6000 with paired-end 250 bp reads.

The Illumina raw read data was processed using the bioinformatic pipeline SECAPR 2.2.5 (Andermann et al. 2018). The bioinformatic pipeline was run on the Sigma2 High-Performance Computing cluster at NTNU, Norway. Raw sequence data was quality checked using FastQC (Andrews 2010) and MultiQC (Ewels et al. 2016) to gain an overview of sequence quality and determine cleaning parameters. Illumina adapters were removed and cleaning of sequences was carried out using FastP 0.23 (Chen et al. 2018). FastP default settings implemented in SECAPR were: i) the read was cut if the accuracy between adapter and read Phred quality score was below 20; ii) the maximum percent of low-quality nucleotides allowed was set to 40 and reads with a higher percentage of unqualified (low quality) nucleotides were discarded; iii) size of sliding window for quality trimming was set to 5 nucleotides; iv) trimming from front and tail if quality value was lower than 10; v) reads below complexity threshold of 10 removed; vi) trim poly repeats at end of read of length 7; vii) low complexity filtering was enabled; and viii) length filtering was disabled. Quality of cleaned reads was checked, using FastQC, MultiQC, and the plotting function in SECAPR.

De novo contig assembly was performed on cleaned reads using Spades 3.15.2 (Bankevich et al. 2012). Overlapping sequences were combined into contig sequences using kmer values 21, 33, 55, 77, 99, and 127. The minimum contig length was set to 200 and contigs under this threshold were discarded. Contigs belonging to target loci were identified by using Blastn (Camacho et al. 2009) to match the contig sequences with a set of reference sequences for each locus. The reference sequences used were the Gardenia philastrei Pierre ex Pit. Davis, A.P. 4055 (K) sequences from the Royal Botanic Gardens Kew PAFTOL project (Baker et al. 2022). A custom target file from Gardenia philastrei, a species close to Genipa, was used as this method has been shown to maximise gene recovery (McLay et al. 2021). A sequence-match was identified if the sequence matched with at least 80% identity across at least 80% of the contig length. Loci with multiple contig matches were discarded as they may represent paralogous sequences. A multiple species alignment (MSA) was created from the contig data using MAFFT 7.490 (Katoh et al. 2019) for each locus that was recovered across at least three samples with the addition of the “no trim” parameter to keep full contig sequence length. Next we repeated the read assembly using the consensus sequence of each locus' MSA as a genus-specific reference library. This additional reference assembly leads in general to a more efficient and less biased retrieval of DNA reads across all samples for each locus (Andermann et al. 2018), as opposed to using the recovered contig sequences for each sample. The minimum coverage parameter was set at four reads. Consensus sequences were generated from the reads mapping to the genus-specific reference at each locus for each sample and from these consensus sequences multiple sequence alignments were computed for each locus using MAFFT 7.490 (Katoh et al. 2019).

Phylogenetic Analysis—Two different phylogenetic methods were used. The first method employed was ASTRAL-III (Zhang et al. 2018), which produces a species tree that shares the maximum number of quartet topologies with the input gene trees and the lengths of the internal branches are inversely proportional to the number of quartets concordant to the split. Gene trees were created using IQ-TREE 2 (Minh et al. 2020). A set of 245 bootstrap consensus maximum likelihood gene trees were created using 1000 bootstrap replicates with UFBoot2 (Hoang et al. 2018) and automatic substitution model selection with ModelFinder (Kalyaanamoorthy et al. 2017) implemented in the IQ-TREE 2 software package. The tree was visualised using Figtree v.1.4.3 (Rambaut 2017).

The second species phylogeny was produced using Bayesian inference, created with Species Tree And Classification Estimation, Yarely (STACEY; Jones 2017) in BEAST2 (Bouckaert et al. 2019) on the CIPRES Science Gateway web portal (Miller et al. 2012). This method simultaneously estimates gene trees and species trees using a birth-death collapse model. The input data was a subset of six loci from the de novo contig assembly dataset. The subset selection consisted of the first six loci in the de novo assembly dataset (5, 9, 20, 43, 55, and 62), with the exception of locus 59, it was excluded from the analysis as it only had seven out of 29 samples. The xml input was generated in BEAUTi 2.6 (Bouckaert et al. 2019). The samples were not preassigned to species and no partitions were selected. The following parameters and priors were selected: species tree model collapse height: 1e^-5; strict clock model: each locus was set as relative to each other; JC69 substitution model; bdcGrowthRate: lognormal (M = 5, S = 2); collapseWeight: beta (alpha = 2, beta = 2); population prior log normal (M = -7, S = 2); relativeDeathRate: beta (alpha = 1, beta = 1). The MCMC was run for 100 million generations and Tracer version v1.7.1 (Rambaut et al. 2018) was used to explore convergence of parameters and effective sample size (ESS). The species tree was generated using TreeAnnotator 2.6.3 (Drummond and Rambaut 2007), after discarding 10% as burn-in, and then visualised using Figtree v. 1.4.3 (Rambaut 2017).

All sequence data generated for this study are available at the GenBank Sequence Read Archive (SRA) under BioProject ID PRJNA1029819 and individual sample accession numbers can be found in Appendix 1. The assemblies and individual gene trees generated in this work are deposited in Dryad (Ridley et al. 2024). The scripts used to create the data in this paper are available on Github at  https://github.com/AntonelliLab/seqcap_processor.

Geographic Distribution—An analysis of the geographic distance between samples was undertaken to determine to what degree geographic isolation assists in the interpretation of our phylogenies. The distance between samples was calculated using the open source GIS software QGIS. Those samples found within 50 km of each other and 51–100 km of each other were noted.

Results

Phylogenomic Analyses—The mean number of raw reads for the samples was 1,126,098, the maximum was 2,183,270 and the minimum 535,602. After cleaning, the average number of raw reads per sample was 1,108,523. The maximum percentage reduction after cleaning was a reduction of 4.48% and the minimum 0.57%. The mean number of de novo contigs that could be identified as being part of the set of target loci was 198 of which 28 loci had matching de novo contig sequences in all samples. An additional step of remapping the reads to a Genipa specific reference library created from the MSAs of the recovered contigs led to the recovery of more loci for more samples. In this approach, sample-specific sequences for 245 loci were recovered on average, 240 of which contained all 29 samples (28 Genipa and one Tocoyena pittieri outgroup). The following alignment summary statistics were calculated using the AMAS tool (Borowiec 2016), mean alignment length was 1508 bp, the maximum was 5674 and the minimum 190, the missing data per alignment mean was 21.55% (min = 2%, max = 96%), the mean proportion of variable sites was 0.12 (min = 0, max = 0.58) and the mean number of parsimony informative sites was 81 (min = 0, max = 764).

The ASTRAL-III phylogeny for the MSA from the reference-based phylogeny containing all 245 loci is shown in Fig. 2. It identified G. infundibuliformis and G. spruceana as fully-supported clades with a local posterior probability (LPP) of one. Within G. americana there are three subclades that are fully supported with a LPP of one: clade A, which contains eight samples: three from Colombia, four from Ecuador, and one from Panama; clade B is comprised of three Bolivian, one Peruvian, and one Colombian sample; and clade C is comprised of six G. americana var. caruto samples: one from Costa Rica, two from Panama, two from Guyana, and one from French Guiana. One G. americana var. caruto from Ecuador is found in clade A and one G. americana var. caruto from Bolivia is found in clade B. The clades within G. americana received full support.

The ESS for all parameters in the STACEY anlaysis was > 241. The STACEY phylogeny, shown in Fig. 3, supports G. infundibuliformis and G. spruceana as monophyletic clades. Three clades are present within G. americana A, B, and C. However, one Peruvian sample (G_am6) was placed within G. americana clade A, whereas in the ASTRAL-III tree it is in clade B. Genipa americana clades B and C received maximum posterior probability scores in STACEY and 0.98 for clade A. The node bars shown on the tree are the height posterior density which represents the 95% central posterior distribution of species tree split times.

Geographic Distribution—The sample distribution shows that several taxa grow in geographical sympatry. Samples of different taxa that were collected within 50 km of each other are: G. americana var. caruto G_car27 and G. spruceana G_spru18; G. spruceana G_spru11 and G. americana G_am9 clade A; G. spruceana G_spru11 and G. americana G_am14 clade A. In addition the sample pairs: G. americana G_am30 clade A and G. americana G_am26 clade B; G. americana G_am8 clade A and G. spruceana G_spru1; G. americana G_am24 clade A and G. americana var. caruto G_car28 were collected within 100 km of each other.

Fig. 2.

Cladogram produced using ASTRAL-III, of 28 Genipa samples, based on 245 nuclear loci, with ASTRAL local posterior branch support shown. Color is used to highlight the three separate clades in G. americana: clade A yellow, clade B pink, and clade C blue. The tip labels show species abbreviation, sample number, and country. Species abbreviations: G_am = G. americana var. americana; G_car = G. americana var. caruto; G_spru = G. spruceana; G_infun = G. infundibuliformis.

Discussion

We produced two well-resolved phylogenies from Angiosperms 353 target capture data using two methods, both based on the multi-species coalescent model that are consistent in topology. The data support the monophyly of G. americana, G. infundibuliformis, and G. spruceana. This study identifies genomic support for the recognition of G. infundibuliformis and G. spruceana as sister species separate from G. americana. We have considered the effect that sampling regimes can have on species delimitation. In this study several sample locations show sufficient distributional overlap to suggest that phylogenetic structure is not solely the result of the sample locations or geographic distance. Our interspecific sampling is especially high in the north of South America where G. americana s.l. and G. spruceana are sympatric. These samples could potentially interbreed in these contact zones but they form independently evolving lineages even when there is sympatry.

Genipa americana, a widespread species distributed from southern Mexico and the Caribbean to northern Argentina, is divided into three well-supported clades. This pattern of separate lineages within one species is common to many species and it is likely to be an understudied but frequent occurrence in the Neotropics (Antonelli et al. 2018; Finch et al. 2022). The phylogenies show that most G. americana var. caruto samples are found in a single clade however two G. americana var. caruto samples fall outside this clade. Samples determined as Genipa americana var. americana are found in clades A and B, and samples determined as G. americana var. caruto are found in clades B and C. This indicates that the current morphological infraspecific classification dividing G. americana into two varieties (Steyermark and Persson 2004) is not supported by the current study. The geographic analysis of sample locations did not find that the distribution distinguished any of the three G. americana clades.

In an attempt to increase taxonomic stability and not add to the already lengthy list of synonyms in this genus, no taxonomic changes are recommended in the genus until diagnostic evidence other than genomic data, such as morphological differentiation, is acquired. Current phylogenomic species delimitation methods do not readily distinguish between population structure and species (Carstens et al. 2013; Sukumaran and Knowles 2017). This can result in taxonomic inflation whereby previously identified infraspecific taxa or new clades are erroneously recognized as new species (Isaac et al. 2004; Sukumaran and Knowles 2017). The use of genomic data as the only means to delimit angiosperm species is not desirable and should only be considered for truly cryptic taxa, as defined by Struck et al. (2018). Here, further evidence is required to determine if the clades in G. americana warrant species status or if the infraspecific rank variety is more appropriate or if they should be recognized at all.

Fig. 3.

Phylogeny from STACEY (BEAST2 plugin) analysis of six locus dataset, the units of branch length are the number of nucleotide substitutions per site, and node bars show 95% height posterior density. Color is used to highlight the three separate clades in G. americana: clade A yellow, clade B pink, and clade C blue. The tip labels show species abbreviation, sample number, and country. Species abbreviations: G_am = G. americana var. americana; G_car = G. americana var. caruto; G_spru = G. spruceana; G_infun = G. infundibuliformis.

Genipa infundibuliformis is fully supported as a separate species in both our phylogenomic analyses. It can be readily determined by its morphology, namely the long corolla tube, reflexed petal lobes, lobed juvenile leaves, and spherical fruit with a smooth surface, which are all distinct characters only found in this species of Genipa. It has a restricted distribution in southeastern Brazil, however, it is not geographically isolated as G. americana s.l. is also present in the Brazilian Atlantic Forest.

Currently G. spruceana is not universally recognized; for example the recent Flora e Funga do Brasil (2020) treats G. spruceana as a synonym of G. americana. This is an example where species circumscription can have considerable consequences for conservation. The taxonomist's decision of broad versus narrow species circumscription can impede efforts to halt biodiversity loss (for example May 1990; Mace 2004; Garnett and Christidis 2017). In this case lumping G. spruceana in G. americana falsely inflates the abundance and possibly the distribution of G. americana, while G. spruceana goes unrecorded, which likely has considerable conservation implications for both species (Bickford et al. 2007; Adams et al. 2014).

In order to ensure that species in the genus are diagnostic beyond the genomic level, more field studies are required, particularly in Brazil, Guatemala, and Paraguay as it would provide further morphological and ecological data. The results of this phylogeny can be combined with additional lines of evidence such as morphology in an integrated approach. This may elucidate diagnostic features for the three clades in G. americana. Once the above avenues are investigated a decision can be made on the taxonomic rank applicable (if any) to the clades within G. americana. By applying the multi-species coalescence model to detect independently evolving lineages in Genipa, we show support for three species and evidence of infraspecific genomic structure within G. americana s.l. A stable systematic framework for Genipa based on an integrative taxonomy approach is important for conservation of species in areas undergoing unprecedented rates of habitat modification, putting species and its interspecific variation at risk of extinction.