Plant material was sourced from a cloned elite population from the New Zealand Radiata Pine Breeding Company (RPBC). We selected a representative population of the third generation of the New Zealand P. radiata breeding program originated from the Elite germplasm collection (Dungey et al., 2009). The study included 622 genotypes from 63 families (10 genotypes per family) resulting from crossing among 55 parents, with an average of 30 clones per family. The P. radiata genotypes selected for this study displayed a wide range of variation in stem growth rate, wood density and pathogen resistance. Three ramets per genotype of similar height and basal diameter were used, with a total of 1,866 plants included in the study. Genotypes were vegetatively propagated from stool beds as bare-rooted cuttings. In August 2015, trees were transferred to 1-L pots with the industry-standard potting mix at 70% water holding capacity (WHC) and acclimated to the polyhouse conditions for over 9 months. The polyhouse was located at Scion’s nursery in Rotorua, New Zealand (Latitude −38°09′28.8′′S, Longitude −176°16′03.3E), and provided natural light but no temperature or humidity control. Plants in the polyhouse were not exposed to precipitation and thus plants were irrigated. The climate in the trial location is temperate, with the warmest months being in summer (January to March), and the coolest months being in winter (June to August). The mean daily maximum temperature is 27°C (summer), while the mean daily minimum temperature is 5°C (winter) (Chappell, 2014).

Pots were positioned within the polyhouse in an incomplete block design. Irrigation was provided by manual watering of pots twice a week, with the irrigation adjusted to maintain potting medium at 70% WHC before the water stress was applied. The experiment was conducted on three biological replicates of each genotype of similar height and basal diameter. Before the commencement of water stress (May 2016), the basal diameter (BD_I) and height (H_I) of each plant were measured. In addition, four fully developed expanded fascicles were collected from under the flushed buds of each ramet. Fresh needles were oven-dried at 65°C and cryo-ground before analyzing carbon isotopic composition (δ¹³C_I) using continuous-flow isotope ratio mass spectrometry (Europa Scientific Ltd., Crewe, United Kingdom) in the Stable Isotope Unit at the University of Waikato (Hamilton, New Zealand). The δ¹³C values (‰) were expressed relative to the Vienna Peedee Belemnite standard (Craig, 1957):

where R_sa and R_st are the ¹³C/¹²C ratios of the sample and the standard, respectively.

In early May 2016, three ramets per genotype of similar height and basal diameter were exposed to drought stress at 22–25% WHC. The amount of water required to reach 22–25% WHC was estimated based on the pot weight (including plant and medium). Over the course of the experiment, changes in pot weight due to tree growth were corrected using a linear relationship between pot weight and tree height. The water stress treatment was maintained by calculating the average weight change of 20–30 pots and compensating water loss every week for 11 months (late March 2017). After 6 months (November 2016) and 10 months (March 2017) of water stress, basal diameter (BD₆ and BD₁₀) and height (H₆ and H₁₀) of each plant were measured. An additional δ¹³C analysis in leaves was performed 10 months after water stress (March 2017, δ¹³C₁₀). The maximum quantum yield of Photosystem II (PSII) in the dark-adapted state (Y) was measured after 5 months (October 2016, Y₅) and 10 months (March 2017, Y₁₀) of water stress using a pulse-amplitude modulated Portable Chlorophyll Fluorometer (Mini-PAM Photosynthesis Yield Analyzer, Heinz-Walz GmbH, Effeltrich, Germany). For each measurement, three ramets per genotype were assessed over 3 days to reduce the effect of the diurnal cycle. On the first day, maximum (F_m) and minimum (F_o) chlorophyll fluorescence emissions of dark-adapted needles were measured randomly in all ramets. During the two consecutive days, F_m and F_o measurements were repeated in a different order. The Y values for each plant were calculated as (F_m − F_o)/F_m and averaged across 3-day measurements. For each ramet, three Dark Leaf Clips (DLC-8) were simultaneously placed on mature, fully elongated needles. Needles were arranged flat across the opening of the leaf clip and the sliding shutter close 20 to 60 min before the measurement to ensure stable readings. Preliminary testing on a plant subset showed that diurnal variation was not significantly influencing the repeated Y measurements within the 3 days.

