Introduction

Radiata pine (Pinus radiata D. Don) is economically the most important forest plantation species in New Zealand, Chile and Australia, exceeding 1.5 million ha, 1.27 million ha and 0.7 million ha respectively (Ministry for Primary Industries 2022). Radiata pine breeding in New Zealand began in 1953 through the New Zealand Forest Service (NZFS) using plus-tree material from local landraces, selecting primarily for growth and form traits (Shelbourne and Carson 2019). The Radiata Pine Breeding Company (RPBC) was later established in 2002 as an Australasian radiata pine breeding entity with shareholders in both New Zealand and Australia. There are four harvest-age traits that the RPBC focuses on in its breeding objective: stem volume, corewood density, wood stiffness and branching at rotation age (~ 28 years). These are driven by five selection criteria: diameter-at-breast height, corewood density, branching score, predicted modulus of elasticity and stem straightness measured at 5 to 8 years after establishment.

After material has undergone genetic improvement in the breeding cycle, it is available for transfer through two main deployment pathways: (1) commercial seed orchards which then deploy seed progeny to propagation nurseries and subsequently the production forest, (2) clonal varietal deployment programmes. Radiata pine is a monoecious, predominantly outcrossing species that can distribute pollen over large areas. Open wind-pollinated (OP) orchards are planted in sheltered areas and rely on natural pollination. These orchards are simpler to manage but provide opportunity for inbreeding and generally produce half-sib progeny of lower genetic merit and greater stand variation on average compared with progeny from control-pollinated orchards. OP orchards consist of ~ 40 parents grafted onto rootstocks and arranged in a variety of layouts to minimise inbreeding (El-Kassaby et al. 2014). Control-pollinated (CP) orchards have designated male parents for pollination and female parents for cone production and have a much higher labour requirement compared to OP orchards. Pollination involves collecting pollen from either a single pollen parent (full-sib progeny) or multiple pollen parents (polymix/half-sib progeny) and placing a bag over the female strobili before introducing the pollen from designated male parents (Carson 1986).

Clonal deployment relies on somatic embryogenesis (SE), whereby selected totipotent cell-lines can be rapidly multiplied in ideal laboratory conditions to give rise to an embryo and a subsequent plantlet (Loyola-Vargas and Ochoa-Alejo 2016). This can be used in combination with cryopreservation which acts as a clonal archive for storing embryonic material. To shorten the clonal deployment timeline, fresh green cones are harvested from the top parental crosses and SE is performed to convert somatic embryos to plantlets. Over the past decade, the use of clonal varietal deployment has been steadily rising in New Zealand; however, it has yet to reach its full potential, with only 20% of deployed material coming from clonal deployment methods, and the remaining 80% coming from seed orchards (Ministry of Primary Industries 2021). The uptake of SE derived clones is limited by the specialised lab skills required to develop SE lines, a limited window of opportunity for embryo retrieval and the recalcitrance of a large proportion of genotypes (Montalbán and Moncaleán 2019). To prevent stem breakage, which occurs when juvenile clone stems topple or break in wind-prone sites with high rainfall (Watson and Tombleson 2002), SE radiata clones must go through a 4-year vegetative cutting stage. This ensures clones have developed robust root biomass to make them resilient to adverse weather conditions. The implementation of bioreactors to scale-up SE production could reduce the deployment time to around 2 years (Bonga et al. 2018), but does not resolve the stem breakage problem. Clonal material often undergoes 6-year field trials to test for growth defects before deployment to the production forest. The main benefit of clonal deployment is the ability to rapidly deploy elite material with high genetic merit and create uniformity in the production forest. Any non-additive genetic component can also be better captured which is not possible in seed production and deployment due to recombination. While there is a current bottleneck in producing SE plants, new technologies are showing promise in being able to scale-up SE-derived clonal production (Carson 2019).

The genetic merit of seed orchards is predominantly affected by the connectivity between the orchard and the breeding program through seed orchard turnover (Funda and El-Kassaby 2012). As the breeding program advances in genetic merit, elite material is selected and transferred to the seed orchard population. Material of lower genetic merit or with lower seed yields is removed (rogueing) and replaced with elite material selected from the breeding program in this process. This could be classified as a two-step selection process: selecting older, lower genetic merit material for rogueing increases the average genetic merit of the orchard and reduces seed yield, while selecting new, elite material from the breeding program for replacement also raises the average genetic merit of the orchard but lowers seed yield until the material matures (Funda and El-Kassaby 2012). Increasing the turnover rate of the seed orchard should hypothetically enhance the rate of genetic gain in the seed orchard, assuming that the breeding program is improving in genetic merit. However, the genetic gain obtained by the seed orchard will be dependent on the accuracy of breeding values in both the material rogued from the seed orchard and the new material deployed from the breeding program. Genetic gain is also negatively correlated with genetic diversity and seed yield. Balancing these three factors in seed orchards has been the subject of numerous studies examining different selection methodologies (Lindgren et al. 2004), management strategies (Funda et al. 2009; Kang et al. 2001), and seed orchard designs (El-Kassaby et al. 2014). However, the impact of turnover rate on these three factors has not been widely studied, particularly in radiata pine.

