2.1 Plant material

In 2007 and 2008, the Institute for Forestry, Agriculture and Livestock Research (INIFAP) carried out an extensive collection of guava sampling germplasm from commercial plantations, home gardens, and native environments in 11 of the 32 states in Mexico (Padilla-Ramírez & González-Gaona, 2010). A germplasm collection was established at INIFAP's research station in Zacatecas, Mexico (21°44.7′N, 102°58.0′W, 1508 masl) using either cuttings or seeds. The climate type for this location, according to García (2004), is BS1(h’)w, or warm semiarid, with maximum, minimum, and mean annual temperatures of 31.0°C, 9.7°C, and 20.3°C, respectively. Precipitation is between 500 and 550 mm per year, with most of the rain happening between July and September. The soil is characterized by an alkaline pH (8.2), with inorganic nitrogen and organic matter of 28.9 ppm and 1.1%, respectively; potassium and calcium content is high, with values over 2300 and 5800 ppm, respectively (Padilla-Ramírez et al., 2014). Materials in the plantation were arranged in rows separated by 3 m, with a 1.5 m distance between trees. Agronomic practices were as described in Padilla-Ramírez and González-Gaona (2010). Forty-eight materials (Dataset S1) were used in this study to explore the genetic diversity of guava, which included 35 accessions collected in Mexico and 8 from other parts of the world, including Bolivia (n = 1), Brasil (2), Colombia (1), Cuba (1), Honduras (1), India (1), and South Africa (1). The other five remaining accessions of unknown origin were included to maximize phenotypic and genetic diversity in the study. Three of the five accessions with an unknown origin were wild species related to guava, Feijoa spp. (also known as Acca spp.), P. friedrichsthalianum (reports of diploids, tetraploids, and hexaploids; Rojas-Gómez et al., 2020), and Psidium cattleianum (2n = 4x = 44; Souza et al., 2015). The ploidy of these three species was not confirmed. The germplasm collection studied here exhibited substantial variation in fruit size, pulp and skin color, shape, seed content, and metabolites (see below for more details; Figure 1).

2.2 Genotyping

Genotyping was performed using GBS (Elshire et al., 2011) following Diaz-Garcia et al. (2020). DNA was extracted from young and healthy leaf tissue at LANGEBIO, Mexico, using an in-house protocol. Extracted DNA was digested with BglII and DdeI restriction enzymes and sequenced with Illumina (San Diego, CA, USA) NextSeq at LANGEBIO. GBS reads were processed in Tassel 5 GBS v2 Pipeline (Bradbury et al., 2007) using the recently published guava genome assembly as reference for SNP calling (Feng et al., 2021); bowtie2 (Langmead & Salzberg, 2012) was used to align filtered reads (“tags”) against reference genome, and only uniquely mapping reads were kept using samtools (Li et al., 2009). The variant calling file (VCF) was imported into R for further processing and analysis using custom scripts. SNP markers with <5X and >1000X coverage were removed. Moreover, markers with excessive missing data (>10%), minor allele frequency (MAF) < 5%, and with more than two alleles, were also discarded.

2.3 Genetic diversity and selection

Guava genetic diversity and population structure were assessed with complementary methodologies. Genetic relationships among accessions were investigated through principal component analysis (PCA) in R. Similarly, an analysis of population structure was conducted using a Bayesian approach with R package LEA (Frichot & François, 2015), where different numbers of clusters (K) were tested based on a cross-entropy criterion. A phylogenetic tree was constructed for the wild species using Euclidean genetic distances and clustering; the dendrogram was then visualized using R.

To identify putative selective sweeps, a whole-genome scan was performed using the cross-population composite likelihood ratio test or XP-CLR (Chen et al., 2010). Putative selective sweeps were evaluated by comparing wild/native and improved germplasm sets, which are further described in Section 9. A sliding window of 1 Mb and steps of 5 Kb were used for XP-CLR scanning. Because no genetic map is available for guava, genetic positions were approximated by multiplying physical positions by 10e-10, which assumes uniform recombination between markers. The top 5% highest XP-CLR values were considered as putatively selected regions, as in Liu et al. (2020). In addition, F_ST estimates were computed across the genome with the same germplasm groups using the pegas package (Paradis, 2010). XP-CLR and F_ST estimates were visualized in R with the ggplot2 package (Wickham, 2011).

