A total of 104 avocado “criollo” “plus trees” were selected in different remote agroecological zones up to 2,400 m asl of the province of Antioquia, northwest Andes of Colombia, during the years 2014 and 2015 (originally enlisted under the first table of Cañas-Gutiérrez et al. (2019b), and summarized here under Supplementary Table 1). Their selection was made based on criteria, such as adaptability to their growing area, longevity, health, and productivity. During this period, trees were in situ characterized by taking morphological data, and carrying out their molecular fingerprinting (section below on “microsatellite markers characterization”). Collection of yield data was additionally carried out in these trees during four consecutive years (2014–2017).

The seeds produced by selected elite avocado “criollo” trees were germinated at greenhouses, and obtained 3-month-old seedlings were grafted with cv. Hass scions. In the year 2015, the cv. Hass grafts were planted at a common garden (i.e., provenance test) at AGROSAVIA’s research station La Selva, located in Rionegro, province of Antioquia, at 2,100 m asl. Morphological data, such as length and diameter of the rootstock, the canopy length and diameter of the scion, and the morphology of the graft scar were collected during the years 2017–2021. The production data of the cv. Hass grafts were collected in the same period from 2017 to 2021. Since half-sib seedling rootstocks grafted with the Hass clone have been in production during five consecutive years in the same environment (La Selva locality), we were able to minimize environmental variance and equal phenotypic variance to the rootstocks’ additive genetic component, which is the main benefit of carrying out common garden or provenance tests. These trials intend to neutralize environmental effects over phenotypic expression, by planting diverse genotypes in a single locality (Lascoux et al., 2016). Our setup has allowed heritability analyses to include a total of 9 years of production data, both for in situ “criollo” “plus trees” (2014–2017), as well as in their OP half-sib seedling progenies used as rootstocks for the Hass scion (2017–2021).

A total of 17 morphological and agronomic traits were measured during the sampling trips of 2014 and 2015 in the avocado “criollo” selected “plus trees,” according to descriptors that are specific for avocado based on the International Plant Genetic Resources Institute (IPGRI). These descriptors included tree, leaf, fruit, and seed traits (Table 1 and Supplementary Table 2). Fruit production was recorded in these trees during four consecutive years until 2017. Pearson’s, Spearman’s and Kendall’s correlation matrices (Figure 1) were visually inspected using customized plots drawn with the cor.test and heatmap functions in R v.3.4.4 (R Core Team). Complete phenotypic measures were possible for a total of 82 “criollo” “plus trees” from the original sampling.

Furthermore, a total of 22 OP half-sib seedling progenies of 82 elite avocado “criollo” “plus trees” were used as rootstocks for the Hass scion. Grafted trees were characterized at AGROSAVIA’s La Selva research station in terms of their early growth traits (i.e., young shoot color, stem color, stem lenticels, lenticels color, average number of leaves at transplant, and average stem diameter at grafting), as well as in terms of their production (i.e., cumulative total number of fruits and fruit’s weight) 5 years after grafting (2017–2021). The analyses focused on yield traits because they allowed for a direct comparison with the recorded production traits in the “criollo” “plus trees,” besides the fact that yield traits are of paramount interest for commercial avocado cv. Hass orchards, especially those aiming for exportation markets.

Total genomic DNA was obtained from avocado “criollo” leaves based on the extraction method standardized by Cañas-Gutiérrez et al. (2015). Thirteen microsatellites [simple sequence repeats (SSRs)], originally designed by Sharon et al. (1997) and Ashworth et al. (2004), were chosen (Supplementary Table 3) for their high polymorphism information content following estimates by Alcaraz and Hormaza (2007). Fingerprinting highly polymorphic SSR markers allowed disclosing pedigree-free marker-based relationships among genotypes, which in turn enabled distinguishing from highly related to completely unrelated individuals. Such contrasting scales of SSR-reconstructed shared relatedness due to recent co-ancestry (i.e., isolation by descendent, IBD) offer the basis to quantify phenotypic clustering patterns across family types (i.e., genetic heritability), as demonstrated in oil palm (Cros et al., 2015) and rubber tree (Cros et al., 2019), and as expanded in the section below.