Genomic selection (GS) is a strategy that has gained use in numerous animal (Meuwissen et al. 2001), crop (Crossa et al. 2017), and tree breeding programs (Lebedev et al. 2020). GS relies on the estimation of the relationship between single nucleotide polymorphism (SNP) markers and phenotypic values to predict breeding values in target populations based on genotype alone. To estimate this relationship, a training population must be formed to train a genomic prediction model. Due to changes in linkage disequilibrium (LD) caused by selection and the need to prevent reductions in relatedness (Habier et al. 2007), it is necessary to continuously retrain this model each generation. In traditional forest tree breeding, extensive field testing resulting in generation intervals of over 15 years is the primary factor hindering progress in genetic gain (Grattapaglia and Resende 2011). The main advantage of GS is the ability to select from large numbers of individuals at the seedling stage, i.e. 2 to 3 years after pollination in radiata pine, leading to accelerated breeding cycles and greater genetic gain compared to traditional selection (Bernardo and Yu 2007). The time it takes for genetic gain to be translated to the production forest can be much longer than this if the entry of new improved germplasm to seed orchards or clonal programmes is slow.

Presently, seed orchards and clonal deployment programs select material based on breeding values from best linear unbiased prediction (BLUP) models using phenotype data collected from extensive field trials and information captured by the pedigree. The introduction of genomics creates a new opportunity for material to be transferred to the seed orchard at an earlier stage of the breeding cycle. As the GS-led breeding program advances, material will be available each year from both the GS seedling population and the combined phenotypic-genomic training population (Dungey et al. 2009). However, it is unknown whether the early selection of material for deployment at the seedling stage with a lower breeding value accuracy will result in more genetic gain than deploying material approximately 5 years later with a higher breeding value accuracy when trees have been phenotyped in test trials.

At present, the accuracy of genomic-estimated breeding values is around 0.6 for most key selection criteria (Klápště et al. 2020). However, when phenotypic/pedigree-based estimated breeding values are used, the accuracy could increase to over 0.8 (Radiata Plant Breeding Company 2022). The increase in accuracy levels could be overestimated since pedigree-based relationship matrices do not consider Mendelian sampling, which is taken into account by genomic relationship matrices (Beaulieu et al. 2022). If breeding programs are increasing the rate of genetic gain using genomic selection, then there is a requirement to transfer this gain more effectively through seed and tissue culture deployment pathways if the benefit is to be fully realised. Furthermore, if marker-based selected material is available at a faster rate, then it is important for germplasm producers to understand how census size and turnover rate in a seed orchard can be optimised to cope with the flow of new germplasm and changes in breeding value accuracy and seed yield.

This study focuses on three deployment options: open-pollinated orchards (OP), control-pollinated orchards (CP) and clonal deployment through somatic embryogenesis. This study aims to utilise stochastic simulation to evaluate and optimise each deployment option under two different selection strategies: genomic selection and phenotypic selection (PS). The first objective of this study is to determine how the genetic gain of deployment option differs under genomic and phenotypic selection strategies when accounting for the different deployment timelines. The second objective is to establish an optimal turnover rate and census size for a seed orchard when considering seed yield, genetic merit, and genetic diversity.

Materials and methods

Timeline assumptions and overview of the simulation

This simulation assumes that breeding activities are occurring every year with overlapping generations and improved breeding material is available for entry into the deployment pathways every year. Four deployment pathways are considered with their respective timelines: (1) phenotypic selection through seed orchards (22 years), (2) genomic selection through seed orchards (17 years), (3) phenotypic selection through clonal deployment with field trials (16 years), (4) genomic selection through clonal deployment without field trials (5 years) (Fig. 1). The timeline for PS with clonal deployment could theoretically be shorter (~ 12 years); however, it takes ~ 16 years for phenotypic breeding values to become available and clones selected for deployment must be multiplied after the clonal field-testing phase which takes ~ 4 years. Archiving genotypes in each deployment pathway both acts as an insurance measure to prevent material loss and provides a source of vegetative material for deployment. Archiving is done through either a clonal archive (seed orchard deployment) and/or cryopreservation (clonal deployment). Seed orchards currently use scions grafted onto rootstocks as parents instead of directly planting SE clones; this is because grafting has a higher flowering rate compared to SE clones. If this was not a factor, utilising SE clones as seed orchard parents would reduce the GS seed orchard deployment pathway by ~ 8 years. It was assumed that bareroot seedlings were used to deploy seedlings from seed orchards, in contrast, some forestry companies in New Zealand use containerised cuttings. While bareroot seedlings are cheaper and shorten the deployment timeline, containerised cuttings increase the seed orchard deployment cycle by around ~ 2 years but extend the length of the planting season. GS-based clonal deployment with no field testing is inherently a high-risk practice that is not currently employed in New Zealand, as certain clones could have growth defects such as resin bleeding or internal checking.

Fig. 1
figure 1

Four deployment pathways and their processes: (1) seed orchards with phenotypic selection; (2) seed orchards with genomic selection, seed orchards with phenotypic selection; (3) clonal deployment with phenotypic selection and field trials; (4) clonal deployment with genomic selection and no field trials. Note: seed orchard deployment pathways consist of OP and CP options

Full size image

The length of the deployment cycle is important to calculate the improvement in net-present-value per year ($NPV/Y) determined through a selection index. Three key components determine the improvement of $NPV/Y at the deployment year: (i) the selection intensity into the deployment pathway, (ii) the length of the deployment pathway, (iii) the structure of the deployment pathway. The selection intensity in the deployment pathway describes the replacement rate of the orchard or the number of clones selected for clonal deployment. The length of the deployment pathway is altered by the number of steps to achieve deployment of material. The structure of the deployment pathway has effects including cell lines lost due to SE failure, clonal testing, and the level of clone age-dependent gamete contribution in the seed orchard.

The simulation begins with four cycles of phenotypic selection to create a founder population. The founder population is then used as a source of selection for the first 9 years of breeding cycles B1-9 (Fig. 2). To replicate the rolling front breeding strategy, there is a 9-year cycle for each breeding population; therefore, B0 is crossed to create B9 and B1 is crossed to create B10. Deployment begins with selection of genetically superior material based on information from genomic markers at B0 to form Di where i is the length of the deployment pathway. Twenty cycles/years of deployment were used in this simulation. To account for the length of the deployment pathways, genetic progress is expressed in a per year basis by subtracting the value of the selection index at the final cycle denominated by the number of years in the deployment cycle plus the number of simulation cycles as follows:

Fig. 2
figure 2

Schematic of overlapping deployment and breeding cycles

Full size image

Equation 1 (Clonal Deployment)

$$\Delta {G}_{SI}/Y=\frac{{{{{GV}_{SI}D}_{i+n}-GV}_{SI}B}_{0}}{{D}_{i+n}}$$

Equation 2 (Seed Orchard Deployment)

$$\Delta {G}_{SI}/Y=\frac{{{{{GV}_{SI}D}_{i+n}-GV}_{SI}D}_{0}}{{D}_{i+n}}$$

where \(\Delta G\) is genetic gain of the selection index per year (also interchangeably expressed as \(\Delta G\) in this paper), \({GV}_{SI}\) is the mean value of the selection index, n is the number of cycles simulated, i is the length of the deployment pathway D and B0 is the first breeding population. The difference between calculating genetic gain for clonal and seed orchard deployment is that the base seed orchard population (D0) in most cases will have a different genetic merit from the base breeding population (B0) and subtracting the genetic merit of the improved orchard (Di +n) from the base orchard (D0) will account for this.

To account for the value of time, the discounted NPV for each deployment pathway was calculated (NZ$). Equation 3 provides the net present value of \(\Delta {G}_{SI}/Y\) for each deployment pathway over given number of deployment cycles or years. A 5% discount rate (d) was assumed in this study which is consistent with industry practice (Manley 2022). Although we did not conduct a sensitivity analysis of the discount rate, it is worth noting that an increase in the discount rate could diminish advantages of deployment methods with shorter timeframes (Chamberland et al. 2020). The NPV is expressed per hectare for revenue only; therefore, it does not account for operating expenses:

Equation 3

$${NPV}_{D}=\sum_{j=0}^{n}\frac{\Delta {G}_{SI}/Y }{(1+d{)}^{i+n}}$$

where \({NPV}_{D}\) is the net present value for deployment pathway D, n is the number of cycles simulated (n = 20), i is the length of the deployment pathway and d is the discount rate (d = 5%).

The founder population contained 500 individuals. Individuals were created based on the genomic parameters from the Markovian Coalescent Simulator which simulates genome evolution with multiple cycles of mutation and drift (Chen et al. 2009). This study used a genetic map that was representative of the Pinus radiata genome, consisting of two chromosomes which has been used in other studies (Denis and Bouvet 2013); each chromosome was 150 cM in length which is representative of the longest linkage group (Li and Dungey 2018), with a mutation rate of 2.9 × 10−9 similar to the closely related loblolly pine (Pinus taeda L.) (Brown et al. 2004), which spanned 1.79 Gb. For each chromosome, 37,500 segregating sites were used based on the number of segregating sites detected in studies (Telfer et al. 2018).

This study focused on a multi-trait evaluation consisting of five selection criteria, diameter-at-breast height (DBH) in mm, stem straightness (STR) on a 1 to 9 categorical scale, branch-cluster frequency (BRH) on a 1 to 9 categorical scale, corewood density (DEN) in kg/m3 and modulus of elasticity (MOE) in GigaPascals (GPa). To calculate genetic gain, we used the total additive genetic value (GV) in AlphaSimR as an approximation of the true breeding value. As performed in other simulation studies (Carson 2019), a subset of 1000 individuals from deployed populations were randomly sampled across all families to account for any differences between deployed clonal and seed orchard progeny sizes. The architecture of the selection criteria was infinitesimal across 1000 QTLs per chromosome that were randomly drawn from a gamma distribution with a shape parameter of 0.5 used in other radiata studies (Li and Dungey 2018). Specified correlations between selection criteria (Table 1) and additive variances were supplied by the RPBC (Table 2).

Table 1 Additive genetic correlation matrix between selection criteria
Full size table
Table 2 Genetic variances and heritabilities used for selection criteria and objective traits
Full size table

It is predicted in radiata pine that selection criteria with lower heritability (growth) tend to have higher levels of dominance variance compared to those with higher heritability (wood quality) (Burdon et al. 1992). However, reliable estimates are not currently available to the RPBC, due to the high variability between trial sites (Apiolaza 2012; Baltunis et al. 2010, 2009). Dominance variance also becomes fairly negligible past the age of seven (Dean et al. 2006). Consequentially, we still included dominance variance in the simulation model due to the benefit it can provide particularly to clonal deployment (Wu 2019). Dominance effects were sampled from a mean of zero and 15% additive genetic variance, which is similar to other studies (Li and Dungey 2018).

Starting with the 500 founder individuals, the simulation began with a burn-in cycle of 100 generations to generate linkage disequilibrium (LD) with 1000 crosses and 2 progeny per-family. Historic generations were then created by performing four selection cycles using 100 parents and 200 crosses. PBLUP and GBLUP accuracies (r) used in this study are in accordance with accuracies determined by the RPBC (Table 2) and were used throughout the study when referring to phenotypic and genomic selection respectively. To reduce computation time, historic generations’ breeding values were generated using pre-specified PBLUP accuracies (r) (Table 2) that were generated by stipulating the corresponding R-value between each selection criteria and the breeding value using the ‘rnorm_pre’ function from the ‘faux’ package (DeBruine 2023). This practice was also used to generate GBLUP values, avoiding the need for a training population. Selection was performed using a multi-trait selection index calculation, using an index I:

Equation 4