All replicate plants were carefully dissected after the 11 months of water stress into shoot and root samples. The potting medium was removed gently from the roots, and the root systems were then rinsed in running tap water. Biomass fractions, including shoot dry weight (SDW) and root dry weight (RDW) as well as total biomass or total dry weight (TDW), were determined after oven-drying the samples at 65°C for at least 72 h.

The genomic information for the 622 genotypes was obtained using the exome-capture genotyping by sequencing (GBS) method (Neves et al., 2013). Further details of SNP discovery and capture probe design and testing were previously described in Telfer et al. (2018, 2019). The total number of single nucleotide polymorphisms (SNPs) markers genotyped was 1,371,123. Markers with minor allele frequency (MAF) < 0.01 and missing marker data across genotypes > 40% were discarded from the analysis. The final number of SNP markers used in this study was 61,418. The individual missing SNP genotypes were substituted with the mean genotype for that SNP as implemented in rrBLUP-R package (Endelman, 2011).

Genetic analysis was computed with the Average Information Restricted Maximum Likelihood (REML) algorithm in the ASReml-R v.3 R package (Butler et al., 2009). Single-trait analysis was performed to estimate variance components and heritability for each trait separately, whereas bivariate analysis was performed to estimate genetic correlations among traits. All genetic parameters were estimated and compared using a pedigree-based model best linear unbiased prediction (ABLUP) and genomic-based unbiased prediction (GBLUP).

Preliminary analyses were used for all traits to compare the performance of a spatial first-order autoregressive mixed model and a mixed model without the autoregressive term. The preliminary analyses showed that the spatial mixed model performed better than the mixed model without the autoregressive term (non-spatial) based on Akaike’s information criterion (AIC). Therefore, only spatial mixed models were used in the current study.

The single-trait pedigree-based analysis for all traits was performed using the following model:

where y is the vector of individual tree observations; β is the vector of fixed effects containing the overall mean; a is the random additive genetic effect which was assumed to be normally distributed ∼N(0,A⁢σa2), where σa2 is the additive genetic variance and A is the average numerator relationship matrix; pe is the vector of the random non-additive genetic effect of the individual ramet within the genetic material ∼N(0,I⁢σp⁢e2), where σp⁢e2 is non-additive genetic variance and I is the identity matrix; b is the vector of random incomplete block effect ∼N(0,I⁢σb2), where σb2 are the incomplete block variance and I is identity matrix; X and Z are incidence matrices of the fixed and random effects of the vector of phenotypes y; e is the random residual deviations of individual trees which was partitioned into a spatial (ξ) and independent (η) residuals as σξ2 [AR1 (p_col) ⊗ AR1(p_row)] + I σ² where σξ2 is the spatial residual variance, ⊗ is the Kronecker product, and AR1(p) represents a first-order autoregressive correlation matrix for rows and columns where ρ is the autocorrelation parameter. The independent residual (η) was assumed to be pairwise independent as Iσ², where I is the identity matrix, σ² is uncorrelated residual variance.

A bivariate analysis was employed to estimate the pairwise genetic correlations between each pair of traits. The variance-covariance structure for this model is:

where icb₁ and icb₂ represent the incomplete block effects for 1st and 2nd trait; σb12 and σb22 incomplete block variances for the 1st and 2nd trait, a₁ and a₂ represent the breeding values for the 1st and 2nd trait; σa12 and σa22 are the additive genetic variances for the 1st and 2nd trait; σ_a1a2 is the additive genetic covariance among traits, pe₁ and pe₂ represent the non-additive genetic effects for the 1st and 2nd trait; σp⁢e12 and σp⁢e22 are the non-additive genetic variances for the 1st and 2nd trait, e₁ and e₂ represent the spatially independent residuals; σe12 and σe22 are the independent residual variances for the 1st and 2nd trait and; σ_e1e2 is the residual covariance among traits, and I is the identity matrix. The genetic correlations between traits (r_a) were estimated as

The GBLUP model was the same as the ABLUP model, substituting the average numerator relationship matrix A by the marker-based relationship matrix G implemented based on (Forni et al., 2011) as follows:

Where M is the marker matrix with genotypes coded 0, 1, and 2 for the alternative allele homozygote, heterozygote and reference allele homozygote, respectively, and P is a vector of twice the allele frequency, trace[(M−P)(M−P)′]is a trace of the matrix defined in the nominator and n is the number of markers (Forni et al., 2011).