Three multiplex PCR amplifications were performed in 10 μl volume containing 16 mM (NH4)2SO4, 67 mM Tris-HCl pH 8.8, 0.01% Tween20, 2 mM MgCl2, 0.1 mM each dNTP,0.4 μM of each primer, 25 ng genomic DNA and GoTaq^® ADN polymerase (Promega, WI, United States). Forward primers were labeled with WellRed fluorescent dyes on the 5′ end (Proligo, France). The reactions were carried out in an I-cycler (Bio-Rad Laboratories, Hercules, CA, United States) thermo cycler using the following temperature profile for multiplex PCR 2 and 3: an initial step of 1 min at 94°C, 35 cycles of 30 s at 94°C, 30 s at 50°C, and 1 min at 72°C, and a final step of 5 min at 72°C. For multiplex PCR 1, the same temperature profiles were used, except for 55°C for annealing temperature that was used instead. The PCR products were analyzed by capillary electrophoresis in the equipment ABI PRISM^® 3130 Genetic Analyzer (Applied Biosystem, CA, United States). Allele sizes were estimated in base pairs with a Peak Scanner (Thermo Fisher Scientific, United States), allowing for a maximum of two alleles per sample.

We implemented a G-BLUP (genomic best linear unbiased predictor) model to quantify genetic parameters (i.e., narrow-sense heritability and breeding values) for key agronomic traits across the avocado “criollo” “plus trees” selected as seed donors for the OP half-sib progenies used as rootstocks for the Hass clone. Genetic parameters are a baseline of any breeding program (Holland et al., 2003) since they guarantee that additive genetic gains are maximized, while breeding cycles are minimized (Dieters et al., 1995). This is because genetic parameters account for the proportion of phenotypic variance among individuals in a population due to genetic effects (Milner et al., 2000; Kruuk, 2004; Berenos et al., 2014). In this sense, the G-BLUP precisely offers a modern pedigree-free marker-based approach (Meuwissen et al., 2001; Crossa et al., 2017) to estimate genetic parameters on populations of mixed ancestry (Frentiu et al., 2008; Wilson et al., 2010), such as OP seedlings.

We first computed a G relatedness matrix, following Lynch and Ritland (1999), to account for the half-sib structure of the OP seedling progenies used as rootstocks. To validate relationships within the G matrix we further computed (1) a distance tree by means of the Neighbor Joining algorithm (Kamvar et al., 2014) in R (R Core Team) with 10,000 bootstrap replicates, as well as (2) an unsupervised Bayesian clustering using the STRUCTURE software (Pritchard et al., 2000) with five independent runs for each K-value from K = 2 to K = 8 and 100,000 Monte Carlo Markov chain replicates with a burn-in of 50,000.

The admixed origin of the avocado “criollo” genepool, due to avocado’s protogynous dichogamy, is known to favor an adequately variable G matrix, embracing various degrees of relationship, as much as the Avocado Genebank (Cañas-Gutiérrez et al., 2019b). The variability in the G matrix (from highly related to completely unrelated trees) enabled the use of the G-BLUP predictor. Hence, we relied on this G relationship matrix to compute the genetic-based G-BLUP model, following Arenas et al. (2021) and Fernandez-Paz et al. (2021), as in:

where y corresponds to the phenotypic trait vector, a is the vector of individual random additive genetic effects with a normal distribution [in other words, a ∼ N (0, G σ_a²), where G is the among-tree pedigree-free SSR-based relationship matrix, and σ_a² corresponds to the additive genetic variance (Müller et al., 2017)], b is the vector of fixed effects (i.e., general mean, intercept), e is the vector of residual effects, and X and Z are the incidence matrices for fixed effects and additive genetic effects (Chen et al., 2018; Gutierrez et al., 2018).

Parameter estimation for the G-BLUP model was carried out using reproducing kernel Hilbert space (RKHS) implemented in the BGLR (Bayesian Generalized Linear Regression) software (Perez and De Los Campos, 2014) under an R v.3.4.4 environment (R Core Team). The multi-dimensional space of parameters was sampled with 10,000 iterations, an initial burn-in of 5,000 steps, and a thinning interval of 10 for data recording. Trace plots were drawn in R to verify convergence in the posterior distributions.