$$I={b}_{1}{X}_{1}+{b}_{2}{X}_{2}\dots +{b}_{n}{X}_{n}$$

where bi is the index coefficient and X is the phenotypic estimated breeding value (EBV) or genomic estimated breeding value (GEBV) for a given trait. The vector of index coefficient (b) used was:

Equation 5

$${\varvec{b}}= {{\varvec{G}}}_{{\varvec{s}}{\varvec{s}}}^{-1}{{\varvec{G}}}_{{\varvec{s}}{\varvec{o}}}{\varvec{v}}$$

where v is a vector of the economic weights, Gss is the additive genetic variance–covariance matrix between selection traits and Gso is the additive genetic variance–covariance matrix between selection criteria and objective traits measured at log harvest, calculated using a selection criteria-objective trait correlation matrix (Table 3) (Schneeberger et al. 1992). The objective traits used were log volume (VOL) in m3 ha−1, branching index (BIX which is the average size of four branches per 5.5 m log, using the largest branch in each quadrant in mm) (Inglis and Cleland 1982), corewood density (DEN) in kg/m3 and modulus of elasticity (MOE) in GPa. Pearson’s correlation coefficient (rIH) was used to determine the correlation between the EBVs in the breeding population with the GV of an individual. The vector of economic weights was provided by the Radiata Pine Breeding Company (Table 2) (RPBC).

Table 3 Correlation between selection criteria and objective traits
Full size table

The genomic breeding strategy was implemented by selecting individuals from the founder population using genomic-estimated breeding values (GEBVs). After a complete genomic selection cycle (year 9), material was selected each year from the generation prior, providing overlapping generations for the deployment of genotypes. Selection of the breeding material was conducted in two steps: (i) selecting the top number of families per generation based on mean EBVs or GEBVs and selection index, (ii) balanced within-family selection across top families. The population size of the breeding program was 10,000 and 100 parents were selected within 25 top families and randomly crossed twice to form 200 families for the subsequent generation. This was performed over twenty years/cycles. The numbers of families were only modified as a sensitivity analysis for clonal deployment. One hundred iterations were used for each scenario to reduce the computation time. All simulations were run on a 6-core (12 thread) Intel Xeon processor with 128 GB of RAM and an Intel Math Kernel library to enable multithreading.

Deployment

Seed orchards

AlphaSimR was used to create seed orchard populations. The initial parents were comprised of elite material selected using PBLUP values from the founder population and were assumed to be at full seed-producing maturity (7 years of age). Two GS-based scenarios were tested: (1) both new parents from the breeding program and material to be removed from the seed orchard were selected based on GEBVs only (GS scenario); (2) new material is selected from the breeding program based on GEBVs; however, once clones reach the age of five PBLUP values with higher breeding value accuracies were assigned to inform rogueing decisions (GSPS scenario).

It was also assumed that the ramets of each clone were equally distributed through the orchard; therefore, the proportion of the individuals represented in the genetic value of the orchard is equal to their reproductive capacity. Based on current orchard practices, the base assumption was a turnover rate of 10% per annum and a census size of n = 40 in the open-pollinated orchard, and equally the control-pollinated orchard contained n = 20 males and n = 20 females; these base assumptions were modified in a sensitivity analysis. When new clones are introduced to the orchard population, it is assumed that every ramet of a clone selected for replacement would be rogued. There is a stepwise increase in male and female cone production, dependent on clone age. The age of each clone was tracked in the simulation, and this was linked to yield and production data (strobili [flowers] per ramet, year of cone development, number of seeds per kg, stems per hectare) from a New Zealand radiata pine seed orchard. To estimate yield, we used the following equation from the seed orchard data summed over ramets for each year:

Equation 6

$${Y}_{{X}_{i}}={F}_{\left({X}_{i}-2\right)}/(RS)$$

where Y is the yield per hectare (ha) expressed in kilograms (kg) per ha, F is the flowers per ramet, Xi is the ith year of cone development, R is the fixed number of seeds per kg and S is the fixed number of stems per hectare. Therefore, this yield estimate assumes a 2-year development for cones from flowering and a penalty for cone loss.

Male gametes in CP and OP orchards began contributing gametes to the orchard earlier than female cones in accordance with seed orchard data. It was assumed that female and male coning traits are independent of other traits because not enough information is available for correlations between all traits used in the index. In the control pollinated scenarios, we assumed that male pollination contributions were equal, regardless of clone age past the age of two. In contrast, male pollen contributions from OP orchards were estimated based on the link between the number of strobili per ramet and clone age from empirical seed orchard data. Each year, 2250 crosses were performed using the makeCross2 function in AlphaSimR which incorporates a crossing plan, with 15 progenies per cross, approximately equivalent to a 1-hectare seed orchard block. Parents were randomly sampled into a crossing plan based on their weighted contributions to the seed orchard; selfing in the OP orchard was assumed to be 4% (Kumar et al. 2007). This is in accordance with the following modified formula for estimating the genetic merit of an OP orchard from (Funda et al. 2009; Kang et al. 2001):

Equation 7

$$GV=\sum_{i=1}^{N}{I}_{i}\left(\frac{{(f}_{i}+{m}_{i})}{2}1-M\right)+M\sum_{i=1}^{N}{I}_{i}\left(\frac{{(f}_{i}+{m}_{i})}{2}1-M\right)C$$