The individual narrow-sense heritability (h²) is the proportion of phenotypic variance that is explained by the additive genetic variance (Falconer and Mackay, 1996) and was estimated as:

where σa2 is the additive genetic variance, σp⁢e2 is non-additive genetic variance, and σe2 is the spatially independent residual variance. The two heritability estimates are reported in this study regarding implemented methods: ha2, a narrow-sense heritability using variance components estimated in ABLUP and hg2, a narrow-sense heritability using variance components estimated in GBLUP.

The cross-validation was conducted by examining the accuracy of breeding values obtained by the ABLUP (ABLUP-EBV) and GBLUP (GBLUP-EBV) models. The accuracy of both models was evaluated using 10-fold cross-validation with ten random replicates. All data were randomly divided into ten equal groups. In each replication, phenotypes from one group were masked and used as a validation dataset. The accuracies of breeding values obtained by estimating the average Pearson’s correlation between predicted breeding values from masked phenotypes and the breeding values estimated using the full dataset obtained from the ten replicates. These breeding values were obtained from the single-trait analysis (Eq. 2).

Response to selection was calculated as the ratio between selection accuracy and breeding cycle length in years (Chen et al., 2018). Thus relative efficiency (RE) of GBLUP over ABLUP is:

Therefore, the relative efficiency of GBLUP over ABLUP per year is:

Where GBLUP-EBV_V is the genomic breeding values for the trees in the validation dataset, ABLUP-EBV_F is the estimated breeding values from the full dataset, ABLUP-EBV_V is the estimated breeding values for the trees in the validation dataset, L_TS is the length of the breeding cycle for selection based on ABLUP which is approximately 17 years, and L_GS is the length of the breeding cycle for selection based on GBLUP, which is approximately 9 years (Li and Dungey, 2018).

The effective population size Ne (Wright, 1931, 1933) in the current study was 23 and was estimated based on the status number concept (Ns) of Lindgren et al. (1997) as:

where Ns is the status effective number and f is the average coancestry of the population including the coancestry of individuals with themselves.

One-way analyses of variance (ANOVAs) were conducted, after confirming data met assumptions for normality and homogeneous variance, to test whether means for drought-tolerant traits (i.e., δ¹³C and Y) differed between tree families and genotypes. A correlation was also assessed to test for the relationship between δ¹³C and Y after 10 months of water stress.

Trees showed an average 39% increase in the growth of basal diameter and a 63% increase in height growth after 10 months of water stress. The mean δ¹³C value before water stress was −29.67‰, with a range of −26.16‰, to −32.29‰, and a standard deviation of 1.29‰ (Table 1). After 10 months of water stress, the mean δ¹³C value 26.50‰ with a range from −23.27‰ to −31.53‰ (Table 1). The mean Y-value for month five (Y₅) was 0.79 with a range of 0.58–0.87 (Table 1). The mean Y-value at month 10 (Y₁₀) was 0.77 with a range from 0.30 to 0.84 (Table 1).

The measured traits showed a wide range of heritability estimates based on either ABLUP (ha2) and GBLUP (hg2) approaches, ranging from 0.07 to 0.44 (Table 2). Using the ABLUP, ha2 for BD decreased from 0.17 (BD_I) to 0.12 (BD₁₀). In contrast, water stress increased ha2 for height from 0.38 (H_I) to 0.44 (H₁₀). The ha2 for δ¹³C_I was 0.17 compared with 0.07 for δ¹³C₁₀. The ha2 of chlorophyll fluorescence after 5 months of water stress (Y₅) was close to zero due to the small number of observations recorded for this trait (N = 607). However, Y₁₀ was 0.06. The ha2 for dry weight measures (i.e., SDW, RDW, and TDW) were similar of approximately 0.13.

Heritability estimates based on GBLUP were generally lower than the ABLUP estimates except for δ¹³C_I and H_I, which showed slightly higher hg2values than ha2. The greatest difference in hg2 after water stress was recorded for tree height, in which the estimate decreased from 0.41 (H_I) to 0.29 (H₁₀).