2.4 Phenotypic evaluation

The germplasm collection was evaluated during the 2017 and 2018 seasons for 13 phenotypic traits, including fruit weight, pericarp weight (after removing seeds), polar diameter, equatorial diameter, sepal size, normalized sepal size (sepal size/equatorial diameter), titratable acidity, Brix, number of seeds, total and individual seed weight, pericarp-to-total seed weight ratio, and aspect ratio (polar/equatorial diameter). Harvest and fruit phenotyping took place during September, October, and November. Measurements were conducted on five fruits per accession (to ensure all materials had a consistent number of fruits) in an unreplicated design. Best linear unbiased estimates (BLUEs) were calculated with a linear model in R; only genotype and year effects were included in the model.

3.1 SNP discovery and genotyping in guava

A collection of 48 Psidium accessions, including 45 P. guajava and 3 wild species (Feijoa spp., P. friedrichsthalianum, and P. cattleianum) collected from Mexico and other parts of the world, were genetically characterized using GBS (Figure 2a,b). An average of 2.98 million filtered reads were generated for each accession, resulting in 19,533 unfiltered, unique SNPs (average read depth = 34.80). The genome assembly for the “New Age” cultivar (Feng et al., 2021), widely grown in China, was used as a reference for calling SNPs. Overall, 78.23% of the GBS tags aligned to the reference (58.35% aligned one time and 19.85% aligned >1 time; only uniquely mapping reads were used for SNP calling), which might indicate a recent divergence between American and Asian germplasm. For comparison purposes, the GBS tags were also mapped against the related P. friedrichsthalianum genome (Rojas-Gómez et al., 2022), a diploid Costa Rican native species (~1 Gb), with only a 37.74% mapping rate. Although this low mapping rate might be due to genome misassembly, it also may suggest a closer relationship between American and Asian P. guajava than P. guajava and P. friedrichsthalianum. After filtering sites by coverage (considering SNPs only with >5x and <1000x depth), missing data (<10%), MAF (>5%), and biallelic, only 6712 SNP markers remained for further examination. As expected, on this marker set (which considers all 48 accessions), all three wild species were filtered out, likely because of the genetic dissimilarity between species (Datasets S2 and S3). To further analyze genetic diversity in the wild species, the marker filtering (without applying MAF > 5% because of the small sample size, n = 3) was repeated separately, which resulted in 1155 markers. These 1155 markers included both polymorphic (157) and non-polymorphic markers (998) within the three wild species. However, the 998 non-polymorphic sites in the wild species were polymorphic in the guajava accessions and were useful to assess the relationship between guajava and wild species. From the 1155 markers found across Psidium species, only 297 (25.71%) were also included in the 6712 filtered markers found for the P. guajava samples (Figure 2c). The relationship between Feijoa spp., P. friedrichsthalianum, P. cattleianum, and P. guajava was explored through a phylogenetic tree using the set of 157 markers that were polymorphic in the wild species (Figure 2d); all three species within the Psidium genus were grouped as in previous studies (Flickinger et al., 2020; Vasconcelos et al., 2017). In general, markers were well distributed across the chromosomes in both P. guajava (one SNP per 65.37 Kb; Figure 2e) and wild species (one SNP per 378.57 Kb; Figure 2f).

3.2 Classification of native and improved P. guajava germplasm

In many crop germplasm collections, there is a priori knowledge of the improvement status of the accessions—wild, landrace, and improved/elite (Cao et al., 2019; Diaz-Garcia et al., 2020; Iorizzo et al., 2013; Jeong et al., 2019; Julca et al., 2020; Li et al., 2020; Wu et al., 2018). Guava has a long evolutionary and domestication history; however, the improvement status of much of the cultivated germplasm and germplasm within collections remains unknown or ambiguously assigned. Consequently, the relationships between improved/cultivated and wild/native P. guajava germplasm are poorly understood (Arévalo-Marín et al., 2021). Standard domestication syndrome traits for fruit species—fruit size and sugar/acid balance—can be used to assign an improvement status (e.g., wild vs. cultivated) (Denham et al., 2020). However, guava trees in natural habitats or home gardens of unknown pedigree or origin often have fruit quality and morphological characteristics comparable with improved materials, further complicating the classification of guava germplasm.