Narrow sense rootstock-mediated heritability (h²), which is the proportion of the overall phenotypic variance accounted for additive genetic variance, was then calculated by retrieving the additive (σ_a²) and residual (σ_e²) variances from the RKHS algorithm, following de los Campos et al. (2015):

The vector of h² scores was summarized using the median and the 95% CI of the BGLR’s posterior distribution. Statistical accuracy (a.k.a., prediction ability, r_y) was computed per trait using Pearson’s correlation between observed and predicted (μ + GEBV) trait values (Müller et al., 2017; Zhang et al., 2019). Predicted trait values were recorded and subtracted to the overall phenotypic mean (μ), to gather the vector of genetic-estimated breeding values (GEBVs), which is the deviation from the overall trait mean attributed to additive genetic effects. GEBVs were summarized throughout a biplot diagram using R’s princomp and biplot functions for the principal components analysis (PCA).

Finally, we implemented mother-offspring regressions to estimate h² scores for key agronomic traits in the OP half-sib seedling progenies of avocado “criollo” “plus trees” grafted with the Hass scion. This heritability analysis included 5 years of data in the OP half-sib seedling progenies of avocado “criollo” “plus trees” used as rootstocks for the Hass scion via the mother vs. half-sib regression analysis (Falconer et al., 1996). These h² estimates were gathered for cumulative number of fruits and weight of fruits from the Hass clone grafted on top of the OP half-sib seedling rootstocks derived from the selected avocado “criollo” “plus trees.” Mother-offspring regressions were carried out in R using the function lm. Regressions were plotted for validation in ggplot, under the same R environment. The narrow-sense heritability score was calculated as twice the slope from the mother-offspring regressions, following Lynch and Walsh (1998).

The 13 microsatellite markers that were used for the molecular characterization of 104 avocado “criollo” “plus trees” amplified a total of 147 allelic fragments. This indicated sufficient polymorphism to reconstruct the pedigree-free marker-based relationship matrix (Figure 2), while distinguishing highly related from completely unrelated admixed individuals (Figure 3). The phylogenetic tree (Supplementary Figure 1) corroborated the latter point. Therefore, the segregation of 147 allelic SSR fragments allowed to effectively quantify phenotypic clustering patterns across family types via de G-BLUP model (i.e., SSR-based genetic heritability estimates, Cros et al., 2015, 2019) for key morphological and agronomic traits in “criollo” avocado “plus trees.”

The heritability parameters for a total of 11 morphological traits (Figure 4) and 6 agronomical traits (Figure 5) were computed in 82 of the 104 avocado “criollo” “plus trees” selected as seed donors for rootstocks of the cv. Hass. Estimates of avocado “criollo” “plus trees” narrow-sense heritabilities (h²) were significant for the 11 measured morphological traits and for the 6 agronomic traits, according to the 95% CIs from the posterior distribution, all of which were above absolute zero (Table 2).

The estimated heritabilities (h²) for the 11 morphological traits ranged from 0.28 to 0.51. The tree shape and branching traits presented the h² lowest scores, 0.28 and 0.29, respectively. The highest h² values were presented for ribbed petiole and adaxial veins with 0.47 and 0.51, respectively. The 6 agronomic traits evaluated presented very similar h² values that ranged from 0.34 to 0.39. Seed weight was the trait with the lowest h² score (0.34), while fruit weight and total volume traits exhibited the highest values, both with 0.39 (Table 2).

Confidence intervals for the heritability scores were generally narrow for the 17 examined traits. The traits that presented the widest confidence intervals were: number main veins (0.1–0.7, r = 0.64) pulp weight (0.1–0.7, r = 0.66), adaxial veins (0.2–0.8, r = 0.78), and ribbed petiole (0.2–0.8, r = 0.75). Meanwhile, the trait that presented the narrowest confidence interval was the shape of the tree with 0.1–0.4 and r = 0.69 (Table 2). Adaxial veins and ripped petiole traits exhibited the widest 95% CIs, as well as the highest heritability scores (0.51 and 0.47, respectively). On the contrary, the tree shape trait had the lowest 95% CI, and also showed the lowest heritability score (0.28). Broadly speaking, the morphological and agronomic traits measured in avocado “criollo” “plus trees” selected as seed donors for rootstocks of the cv. Hass presented heritabilities in a medium range. This suggested a relative contribution of additive genetic effects in the total variation of these traits, which may enable moderate genetic gains via recurrent selection.