The genetic correlation between BD_I and BD₁₀ was very strong (Figure 1). Similarly, a strong genetic correlation was found between H_I and H₁₀ (r = 0.92). There were a high genetic correlations with H and BD at the beginning and end of the experiment. Genetic correlations between H and BD were moderate to large and ranged from 0.48 to 0.90. Moderate to high genetic correlations were found between biomass and tree growth traits, ranging from 0.54 to 0.99. The genetic correlation between δ¹³C_I and δ¹³C₁₀ was 0.73. In addition, genetic correlations between δ¹³C and growth traits were generally positive, with the strongest correlation found between δ¹³C_I and BD_I (r = 0.79). Moderate genetic correlations were observed between δ¹³C and biomass traits, ranging from r = 0.45 to 0.60 regardless of the point in time. The genetic correlation between Y₅ and Y₁₀ was strong (r = 0.90). Furthermore, strong and negative genetic correlations were observed between Y₅ and BD_I (r = −0.92), and between Y₅ and H_I (r = −0.88). Genetic correlations between Y₁₀ and growth traits, and δ¹³C were moderate, ranging from r = −0.41 to −0.64. Genetic correlations between traits based on GBLUP showed a similar pattern but were of weaker strength compared with genetic correlations based on ABLUP (Figure 2).

The cross-validation estimated moderate to large accuracies of breeding values for all traits using all models (Table 3). The accuracy ABLUP-EBV for diameter ranged from 0.77 for BD_I to 0.83 for BD₆, and height ranged from 0.67 for H₁₀ to 0.70 for both H_I, and H₆. The accuracy for biomass traits reached 0.80. The accuracies for carbon isotope δ¹³C were 0.80 for δ¹³C₁₀ and 0.78 for δ¹³C_I, and for Y₁₀ was 0.79. The accuracy of GBLUP-EBV was 7–26% lower than ABLUP-EBV. The accuracy of GBLUP-EBV for diameter ranged from 0.61 for BD_I to 0.65 for BD₁₀, and height ranged from 0.59 for H₁₀ to 0.65 for H₆. The accuracy for biomass traits ranged from 0.66 for RDW to 0.69 for TDW. The accuracies for carbon isotope δ¹³C traits were 0.66 for δ¹³C_I and δ¹³C₁₀, and for Y₁₀ was 0.52.

The RE % of the models showed increased efficiency of GBLUP over ABLUP models (Table 3). For diameter traits, GBLUP showed an increased efficiency that ranged from 75% for BD₆ to 83% for BD₁₀. For height traits, the efficiency ranged from 88% for H₁₀ to 93% for H₆. For carbon isotope δ¹³C traits, the efficiency was 84% for δ¹³C_I and 82% for δ¹³C₁₀, and for biomass traits, the efficiency ranged from 83% for RDW to 87% for TDW. The lowest relative selection efficiency using GBLUP was obtained for Y₁₀ (66%). On the other hand, the relative efficiency in selection per year (RE year⁻¹) resulted in the overall superiority of GBLUP over ABLUP. When assuming the reduction of time for progeny trials from 17 to 9 years by the genomic selection, the lowest relative selection efficiency per year using GBLUP was obtained for Y₁₀.

Significant differences in δ¹³C_I were only observed among families (p < 0.001), but not among genotypes (p = 0.193). In contrast, differences in δ¹³C₁₀ were significantly influenced by both family (p = 0.002) and genotype (p = 0.005). On the other hand, Y₅ differed among families and genotypes (p < 0.001 and p < 0.001, respectively), as well as Y₁₀ (p < 0.001 and p < 0.001, respectively). Nonetheless, no significant linear relationship was found between δ¹³C₁₀ and Y₁₀ (p > 0.05) by family (Figure 3A) or genotype (Figure 3B).

In this study, growth traits such as tree diameter and height presented moderate heritability after water stress. In case of tree diameter, heritability estimates decreased with time after water stress (0.17–0.12). This decrease was mainly due to the increasing residual variance from BD_I to BD₁₀, even though the additive variance increased during the experiment. In contrast, heritability for height increased with time after water stress (0.38–0.44), likely due to increasing the additive genetic variance in balance with increasing residual variance. The observed heritability estimates for BD_I and H_I were in agreement with previous in situ field genetic studies in P. radiata across Australia and New Zealand sites (Kumar et al., 2008; Gapare et al., 2012), which reported heritability estimates between 0.10 and 0.48 for tree diameter, and 0.34 and 0.55 for height. Another study using two populations of water-stressed P. radiata seedlings under glasshouse conditions reported similar initial height estimates but lower initial diameter estimates (Espinoza et al., 2014).