The germplasm collection studied here included accessions collected from commercial plantations, natural habitats, and native settings or home gardens (see Dataset S1 for more information). Instead of using this passport information to classify accessions, two groups were formed based on two uncorrelated domestication syndrome traits (Dataset S4), sugar content (Brix) and size (polar diameter; Figure 3a). The two groups of accessions, called herein (1) poor quality (PQ) and (2) good quality (GQ), were defined as follows. The PQ group accessions were defined as those with polar diameter < 58 mm and Brix below a regression line with y-intercept = 14 and slope = −0.05 (Figure 3b). The rationale for the non-zero slope threshold was that early selection efforts in guava might have targeted fruit size first, independently of fruit sweetness as has been observed in other fruit crops (Cao et al., 2019; Guo et al., 2019; Liao et al., 2021). The GQ germplasm was defined as accessions with polar diameter > 58 mm and Brix above the y-intercept = 14 and slope = −0.05 regression line. Consequently, the PQ germplasm had smaller and less sweet fruit, whereas GQ germplasm had larger and sweeter fruit. The number of accessions in the GQ and PQ groups was 15 and 11, respectively. A previous report (Jiménez-Lozano et al., 2009) found that fruit polar diameter on exclusively wild guava germplasm (n = 22) was between 25.8 and 63.8 mm, which is consistent with the threshold used here. In addition, the same study also described Brix to range between 2.4% and 12.4%, which, in general, is the same interval occupied by the PQ group described here. Based on our classification, fruit and pericarp weight, equatorial diameter, and pericarp-to-seed ratio clearly separated PQ from GQ germplasm in the same magnitude and direction as the polar diameter and Brix domestication syndrome traits (Figure 3c). Other traits such as seed number, total seed weight, pericarp-to-seed ratio, and aspect ratio moderately separated both germplasm groups in the same direction (GQ germplasm had greater values). Conversely, sepal size, normalized sepal size, and titratable acidity separated PQ and GQ germplasm in the opposite direction (the PQ group had greater values). The PQ and GQ groups might be associated with wild/native and domesticated/improved germplasm, respectively (as in the same way other diversity studies name wild and domesticated germplasm), particularly because of the overlapping of the fruit quality and domestication domains. Importantly, we are aware of the influence of seasonality on many agronomic traits in fruit crops; therefore, more replication might be required to estimate more accurate genotypic values.

3.3 Genetic diversity and population structure of guava

PCA using all 6712 SNPs in the P. guajava accessions revealed three partially defined groups composed of 6, 8, and 29 accessions—henceforth, Groups A1, A2, and A3, respectively (Figure 4a). Two accessions were not clearly associated with the three groups. Groups A1 and A2 were separated from Group A3 through PC1, whereas all groups were separated through PC2. Additionally, Group A3 represented a large portion of the variation explained by PC1 (i.e., many accessions scattered through PC1). Group A2, which included eight accessions, was composed mostly of GQ germplasm, except for one unclassified accession (however, its polar diameter and Brix were close to thresholds used to define PQ and GQ groups). Both A1 and A3 had a mix of GQ, PQ, and unclassified accessions; however, A3 included a larger number of PQ accessions. Although PC3 explained 6.6% of the variation, this axis did not help identify any additional genetic groups (Figure 4b). PCA of phenotypic data (13 quantitative traits) did not show clear clustering of accessions; however, a gradient through PC1 (70.7% of the variance) and PC2 (28.8% of the variance), with larger fruit in one extreme and smaller fruit in the other, was observed (Figure 4c).