To assist the ranking of avocado “criollo” “plus trees” as seed donors for OP half-sib seedling rootstocks of the cv. Hass, we summarized the breeding values for the 11 morphological and 6 agronomical traits across 63 of the 104 avocado “criollo” “plus trees” selected as seed donors for rootstocks of the cv. Hass using a biplot diagram from a PCA (Figure 6). The length of each vector in the biplot diagram allowed ranking the traits by their additive variance. In this sense, fruit diameter, total volume, fruit weight, pulp weight, fruit length, and branching exhibited the most notorious additive genetic contributions, reinforcing the heritability results described in the previous section. On the contrary, the ribbed petiole and trunk circumference traits, showed lower vector lengths, indicating a lower additive genetic variance for these traits. These results suggested that fruit traits tended to exhibit more potential to respond to selection, indicating that production at commercial avocado cv. Hass orchards may be leveraged by a proper selection of avocado “criollo” “plus trees” as seedling rootstock donors, based on their breeding values.

The biplot diagram also allowed pinpointing the elite avocado “criollo” trees that may be promising for selection given their heritable superiority in terms of agronomical traits. For instance, the trait vector for fruit diameter (associated with the fruits’ caliber, an important factor for its commercialization) pointed toward “criollo” trees ANSS9 and ANNE55. Similarly, the trait vectors for total fruits’ volume, fruits’ length, and fruits’ weight, respectively pointed toward the avocado “criollo” “plus trees” encoded as ANDA62, ANCI75, and ANCH65. Meanwhile, the avocado “criollo” tree encoded as ANTU57 was closed to the vectors for the traits fruit diameter and pulp weight, which also makes ANTU57 a promising donor for seedling rootstock families (Figure 6). On the other hand, the morphological traits branching, leaf shape, trunk surface, tree shape, and adaxial veins presented the longest vectors, and consequently greater potential for selection. Some of these morphological traits may be of interest for half-sib rootstock family selection. For example, branching and tree shape may contribute desired tree architectures optimum for light penetration to boost flowering and harvesting.

A total of 22 “criollo” “plus trees” were selected as seed donors for rootstock establishment according to the trait criteria enlisted in the previous section. Their corresponding OP half-sib seedling progenies were grafted with Hass clonal scion to estimate rootstock-transferred h² scores for key agronomic traits, such as total number of fruits and fruits’ weight after 5 years of production. The estimation of overall h² scores (spanning both the mother-offspring and the rootstock-mediated inheritance) for both agronomical traits followed a mothers vs. grafted half-sib offspring linear regression. In this sense, the dependent variable corresponded to the average half-sib family trait values for total number of fruits and fruits’ weight measured in avocado cv. Hass trees grafted on top of OP seedling rootstock families obtained from the avocado “criollo” “plus trees.” The trait values for the total number of fruits and fruits’ weight measured in the original mother “criollo” “plus trees” were depicted in the X-axis (Figure 7).

The trait total number of fruits presented a slope of 0.18 and an r score of 0.16. Since the mothers-offspring heritability corresponded to twice the slope, the h² equaled 0.36 (± 0.23). The trait fruits’ weight exhibited a slope value of 0.056 with an r score of 0.18, and therefore an overall h² estimate of 0.11 (± 0.09, Figure 7). According to the results obtained, the trait total number of fruits presented a moderate seedling rootstock-mediated heritability of 0.36, indicating the possibility of harnessing natural variation from “criollo” plus trees via OP half-sib seedling rootstock families for cv. Hass orchards. On the contrary, the trait fruits’ weight presented a lower seedling rootstock-mediated heritability of 0.11, indicating that this trait may be influenced by environmental factors, non-additive genetic factors, such as dominance and epitasis, crop management, and the rootstocks-scion interaction, rather than additive genetic contributions from the OP half-sib rootstock families.