where fi and mi are the proportion of gamete contributions from female and male sources respectively for the ith clone. The parameters fi and mi are calculated by the sum of female cones or male strobili for the ith clone divided by the total number of female cones or strobili in the orchard, so 0 ≤ fi ≤ 1, 0 ≤ \({m}_{i}\)  ≤ 1, \({\sum }_{i=1}^{n}{f}_{i}\) =1 and \({\sum }_{i=1}^{n}{m}_{i}\)=1. I is the selection index for the ith clone and N is the number of clones that have a corresponding selection index value in clone i. M is the assumed proportion of gene migration due to pollen contamination and C is the assumed mean I value of the contaminating pollen. Pollen contamination rates, quality, and diversity in New Zealand radiata pine seed orchards are relatively unknown. A meta-study found that pollen contamination rates in conifer orchards averaged around 45%; therefore, in this study we assumed pollen contamination is 50% or which is the current assumption used by the industry (Adams and Burczyk 2000). Therefore, in Eq. 7 M = 0.25 because pollen is half of the total gametic contribution to the orchard. To standardise pollen contamination, we assumed that the pollen had a genetic value that was 10% inferior compared (C = 0.90) to the value of the seed orchard, which is in line with similar studies of different conifers (Kang et al. 2005, 2003). It was assumed that there was no pollen contamination in the CP orchard.

In the context of a CP orchard without pollen contamination, Eq. 7 simplifies to:

Equation 8

$$GV=\sum_{i=1}^{N}{I}_{i}\frac{{(f}_{i}+{m}_{i})}{2}$$

Genetic value of the deployment year was determined by randomly sampling 10,000 progenies from the crosses and calculating the mean selection index value.

The rate of inbreeding (ΔFt) is the estimation of the rate in which homozygosity increases between one population (seed orchard) to the next (seed orchard progeny). ΔFt was estimated using molecular markers between the seed orchard parents and the progeny. SNPs were randomly distributed across the genome with 3000 SNPs per chromosome, equivalent to the 36 K radiata pine chip (Graham et al. 2022). Seed orchard progeny were sampled at a rate of n = 1000, and monomorphic SNPs were filtered out; the ‘popkin’ package was used to calculate the inbreeding coefficient (F) for the parents and the progeny (Ochoa and Storey 2021). The rate of inbreeding was calculated as follows:

Equation 9

$${\mathrm{\Delta F}}_{t}= \frac{{(F}_{t+1}-{F}_{t})}{(1-{F}_{t})}$$

where Ft is the inbreeding coefficient of the parental population and Ft+1 is the inbreeding coefficient of the progeny.

Clonal

Somatic embryogenesis (SE) failure rates in radiata pine can be as high as 80%, depending on the methodology used (Montalbán and Moncaleán 2019). In loblolly pine, the genetic control of SE failure is highly variable among families; however, this has not been confirmed in radiata pine (MacKay et al. 2006). To simulate SE failure, random sampling was used to remove 50% of families, equivalent to somatic embryogenesis failure, unless otherwise mentioned. The number of families selected for the deployment process is limited by the narrow-window in which green cones can be harvested from top parental crosses and the SE resources required to convert somatic embryos into plantlets. To simulate the selection of green cones, the top 300 individuals from 25 families were selected using mid-parent breeding values in both the genomic and phenotypic deployment pathway scenarios. GEBVs are not available until ~ 1 year after green cones have been harvested. Therefore, the top 2 clones from the remaining 150 individuals were selected based on GEBVs in the clonal genomic deployment pathway. In the phenotypic deployment pathway, the top 2 clones were selected based on PBLUPs from the remaining 150 individuals.

Results

Simulation 1: Genomic vs phenotypic selection and deployment

Three selection types for this simulation were considered, genomic selection (GS), phenotypic selection (PS) and genomic selection in combination with phenotypic selection (GSPS) (seed orchards only). GS and PS selection types were selected and rogued (seed orchard only) based purely on GS and PS values respectively, whereas the GSPS selection type was selected based on GS values from the breeding program and clones were enhanced with PBLUP breeding value accuracies once the clone age reached 5 years of age. Overall, when deployment time (Di) is accounted for, clonal deployment achieved significantly higher ΔGSI/Y compared to seed orchard options (

Table 4). Genomic clonal deployment options provide greater genetic gain than the breeding program due to the deployment pathway (5 years) being shorter than the breeding cycle (9 years). Genomic clonal deployment provides 10% greater ΔGSI/Y compared to phenotypic clonal deployment. CP orchards provide only a 4.3–6% improvement in ΔGSI/Y compared to OP orchards. GS in seed orchard deployment gave only a slight increase in ΔGSI/Y compared with phenotypic selection; however, GSPS offered a significant improvement in ΔGSI/Y for both CP and OP deployment types. The numbers of families selected for entry into the clonal SE pathway and the simulated loss of families did not significantly affect the ΔG in both genomic and phenotypic pathways (results not shown). When accounting for the discounted value of time in net present value (NPV), GS outperforms PS in all scenarios.

Table 4 Deployment type and selection type change the genetic gain in deployment pathways
Full size table

Simulation 2: Seed orchard turnover in OP and CP orchards

Overall, despite turnover effects, CP orchards only marginally improve ΔG compared to OP orchards regardless of the turnover rate used; the difference in genetic gain is likely due to the additive effect of contamination rates, selfing rate and contamination penalty imposed on the OP orchard (Fig. 3). The level of yield and ΔG were closely negatively correlated, whereby, increases in lower turnover rates (5 to 10%) increased ΔG and decreased yield more significantly compared to higher turnover rates (e.g. 20 to 25%) (Fig. 3 and Fig. 4). As the contamination rate (M) decreases, the level of ΔG increases, and when there is no contamination (M = 0), the OP orchard and the CP orchard have similar levels of ΔG due to random selection of parents in both orchard types (Fig. 3).

Fig. 3
figure 3

The relationship between genetic gain and seed orchard turnover rate with different deployment types (OP and CP) and contamination rates (M). Note: The contamination penalty (C) is assumed to be 0.9