A similar grouping pattern was observed in the structure analysis. Accessions in the GQ and PQ germplasm sets had different admixture profiles (Figure 4d). For K = 2, each germplasm group had a different dominant ancestry, as illustrated in the boxplot of Figure 4d; unclassified materials showed a similar admixture pattern to PQ germplasm. With larger K values (3 and 4), additional groups were observed that did not completely agree with the proposed classification (GQ vs. PQ) based on the phenotypic evaluation. For example, based on K = 4, three admixture patterns appeared: (1) accessions with a clear dominant ancestry (>90%), (2) accessions with a partially dominant ancestry or mixed proportions of the other three ancestries, and (3) accessions with equally mixed proportions of two, three, or four ancestries. Based on a cross-entropy criterion (Figure 4e), the most likely number of groups was four, which might be consistent with Groups A1, A2, and A3 (considering A3 members were scattered through PC1, this group might be split in two) identified through PCA, and that some level resembles the differentiation between PQ and GQ.

3.4 Correspondence between geographic origin, phenotypes, and genetic differentiation

Previous molecular studies of Mexican guava germplasm have reported the lack of association between genetic grouping and geographical origin, mostly when including only local germplasm (Hernández-Delgado et al., 2007; Sánchez-Teyer et al., 2010). Likewise, associations between geographic origin and phenotypic diversity in Brazilian germplasm have not been observed (Fernandes-Santos et al., 2010). Only a few studies have included accessions from different countries, in which case, there are clearer relationships between genetic diversity and geographic origin (Mehmood et al., 2016). One of the problems of working with local, understudied germplasm is the lack of confidence regarding the true origin (e.g., an accession of unknown origin might be introduced to a new location) and nature (e.g., wild, native, and improved) of collected germplasm, which might obscure the true relationship between genetic differentiation and geography.

In our study, the genetic and phenotypic variation between accessions was not clearly associated with geographic origin (i.e., latitude or longitude; Figure 5a,b). Accessions from Nayarit and Sinaloa, two coastal states in the Pacific with similar weather patterns, were separated from the rest of the Mexican accessions, which might suggest a latitudinal gradient of diversity across Mexico; however, larger sample sizes and more intense collection through central Mexico are required to validate this finding. By analyzing groups of accessions based on the state they were collected, molecular markers discriminated slightly better (Figure 5c) compared with phenotypic grouping (Figure 5d); however, PC1–3 using phenotypic data explained close to 100% of the variation compared with only 37% for the molecular markers.

3.5 Candidate selective sweeps for fruit size and Brix

By comparing PQ and GQ germplasm sets, 11 potential XP-CLR selective signals were identified (Figure 6a). XP-CLR significant regions were located in the first seven chromosomes, with sizes ranging from 10 kb to 1.69 Mb (one out of the 11 had just one significant marker; therefore, its size was 0), and a total length of 4.96 Mb (~1% of the genome size). F_ST values were also estimated; however, considering that 21/22 of the significant regions did not overlap with any of the XP-CLR signals, the F_ST analysis was not further discussed. XP-CLR significant regions harbored 117 genes with varying functions. Annotation through Gene Ontology Functional Enrichment Annotation Tool (GO FEAT) (Araujo et al., 2018) identified protein products or gene ontology (GO) terms for 130 genes; 56 genes were reported as uncharacterized proteins or with no associated GO term. Protein functions and GO terms were diverse; therefore, further studies with larger population sizes are required to propose functional categories or gene families under selection. Linkage disequilibrium (LD) decay distance varied considerably between GQ and PQ germplasm (Figure 6b); for PQ germplasm, LD decay to r² = .25 occurred over approximately twice the distance required to achieve the same LD decay in GQ germplasm. LD50, which measures the LD decay distance when the r² drops to half its maximum value, was 1.29 Mb for PQ and 1.21 Mb for GQ. Larger sample sizes for both groups might be required to estimate more accurate LD decay patterns.

Because of the small sample size of this germplasm collection, genome-wide association analysis was not performed. Several QTLs for fruit weight, Brix, and other fruit quality characteristics have been reported previously (Padmakar et al., 2016; Ritter et al., 2010; Valdés-Infante et al., 2003). Comparing fruit size and metabolite QTL locations with genomic regions under selection might help to link putative selective sweeps and phenotypes under selection during guava domestication/improvement. However, transferring QTL locations from old linkage maps into the new guava genome assembly was impossible because no information was available regarding the physical identity (chromosome) of the linkage groups (i.e., no correspondence between linkage group and chromosome numbering) and the lack of marker sequences for anchoring points in common.