We observed that average δ¹³C increased after water stress and the magnitude of the shift is similar to previous studies with water-stressed P. radiata (Walcroft et al., 1997; Korol et al., 1999; Barbour et al., 2002; Waghorn et al., 2015). The increase in δ¹³C under water stress had overall low to moderate effect on heritability estimates of growth traits before and after water stress, however, the heritability estimate for δ¹³C_I was twice than that for δ¹³C₁₀. This is presumably due to the greater residual variance obtained for δ¹³C₁₀ which was close to fivefold the residual variance for δ¹³C_I. To our knowledge, no published studies investigated estimates for the δ¹³C in genotypes with narrow-sense heritability, such as P. radiata in New Zealand. In other conifers, low to moderate heritability estimates for δ¹³C have been reported; for instance 0.17 for P. pinaster (Brendel et al., 2002) and 0.09 for P. taeda (Baltunis et al., 2008). The decrease in δ¹³C heritability after the application of water stress has been observed in other studies across different species. Xu et al. (2003) found a reduction of between 0.33 and 0.20 in A. cunninghamii growing in wet sites compared to dry sites. Similarly, Ehdaie and Waines (1994) reported a lower δ¹³C heritability estimate in wheat-growing under drought conditions (0.12) compared to wet conditions (0.57). The authors justified this decrease due to the reduced additive genetic variance estimated under water stress conditions. Additionally, Klisz et al. (2019) found loss of any genetic signal in productivity among different Norway spruce provenances when these were growing at marginal climate conditions. Although it is difficult to make direct heritability comparisons across plant species, the results from these and other studies suggest that the induction of moderate to severe water stress may result in the loss of genetic variance in traits.

The estimates of additive genetic variance and narrow-sense heritability showed greater precision when GBLUP was implemented compared to ABLUP. Such a trend was also found in previous studies across different tree species (Gamal El-Dien et al., 2015, 2016; Chen et al., 2018). In summary, when using the marker-based genomic approach, heritability estimates decreased in most of the measured traits compared with the pedigree-based approach. These differences between the two types of correlations may be due to the fundamental differences between genomic and pedigree relationship matrices. For instance, the relationships in the pedigree-based relationship matrix are based on identity by descent (the probability that 2 alleles come from a common ancestor), whereas, the relationships in the genomic matrix reflect the identity by state probabilities (the probability that 2 alleles are the same) (Powell et al., 2010). Furthermore, differences between pedigree- and genomic-based correlations could also be driven by the source of information used. The former uses expected genetic covariation while the latter uses genetic covariation captured by SNPs (Veerkamp et al., 2011; Momen et al., 2017). Another explanation could be that unknown environmental effects were confounding the pedigree relationship-matrix and inflated the genetic variance (Lee et al., 2010).

This study recorded positive genetic correlations among δ¹³C, chlorophyll fluorescence, basal diameter and height growth, and biomass. The moderate to large genetic correlations between δ¹³C and growth traits, including tree height, diameter and biomass, suggest that fast-growing genotypes under water stress could display greater WUE. However, further studies of genetic correlations between δ¹³C and growth traits in P. radiata are warranted to confirm these results.

The findings of this study confirm previous studies with other species that found moderate to strong genetic correlations (Johnsen et al., 1999; Xu et al., 2000; Baltunis et al., 2008). The moderate genetic correlation obtained between the intital diameter or dimater measured after 6 months with initial δ¹³C suggests that diameter might be an indicator trait for greater WUE of P. radiata breeding populations. The positive genetic correlation between growth and δ¹³C_I or δ¹³C₁₀ provides evidence that the variation in δ¹³C (a surrogate of WUE) between genotypes may be possibly due to the differences in the carboxylation efficiency of photosynthesis, rather than changes in stomatal conductance. Further, more drought-tolerant genotypes (i.e., more negative δ¹³C values) may be less sensitive to water stress and able to maintain photosynthetic capacity required for growth with longer periods of open stomata compared with less drought-tolerant genotypes. Rodríguez-Gamir et al. (2019) found that more water-efficient genotypes may be able to maintain a higher whole-plant hydraulic conductance under water stress, thus maintaining leaf water potential and able to keep the stomata opened for longer for photosynthesis. The same study also found that the down-regulation of root aquaporins is an important factor driving intraspecific changes in hydraulic conductance on P. radiata under water stress. Katul et al. (2003) and Rodriguez-Gamir et al. (2020) hypothesized that an equilibrium exists between maximum carbon demand by photosynthesis and maximum water supply, which might result in constant long-term mean intercellular CO₂ concentration. This may be due to the co-evolution of hydraulic and biochemical plant attributes at the genotype level.