Full size image
Fig. 4
figure 4

The relationship between yield and turnover rate in an OP orchard

Full size image

Increasing the turnover of males in CP orchards to high levels (25%) while maintaining a low female turnover (5%) returned the same level of ΔG as a turnover rate of 10% for male and females but with 15% greater yields (Fig. 5). Consequentially, increasing the turnover rate of male and female CP parents also dramatically increases inbreeding (ΔFt); however, if a low turnover rate is maintained in the females (5%) and a high turnover rate is maintained with males (25%), then the inbreeding remains relatively stable.

Fig. 5
figure 5

The relationship between inbreeding rates, yield and genetic gain in a CP orchard with different male and female turnover rates

Full size image

Interestingly, withholding from rogueing clones in a GS or PS orchard until they reach full seed yield maturity (age 7) is counterproductive to the total yield of an orchard as it reduces the average clone age of the orchard; therefore, rogueing decisions should be based on breeding value unconstrained by the age of the clone (Fig. 6). In contrast, in a GSPS scenario, the age of rogueing constraint did not affect seed yield or ΔG (results not shown).

Fig. 6
figure 6

The relationship between clone age constraint and yield in an OP orchard. Note: The clone age constraint is the minimum age that a clone in an orchard is eligible for rogueing. Simulation was performed on a GS-only OP orchard with 15% turnover per annum

Full size image

The genetic penalty (C) of the contaminating pollen had a greater effect on ΔG compared to the contamination rate (M) (Fig. 7). Every 0.1 increase in M represents a 20% increase in contamination because pollen is half of the total gametic contribution to the orchard. Therefore, if M = 0.25 then the pollen contamination rate is 50%. The value ΔG decreased by − 2.6 ΔG for every 0.1 increase in M. Every 0.1 increase in C represents a 10% penalty reduction penalty applied to the total genetic value of the seed orchard, relative to the level of pollen contamination (M). If M = 0.25, every 0.1 increase in the unit C reduced ΔG by − 7 ΔG.

Fig. 7
figure 7

Comparing the relationship between contamination penalty (C), contamination rate (M) and genetic gain in an OP orchard. Note: The dashed line indicates the base assumption used in this study (M = 0.25, C = 0.9)

Full size image

Increasing the census size of seed orchards only slightly reduces ΔG, although this decrease is more substantial between n = 20 and n = 40; it is estimated that for every n = 20 increase in census size, ΔG is reduced by only − 1 ΔG (Fig. 8A). Conversely, the inbreeding rate between the seed orchard parents and the progeny is high altogether in this study and is closely correlated with ΔG, whereby increasing census size decreases the inbreeding rate (Fig. 8B). Increasing the turnover rate of open pollinated orchards also increases the inbreeding rate; however, there is a rapid increase in ΔFt between lower turnover rates (5–10%) which eventually plateaus in higher turnover rates (20–25%) (Fig. 9).

Fig. 8
figure 8

The relationship between census size, genetic gain (plot A) and inbreeding rates (plot B). Note: Simulation performed in OP orchard with a turnover rate of 10%

Full size image
Fig. 9
figure 9

The relationship between inbreeding rate and seed orchard turnover rate. Note: Simulation performed in OP Orchard with a census size of n = 40

Full size image

Discussion

This study shows that selecting material for deployment with genomic selection provides the same benefits as phenotypic selection. However, when genomic selected orchard clones above the age of five are enhanced with phenotypic data (GSPS) obtained from training/validation test trials, it makes genomic deployment significantly more beneficial compared to phenotypic deployment ΔG. This is because seed orchards undergo a two-step selection process when turning over an orchard, firstly to select material from the breeding program and secondly to rogue material from the seed orchard. The addition of phenotypic-enhanced breeding values at a later stage will improve the confidence of seed orchardists to make early selection decisions with genomic information, knowing that selection accuracies will be improved once phenotypic information is released. Furthermore, this method requires the breeding program to ensure that material deployed in seed orchards is included in test trials. Interestingly, this study found in the GSPS scenario, having a clone age constraint for rogueing made no difference to yield or ΔG. In contrast, when using either GS or PS alone, removing the clone age constraint significantly increased yield and did not impact ΔG. Theoretically, removing only older parents could reduce the overall clone age of the orchard and, therefore, yield. However, this theory is limited by the fact that this study assumed that clones produced the same yields dependent on clone age. Including flowering data as a trait could be included in future simulations if there is more data available on its variance, heritability, and covariance with other selection criteria.

The length of the deployment cycle is one of the main contributors to the genetic gain per-year of a deployment program. While using a genomic selection strategy is a potentially feasible way to reduce the length of deployment timelines, other areas of a research could reduce these timelines further. One of the fastest ways to reduce the deployment timeline in a genomic seed orchard scenario is using SE plantlets to establish new seed parents instead of scions grafted onto rootstocks. This could reduce the genomic seed orchard deployment timeline from ~ 17 years to ~ 9 years. However, SE plantlets are not typically used as seed orchard parents in New Zealand. This is because SE plantlets have not been available to seed orchard producers and anecdotally have lower flowering rates and seed yields compared with parents grafted onto rootstocks. To our knowledge, there have been no published studies in radiata pine in this area. Research to improve flowering in non-grafted clones is required to reduce the deployment lag further. It is notable that SE plantlets would only be beneficial as parents in a genomic seed orchard deployment scenario. This is because a phenotypic-based seed orchard scenario would have to wait ~ 15 years for breeding values to become available before making selections, which is longer than the time it takes for scions to become available (~ 10 years). To counteract the reduction in flowering from SE plantlets, a greater number of pollen parent ramets could be used in a CP orchard to enhance genetic gain and counteract the loss of flowering.