A long-standing debate in the breeding of fruit tree crops is whether seedling rootstocks are capable to enhance productivity in the absence of more uniform clonal rootstocks, which in turn are thought to lack sufficient genetic variability for long-term adaptation. As an intermediate solution, replacing seedlings genotypes by a diverse panel of elite clones is expected to convey major improvements while controlling adaptive genetic erosion (Ingvarsson and Dahlberg, 2019). Nonetheless, avocado cv. Hass plantations in northwest South America have grown exponentially in the last decade without the availability of clonal rootstock genotypes. Instead, they have had to rely on cheaper seedling rootstocks from native “criollo” trees without any sort of selection or traceability. Therefore, our narrow sense heritability quantifications are encouraging guides to additive genetic selection of “criollo” “plus trees” as seed donors for rootstock production, conveying up to 35% (h²_T) of potential genetic gain for a key yield trait as is the total number of fruits. This overall effect totalizes the additive genetic contribution from OP half-sib seedling rootstock families to the Hass clonal scion (i.e., h²_R = 0.44, as computed by Reyes-Herrera et al., 2020), as well as the expected additive genetic variance putatively inherited from the “criollo” “plus trees” to their OP half-sib F1 seedling offspring (i.e., h²_F1 = h²_T/h²_R = 0.7). Ultimately, these estimates reinforce the feasibility of utilizing selected avocado “criollo” “plus trees” as seed donors for rootstocks.

Furthermore, the GEBVs, computed here by weighting and standardizing the overall trait superiority by the corresponding heritability scores, are already a robust selection index to target in situ superior “criollo” “plus trees” as donors of seedling rootstocks. Alternatively, buds from these selected “criollo” “plus trees” may be grafted to establish ex situ seed orchards. In any case, avocado cv. Hass nurseries may still face uncertainty while sorting the OP half-sib offspring from those “criollo” “plus trees” families in the absence of age-age correlations at the seedling stage. To assist early within-family selection of superior half-sib progenies, heritability estimates above 35% (e.g., for total number of fruits) also invite implementing last-generation high-throughput genomic selection (GS) schemes (Desta and Ortiz, 2014; Crossa et al., 2017). Genomic prediction platforms calibrate infinitesimal polygenic additive models (Cortés et al., 2020) capable to forecast breeding values for low-heritability complex traits (Klápště et al., 2020) at early life stages (Grattapaglia et al., 2018). Therefore, GS simultaneously makes indirect selection more precise, while allowing speeding up breeding cycles. Its potential in perennials with prolonged juvenile phases has already been industrially scaled up by the forest tree sector (Resende et al., 2012, 2020; Tan et al., 2017; Cappa et al., 2019; Thistlethwaite et al., 2019), with some prospective pilot studies carried out in temperate fruit trees (Kumar et al., 2012, 2019, 2020; Muranty et al., 2015; Iwata et al., 2016; Ferrão et al., 2021). Based on this, avocado cv. Hass nurseries might exploit the potential of GS to choose the superior progenies from the “criollo” “plus trees,” even before grafting. Such innovation will optimize timing and monetary resources at nurseries, while ultimately translating into reduced tree mortality and increased productivity/adaptability during the establishment and production phases.

Meanwhile, efforts must not be spared to breed a diverse panel of elite clonal rootstock genotypes resistant to root rot (P. cinnamomi), a major fungal threat for avocado plantation worldwide (Toapanta-Gallegos et al., 2017; Sánchez-González et al., 2019), while keeping high productivity and adaptation across heterogeneous mountain climates (Cortés and Wheeler, 2018). However, breeding elite clonal tree genotypes with conventional direct phenotypic selection doubles the breeding cycle length compared to gradual population recurrent selection due to longer progeny-testing and more clonal trials (Neale and Kremer, 2011). Therefore, GS may once more assist the identification of elite clonal rootstocks for the avocado cv. Hass industry, as already exemplified by the forestry sector by reducing the length of the clonal breeding cycle to 50% (Resende et al., 2012). Parallel efforts must be conveyed to better standardize the propagation of clonal rootstocks [e.g., via genotype-dependent micro-cloning (Ernst, 1999) and double grafting (Frolich and Platt, 1971) protocols], which is still a major bottleneck in regions with high availability of native avocado trees capable to source nurseries with cheaper seedlings.