The genetic correlations between diameter or height before and after water stress were strong (0.88–0.99), indicating that genotypes were ranked consistently before and after imposing water stress. These strong genetic correlations indicate that the level of water stress in the current experiment was not strong enough to change the ranking of genotypes. This may be due to the experiment not applying severe enough water stress. Indeed, it is possible that genotype rankings could have changed using lower WHC, which would place the genotypes under more severe water stress deficit. Another theory is that only a few genotypes were sensitive to the water stress. However, results indicate that the best performing genotypes for growth traits will remain the best under water stress conditions, which is important for P. radiata breeding programs.

Negative genetic correlations were obtained in this study between δ¹³C and Y. That is, genotypes that performed better under water stress presented more negative foliar δ¹³C values and larger Y values. Physiologically, these results suggest that more drought-tolerant genotypes are more efficient at photosynthesising and have more intact chloroplasts and chlorophyll (i.e., Photosystem II) under water stress, as suggested by recent studies in P. radiata genotypes (Rodríguez-Gamir et al., 2019; Rodriguez-Gamir et al., 2020). Further, Gallart et al. (2019) suggested that the superior growth performance of some P. radiata genotypes under dry field conditions could be due to genotype-specific responses to water stress. However, future research is required to test whether greater drought tolerance also involves greater photosynthetic rates over a wide range of genotypes (Maxwell and Johnson, 2000), to determine the strength of the relationship between Y and δ¹³C.

We did not observe a significant linear relationship between δ¹³C and Y traits. Although previous studies reported δ¹³C to be negatively associated with Y (Mena-Petite et al., 2000, 2005; Dias and Brüggemann, 2010), Y was a poor predictor of water stress tolerance across the 63 families we tested. The poor predictability of Y in our study could be due to several reasons. First, plants were under moderate rather than severe water stress. Unlike other studies, this study maintained the potting medium to a constant 22–25% WHC. The strong relationship between δ¹³C and Y determined by Mena-Petite et al. (2005) was obtained from Y measurements collected over a wide range of foliar relative water potential (∼25–85% range) after withholding irrigation over a short-term period. However, a small number of Y values in this study were below 0.70. Interestingly, our study showed a range of relative water content (%) in needles between 70 and 90% (data not shown), which suggests that the application of more severe and/or more prolonged water stress treatments could have led to stronger Y differences among genotypes. Second, Y is a measure of maximum potential quantum yield of chloroplasts under Photosystem II, which may differ of the actual achievable quantum yield and photosynthesis (e.g., electron transfer rate) under ambient light and water stress conditions (e.g., Subrahmanyam et al., 2006; Baker, 2008; Wu et al., 2008). Although the use of light-adapted chlorophyll fluorescence may be more appropriate to measure quantum yield, the time required to conduct the measurements was considered unsuitable for this study.

The accuracies of GBLUP-EBV for growth traits under the conditions used in this experiment were consistent or greater than the prediction accuracies obtained for growth traits in other tree species, such as Picea glauca (Moench) Voss (Beaulieu et al., 2014; Gamal El-Dien et al., 2016), P. pinaster (Bartholomé et al., 2016; Isik et al., 2016); and P. tadea (Resende et al., 2012). Accuracies of δ¹³C_I and δ¹³C₁₀ could not be compared with literature because this is the first study that evaluates the accuracy of genomic prediction for WUE assessed by δ¹³C for P. radiata. Although the accuracy of estimated GBLUP-EBV was lower than the ABLUP-EBV, the relative efficiency of selection per unit of time (RE year^–1) was greater, reaching half of the breeding cycle (Isik, 2014; Li and Dungey, 2018). This is consistent with the previous studies for height and diameter in P. taeda and Picea abies Karst, in which the relative efficiency of genomic selection per unit of time ranged from 53 to 181% (Resende et al., 2012; Chen et al., 2018).

In the current study, we used the status number concept (Ns) to estimate the effective population size (Ne). Tambarussi et al., 2019 found that the status number is equivalent to the effective population size only when the size of the population is constant over generations.