Increasing the number of genotypes (census size) in a genomic-selected orchard has negligible effects on ΔG but larger effects on inbreeding. In the context of genomic selection with low to moderate accuracy, it could be thought that increasing the population size could spread the risk of bringing new material into the orchard. However, when assuming that breeding values are normally distributed, the variance in genetic value between elite individuals across all traits appears to be low, which could have contributed to the only slight decrease between small (n = 20) and large (n = 100) census size populations. Increasing the census size from a breeding program with a larger effective population size could decrease the inbreeding between the seed orchard parents and progeny further due to the increase gene diversity contained within the breeding program. However, this was not tested in this study. Future studies could include other factors such as family selection, fecundity variation and spatial dynamics that would provide a more accurate estimate of inbreeding rates (Gömöry et al. 2007). It is likely that radiata pine seed orchards will opt to use small census sizes because genetic diversity is less of a concern in radiata due to the short rotation times (Burdon and Aimers-Halliday 2003). Moreover, radiata pine is exotic to New Zealand so there is no legislation or requirement to maintain a certain level of diversity in the plantation forests, as there is in some conifers in certain countries. Smaller census sizes could also be used in conjunction with optimised contribution selection in the case of CP orchards (Sonesson and Meuwissen 2000) or optimum neighbourhood seed orchard design with OP orchards (Chaloupková et al. 2016; Funda et al. 2009) and future simulations could investigate these practices. Therefore, while increasing the census size to increase the genetic diversity could potentially be done without significantly hampering the genetic value of the population, there is arguably little benefit unless genetic diversity and inbreeding depression is of concern.

Utilising a selection index in the context of a high turnover seed orchard or clonal deployment program could prove to be imperative in maintaining the connection between selection in the breeding program and advancement of the realised breeding objective after deployment. Strategies to utilise selection indices have been reported internationally in Eucalyptus dunnii (Zelener et al. 2005) and Eucalyptus nitens (Johnson 1996) seed orchards. South American Eucalyptus spp. breeding/silvicultural programs utilise selection indices in both breeding and deployment (Rezende et al. 2014). However, in New Zealand, seed orchards often provide germplasm based on customer preferences and therefore deployment choices are not necessarily based on objective economic decisions. In a deployment program with a slow turnover rate this can remain the status quo over long time frames, suggesting that deployment gain, in turn, will be slow. With the implementation of faster deployment and potentially higher seed orchard turnover through genomics, the lag between tree selection decisions and log harvest could be shortened considerably. As used in this study, employment of selection indices in both the breeding program and deployment would maintain a sense of vertical integration, where economic information will help influence deployment selection decisions.

The advantages of clonal deployment depend primarily on timing and the genetic components of the species, including non-additive effects which were modelled in this study. GS clonal deployment provides a significant advantage compared to PS clonal deployment in terms of ΔG. However, GS clonal deployment as described in the current study is unlikely to be implemented operationally as it carries too much possible risk. There is a requirement to check for growth and other defects in clonal trials, where a limited number of clonal varieties are being deployed at scale. Internal checking and resin bleeding are an example of growth defects found in clonal trials with genomic-estimated breeding value accuracies of rGEBV = 0.55 and rGEBV = 0.56 respectively (Li et al. 2019). This makes it unlikely to be able to screen for these traits at an early age. However, discovery of QTLs with significant effects, advancement of early-screening for growth defects and increasing the number of clones deployed could make GS clonal deployment plausible. But this approach would not account for the unmasking of deleterious recessives that has been seen in bull breeding (Pryce et al. 2014). Li and Dungey (2018) found that clonal deployment using genomic selection had gains of 3–12% over deploying progeny which was less than the estimate from this study. This is likely because their simulation relied on the genetic value of the individuals and did not consider rogueing and rotation of individuals in a mature seed orchard. The seed orchards also contain individuals from historic generations to maintain adequate seed yields which inevitably lowers the overall genetic value. Interestingly, this study found that the loss of individuals due to SE failure did not impact ΔG; as mentioned, this is probably because the variance in genetic value between elite individuals across all traits appears to be low, negating the loss of families. The interaction between the number of families selected and the rate of family loss due to SE failure was also non-significant. This means that the breeding program could reduce the number of families selected for SE as a cost-saving measure without losing ΔG in the deployment pathway. It is not known if the SE recalcitrance observed in operational laboratories has a heritable component.

When comparing the CP and OP orchards, random mating was assumed for both orchard types, and while dominance effects were included, specific combining ability (SCA) was not considered. As a result, the differences in genetic gain between the orchard types could only be attributed to contamination rates and penalties, as well as the contribution of male pollen from older clones in the OP orchard. It was assumed that the contribution of male gametes in the CP orchard was equal regardless of clone age past the age of two, while in the OP orchard, male gametes from older clones had a greater gametic contribution. CP orchards can exploit dominance effects by using progeny data from crosses to incorporate SCA in the selection of parental combinations (Li et al. 2021). In contrast, dominance effects are much easier to exploit in clonal deployment due to the direct capture of non-additive variance and the total additive variance as clones do not undergo recombination (Wu 2019). Further studies could explore the impact of dominance variance on genetic gain in both OP and CP orchards by modelling varying levels of dominance variance in the trait architectures and implementing selection models that account for SCA.