The h² and GEBVs estimate gathered as part of this study encourage the implementation of more systematic genetic-guided indirect selection platforms both as part of early (before grafting) seedling screening at nurseries of avocado cv. Hass, as well as within the avocado rootstock breeding cycles per se. Nonetheless, a question persists, which is whether the studied SSR molecular markers are enough to capture sufficient additive variance as to make predictive breeding reliable. This potentially remains open judging by the r scores (> 0.54) of our pedigree-free SSR-inferred estimation of genetic parameters (Frentiu et al., 2008; Wilson et al., 2010), which embraced sufficient variation across distinct kinship levels (Milner et al., 2000; Kruuk, 2004; Berenos et al., 2014). Therefore, avocado rootstock breeders would be able to leverage SSR characterizations across a diverse panel of avocado accessions (Alcaraz and Hormaza, 2007; Cañas-Gutiérrez et al., 2015; Ferrer-Pereira et al., 2017; Boza et al., 2018; Cañas-Gutiérrez et al., 2019b; Reyes-Herrera et al., 2020; Sánchez-González et al., 2020). SSR-inferred relationships also allow bridging the difficulty to generate complete pedigrees in tropical avocados. It is challenging to carry out control crossing in tropical “criollo” avocados because their pollination is open, and exhibit dichogamous protogyny. It would also be unviable to gather complete pedigrees from in situ “criollo” avocado trees because they exhibit disparity in the phenological phases, and are situated in remote locations where access is limited to perform controlled monitoring. This is why our team is currently bringing buds from the “criollo” trees for grafting on their own half-sib rootstocks. This strategy would allow cloning them at a single research station to better study their phenology and control their pollination. Meanwhile, molecular markers offer a realistic pathway to reconstruct relatedness among avocado “criollo” trees. The prospect to utilize highly polymorphic SSRs for predictive breeding across a diverse panel of genotypes spanning various levels of relatedness is also in line with previous proof-of-concept pilot studies carried out in oil palm (Cros et al., 2015) and rubber tree (Cros et al., 2019). This is because indirect selection works either on the basis of shared relatedness (typically measured as kinship relationships due to recent co-ancestry, Sedlacek et al., 2016) or on the basis of linkage disequilibrium (LD) between marker loci and the genetic variants that underlie phenotypic variation (Thistlethwaite et al., 2020). The high mutation rate of SSRs (Ellegren, 2004), and consequently, its unsurpassed polymorphism information content (Cortés et al., 2011; Blair et al., 2012), efficiently disclose the shared co-ancestry among genotypes (Lynch and Ritland, 1999) and its ultimate potential for predictive breeding.

Nonetheless, predictive breeding would be incapable to rely on LD (Blair et al., 2018) between markers and the genetic variants that underlie trait variation because SSRs are insufficient to reconstruct the polygenic architecture of complex (Hirschhorn and Daly, 2005) rootstock-mediated traits. After all, SSR markers, despite polymorphic, are rare across the genome, which makes them unlikely to be found in LD with any causal variant (Slatkin, 2008). To bridge this gap, future research will have to couple SNP screening of avocado genotypes (Kuhn et al., 2019; Rendón-Anaya et al., 2019; Rubinstein et al., 2019; Talavera et al., 2019), via genotyping-by-sequencing (GBS, Elshire et al., 2011; Cortés and Blair, 2018), whole-genome re-sequencing (WGR, Fuentes-Pardo and Ruzzante, 2017) and RNAseq (Jensen et al., 2012; Sun, 2012; Reeksting et al., 2016), with genome-wide association studies (GWAS, Khan and Korban, 2012). These innovative “omic” approaches should also target diverse avocado “criollo” “plus trees,” their OP half-sib progenies used as rootstocks for the Hass scion, and different tissues of the grafted tree, for instance through single-cell sequencing (Tang et al., 2019). Combined genetic mapping efforts will in turn inform on the upstream regulatory gene networks that contribute shaping the rootstock—scion compatibility (Loupit and Cookson, 2020) and interaction (Albacete et al., 2015; Warschefsky et al., 2016; Migicovsky and Myles, 2017), as well as on the corresponding downstream mechanistic physiological processes [e.g., transport of water and nutrients from the rootstocks, and large-scale movement of hormones, proteins, mRNAs, and sRNAs between the rootstock and the scion (Wang et al., 2017; Cookson et al., 2019)].