Contamination control in OP seed orchards is sometimes practised, primarily through the planting of non-orchard species in borders or situating OP orchards in a location that is distant from production forests (Codesido and Fernández-López 2014). In a Douglas fir (Pseudotsuga menziesii) seed orchard, contamination reduced benefits by up to CAD$748,500 per year at the highest site index with a planting density of 1111 trees/ha (Stoehr and Webber 2018). Modern techniques such as isolation tents or covers could be employed which have been shown to reduce contamination rates by 23% (Torimaru et al. 2013). While increasing contamination rates is likely to decrease genetic value, it could also inadvertently increase the effective population size of the seedlots, depending on several factors, including the genetic variability of the pollen donors and the genetic distance between the pollen donor and orchard parents (Sønstebø et al. 2018). Gaining an accurate estimate of effective population size including pollen contamination rates in radiata pine would require a study using genetic markers. There is current research underway in North Canterbury, New Zealand, on OP radiata seed orchards to estimate the value and quantity of contaminating pollen using pollen traps and genetic markers. While contamination rates are easier to control compared to the value of contaminating pollen, seed orchards could subsidise and encourage the planting of stock with a higher genetic value in surrounding forests.

This study found that increasing the turnover in an orchard positively affected the genetic value of the orchard, but significantly reduced the yields. To our knowledge, it is not common practice in New Zealand CP seed orchards to turnover male pollen parents at greater rates compared to female parents. Consequently, we examined whether increasing the turnover of males in the CP orchard would be of benefit. Based on the data we have presented, a CP seed orchard with a low census number (20 males, 20 females), low female turnover (5% per annum) and a high male turnover (25% per annum) would generate high yields with high ΔG while moderating rates of inbreeding ΔFt = 0.026). Increasing the turnover in OP orchards resulted in a significant increase in ΔG at the cost of losing yield and increasing inbreeding rates. The increase of inbreeding at higher turnover rates could be caused by the reduction of gamete contributions due to the greater number of younger parents in the orchard producing less pollen and cones. To counteract this, higher turnover rates could be utilised with higher census sizes. Due to a lack of empirical data, this study did not address the effect of turnover on pollen yields which is a consideration in CP seed orchards. Higher rates of male turnover in a CP orchard would likely require greater land space to accommodate the extra number of male parents established in a pollen archive to satisfy pollen demands.

While this study analyses the relative benefits of each deployment pathway, it does not account for their respective costs. Cost–benefit analyses have been performed with GS in combination with seed orchard and SE-based clonal deployment in white spruce (Chamberland et al. 2020). Where despite costs, GS based seed orchard/SE clonal deployment pathways were found to result in higher land expectation values compared to conventional seed orchard deployment pathways (Chamberland et al. 2020). However, in a combined white spruce and lodgepole pine (Pinus contorta), genomic selection with OP seed orchard deployment would only be justifiable compared to conventional selection at lower discount rates, decreased seedling costs and higher plantation areas (Chang et al. 2019). Future studies will utilise data generated by this study and incorporate costs to investigate the cost–benefit comparison of OP and CP orchards, PS and GS-based deployment and seed orchard turnover rates. A cost benefit-analysis could also identify which factors are the greatest drivers that control the profitability and feasibility of deployed material to the production forest. Because the difference in genetic gain between CP orchards and OP orchards is marginal (Fig. 3), it is uncertain whether CP orchards would yield any benefit when higher labour and management costs are accounted for. Furthermore, there are also added complexities and technological investment to ensure that CP pollination is carried out efficiently. Increasing turnover rates evidently decreases seed yield; therefore, greater land space and higher management costs would be required to meet the demand of industry for new planting stock. Analysing the cost–benefit of each turnover rate could determine the optimum rate for OP and/or CP seed orchards.

Conclusions

This study assumed that phenotypic selection had higher breeding value accuracy compared to genomic selection; however, phenotypic selection-based deployment had a longer timeframe compared to genomic selection. When accounting for the clonal and seed orchard deployment pathway timeframes in a stochastic simulation, genomic selection has the propensity to improve the delivery of genetic gain between the breeding program and production forest in radiata pine. In clonal deployment, material deployed using genomic selection significantly enhances genetic gain per year due to the large reduction in deployment timeframes (16 years for phenotypic selection, 5 years for genomic selection). Yet, in seed orchards, the difference in timeframes between selection strategies is less (24 years for phenotypic selection, 19 years for genomic selection) because the breeding cycle needs to be completed (~ 9 years) to harvest scions. This resulted in similar genetic gain per year between genomic and phenotypic-based deployment selection strategies for seed orchards. However, when the value of time is accounted for through discounted NPV, genomic selection provided significantly greater value compared to phenotypic selection. Additionally, when genomic selected orchard clones above the age of five are enhanced with phenotypic breeding values (GSPS), genetic gain is greater compared to phenotypic deployment alone.

Increasing the turnover rate of seed orchards has a large effect on seed yield, inbreeding rates, and genetic gain. The average clone age of the seed orchard is directly linked to the yield of the orchard. Increasing the average clone age of the seed orchard also increases the number of mature seed producing parents, and this can be achieved by reducing turnover rates. However, while decreasing turnover rates in CP and OP orchards increases yield and reduces inbreeding rates, it also inversely decreases genetic gain. This is because connectivity is reduced between the genetic gain obtained in the breeding program and the deployment of improved material to the seed orchard. CP orchards provide the opportunity to balance yield, genetic gain, and inbreeding rates. This can be achieved by lowering female-parent turnover rates and increasing male pollen-parent turnover rates, resulting in high yield, high genetic gain, and moderated inbreeding rates. This is because the high turnover of male-pollen parents maintains connectivity with the breeding program while a lower turnover of female parents maintains high yields and moderates inbreeding. Increasing census size causes a slight reduction in genetic gain, but significantly reduces inbreeding rates; therefore, the census size of the seed orchard should be dependent on whether genetic diversity or genetic gain is of greater priority.