Leveraging UAV spectral and thermal traits for the genetic improvement of resistance to Dothistroma needle blight in Pinus radiata D.Don

Introduction: Phenotyping is critical in tree breeding, but traditional methods are often labour-intensive and not easily scalable. Resistance to biotic and abiotic stress is a key focus in tree breeding programmes. While heritable traits derived from spectral remote sensing have been identified in trees, their application to tree phenotyping remains unexplored. This study investigates in-situ high-throughput hyperspectral and thermal imaging for assessing Dothistroma needle blight (DNB) resistance in Pinus radiata D.Don.

Methods: Using UAV-based hyperspectral and thermal imaging during a severe DNB outbreak in a clonal trial in New Zealand, we computed narrow-band hyperspectral indices (NBHIs), canopy temperature indices, radiative transfer inverted plant traits, and solar-induced fluorescence. Visual severity scores and remote sensing indices were modelled using spatially explicit mixed-effect linear models integrating pedigree and genomic data in a single-step genomic evaluation. Multi-trait models and sampling simulations were used to evaluate the potential of remote sensing indices to supplement or replace traditional phenotyping.

Results: Remote sensing indices exhibited narrow-sense heritability values comparable to severity scores (up to 0.37) and high absolute correlation coefficients with severity scores (up to 0.79). Carotenoid and chlorophyll-related NBHIs were the most informative, reflecting physiological impacts of DNB. Combining partial visual scoring with NBHIs maintained high estimated breeding value (EBV) accuracy (0.68) at 50% scoring and moderate accuracy (0.59) at 20% scoring. EBV correlation with full scoring was above 0.8 even at 20% scoring. Using solely the most heritable NBHI achieved 0.71 breeding value accuracy and 0.79 absolute EBV correlation with severity scores, suggesting NBHIs can replace visual scoring with minimal precision loss.

Discussion: By utilising UAV-based hyperspectral and thermal imaging to capture single-tree phenotypes related to disease in a forestry trial and pairing the data to genomic evaluation, this study establishes that remote sensing data offers an efficient, scalable alternative to traditional phenotyping. Our approach constitutes a major step towards characterising specific physiological responses, facilitating the discovery of the genetic architecture of physiological traits, and significantly enhancing genetic improvement.

Introduction

Phenotyping refers to the process of observing and measuring an organism’s physical and biochemical traits. It is a fundamental step in the process of breeding to genetically improve species of commercial importance (Isik et al., 2017). It is however a time-consuming process, especially in large and long-generation organisms such as trees. Traditional tree phenotyping mostly relies on human assessment by foot in remote locations, and cannot easily be scaled up to larger trials, therefore limiting the amount of genetic material that can be tested. The precision of measurements is also limited for two main reasons. First, tree breeding trials are established in areas with limited control over environmental variables, creating extra variability in observable phenotypes. Secondly, many quantitative traits are measured subjectively using categorical scores with coarse scales, limiting the consistency of measurements within and across breeding trials, especially when the scoring is performed by multiple people (Poland and Nelson, 2011). The impacts of phenotyping limitations on selection efficacy and speed in tree breeding for forestry purposes have been referred to as the “phenotyping bottleneck” (Dungey et al., 2018). Remote sensing can alleviate many of these limitations by replacing or supplementing traditional measures of tree phenotypes, providing faster and objective measurements to ultimately improve genetic gains. High-throughput phenotyping has been in operation for more than a decade in common agricultural crop breeding systems (Song et al., 2021), with successful recent development of hyperspectral and thermal phenotyping using unmanned aerial vehicles (UAVs) in several crops (Liu et al., 2025; Zhang et al., 2025). However, its application in tree breeding still faces considerable challenges. Our research addresses the use of UAV hyperspectral and thermal measurements in breeding for disease resistance traits in a conifer species widely used in forestry, Pinus radiata D.Don (radiata pine).

Radiata pine is one of the most commercially important tree species in the southern hemisphere and the most planted forestry species in New Zealand. A breeding programme was established in New Zealand in the 1950s and remains very active. As with many pine species, radiata pine is prone to infection and damage by foliar diseases. Dothistroma needle blight (DNB), which is primarily caused by the pathogen Dothistroma septosporum (Dorog.) M. Morelet (Barnes et al., 2004), is characterised by the appearance of brick-red bands (1–3 mm wide) on needles, expanding to tissue dessication and eventually needle loss. The substantial productivity losses associated with the disease in Australia and New Zealand have prompted the inclusion of DNB resistance in radiata pine breeding programmes (Carson, 1989; Dungey et al., 2018). Screening for DNB symptoms is routinely performed in nurseries and breeding trials, and as DNB primarily affects juvenile stands (Bulman et al., 2013), the disease can occasionally occur in genetics trials before or at the usual selection age (6–9 years post establishment). Field phenotyping for resistance to DNB is undertaken using severity scores that are determined by a visual assessment of the percentage of the tree crown that is affected, in 5% increments (Bulman et al., 2004). Breeding for DNB resistance would therefore likely benefit from a more precise, objective and scalable phenotyping approach. Recent research on infected breeding trials reported narrow-sense heritability values for visual severity scores mostly between 0.29 and 0.35 (Li et al., 2018; Ismael et al., 2020; Klápště et al., 2020) but up to 0.43 (Klápště et al., 2020).

Remote sensing research has shown promising results in phenotyping individual tree health within forests and tree crops by detecting changes in physiological functioning caused by biotic stressors (Hernández-Clemente et al., 2019; Smigaj et al., 2023). Paired with UAV technology, remote sensing offers an efficient way to repeatedly collect high spectral and spatial resolution data (Smigaj et al., 2019; Oerke, 2020; Pan et al., 2023). Spores from Dothistroma ssp. infect pine needles through the stomata and create chlorotic lesions that expand into visible necrotic bands caused by the polyketide toxin, dothistromin (Bassett et al., 1970). Dothistromin targets chloroplasts and degrades photosynthetic pigments (Bradshaw, 2004; Kabir et al., 2015). Vegetation indices from spectral imagery can therefore detect DNB infection and characterise severity, as demonstrated in lodgepole pine (Smigaj et al., 2019) and radiata pine (Watt et al., 2023a, b). Water loss is commonly associated with foliar diseases, and a common physiological response is the closure of stomata to reduce transpiration. This leads to a reduction in stomatal conductance and assimilation rate, which can be detected with vegetation indices (Hernández-Clemente et al., 2019), and to an increase in leaf temperature, which can be detected using thermal sensing at the canopy level, as demonstrated in Scots pine (Smigaj et al., 2019). The use of UAV hyperspectral data has been extended to the combination of vegetation indices with plant functional traits inferred from radiative transfer modelling, improving disease severity prediction from the use of vegetation indices alone (Watt et al., 2023b).

The potential for spectral data and vegetation indices to be used as phenotypes in breeding has been explored in a number of tree species. Virlet et al. (2015) used multispectral imagery to identify quantitative trait loci associated with yield in an apple tree hybrid population and found broad-sense heritabilities between 0.60 and 0.77 for three multispectral indices. Tao et al. (2021) assessed the heritability of several UAV-derived vegetation indices to use in tree growth assessments in slash pine (Pinus elliottii) and found narrow-sense heritabilities ranging from 0.07 to 0.26. Čepl et al. (2018) conducted a similar study on Scots pine and found a similar range of heritabilities for vegetation indices (0.02 to 0.22) and up to 0.39 for single wavelength bands at the red edge inflection point, without however relating these metrics to any measured physiological trait. Corbin et al. (2024) used a handheld spectrometer to measure vegetation indices in common gardens of Populus fremontii and calculated broad-sense heritabilities up to 0.48. These exploratory studies showed that some spectral measurements can be heritable enough to be used in breeding. However, further research is needed to relate those measurements to specific adaptive or desirable traits in commercially planted tree species and to implement them within breeding trials in a large-scale phenotyping setup involving UAVs. To our knowledge, there has been no such attempt to date and no genetic assessment of spectral or thermal traits for selection of specific characteristics in tree breeding programmes. Here we address this key research gap by calculating tree-level UAV-collected hyperspectral and thermal measurements in a radiata pine breeding trial experiencing a disease outbreak and using derived indices to directly assess disease susceptibility for genetic selection.

Using UAV collected hyperspectral and thermal data, this study assesses how high-throughput remote sensing data can replace or complement visual severity scoring for DNB resistance. Remote sensing data was collected from a clonal breeding trial located in the central North Island of New Zealand at the peak of a DNB outbreak and linked with severity scores of the disease. We used narrow-band hyperspectral indices (NBHIs), thermal indices and plant traits calculated at the tree level as a means to supplement or replace time consuming measurements of disease severity. We first identify the most promising indices for use in phenotyping DNB severity for breeding. We then focus on two approaches for using the selected indices. The first approach assesses whether high-heritability indices can complement visual scoring, either by enhancing the accuracy of estimated breeding values (EBVs) for severity scores measured on all trees, or by maintaining high accuracy of breeding values for severity scores measured on a reduced proportion of the trial. The second approach determines whether using only remote sensing indices is a viable option to phenotype breeding trials for DNB resistance.

Materials and methods

Genetic material

The study focuses on a clonally replicated breeding trial established by the New Zealand Radiata Pine Breeding Company (RPBC). The trial is in Kinleith forest in the Waikato region of New Zealand (Figure 1A) and comprises 60 6×6 row-column incomplete blocks of 36 single tree plots, with genotype allocation following the “Optimal Design” described by Butler (2013), which takes into account potential spatial correlations in rows and columns within blocks. In total, 1886 genotypes from 163 full-sib families (181 parents) are represented in the trial, with 1–3 ramets each. A pedigree with information for up to 5 generations was available for all tested genotypes. Genomic data was available for 65% of the tested genotypes and their parents, consisting in a set of ∼ 9.5k single nucleotide polymorphisms (SNPs) selected from two arrays: one from exome capture (Telfer et al., 2019) and the second one from a custom radiata pine Affymetrix Axiom array, NZPRAD02 (Graham et al., 2022). Markers in NZPRAD02 have a reasonably even distribution across linkage groups (Graham et al., 2022) and a set of markers almost identical to this used in this study has successfully been used for pedigree reconstruction in radiata pine (Klápště et al., 2022). Genotypes with <0.2 call rate and SNPs with <0.1 call rate were discarded from the analysis. The final genomic dataset had 1453 genotypes and 9570 markers.

Figure 1

Figure 1. (A) Location of the trial site in the Central North Island, Rotorua, New Zealand. (B) Visualisation of the captured hyperspectral data (within the red box) overlaid with thermal data.

Disease assessments

A DNB outbreak was observed in the trial at year 5 after trial establishment. At the peak of the outbreak, severity scores were obtained on 25 September 2023 for all live trees in the trial (n= 2054). The severity was determined as the percentage of the crown affected by the disease, in 5% steps (Bulman et al., 2004). As the disease moves upwards from the base of the crown, this scoring system provides an objective means of characterising disease severity.

Remote sensing captures

Remote sensing data captures were undertaken from a UAV close to the time of visual scoring (Figures 1B, 2). Hyperspectral data were collected on 20 October 2023 using a Nano series HeadWall Photonics sensor (HeadWall Photonics, Inc., Bolton, MA, USA) mounted on a DJI Matrice 600 Pro platform (DJI, Shenzhen, China). The sensor featured a 17 mm focal length, with a 15° angular field of view (FOV). The collected hyperspectral imagers consisted of 270 spectral bands across the visible to near-infrared (VNIR) range from 400–1000 nm with a 6 nm full width at half maximum (FWHM). The UAV was flown at mid-day at 6 m/s and 120 m altitude with 40% side overlap, resulting in a 5 cm/pixel ground sampling distance (GSD). Solar irradiance measurements across a range of 350–2500 nm were recorded concurrently every 15 seconds using a Spectral Evolution RS-5400 spectroradiometer (Spectral Evolution, Haverhill, MA, USA). The raw digital numbers from the hyperspectral data were processed to orthorectified reflectance cubes using SpectralViewTM software (Headwall Photonics, Fitchburg, MA, USA). Radiance was calculated with dark current corrections and calibration files, then converted to reflectance using a 3 × 3 m² ground calibration target located within the flightpath. The reflectance data cubes underwent an orthorectification process to accurately geolocate the imagery and correct any roll, pitch, and yaw imagery distortions. Using post-processed RTK files, previously collected high-resolution digital terrain models (DTM), RGB ortho-mosaics, and the GCPs, the reflectance cubes were orthorectified, mosaiced, and registered to the Universal Transverse Mercator (UTM), Zone 60 South, projection.

Figure 2

Figure 2. The UAV platforms used for data acquisition, namely the DJI Matrice 300 with FLIR Vue^® TZ20-R thermal sensor (top left) and DJI Matrice 600 with HeadWall Nano hyperspectral sensor (bottom left). Also shown are examples of trees exhibiting Dothistroma symptoms (top right) and a UAV pilot in the process of capturing data over the site (bottom right).

Thermal data were collected on 13 October 2023 using a FLIR Vue^® TZ20-R sensor (Teledyne FLIR, USA) mounted on a DJI Matrice 300 RTK platform (DJI, Shenzhen, China). The sensor has a 12 μm pixel pitch, operates within the 7.5–13.5 μm spectral range, and has a reported thermal sensitivity of 85 mK (at F/1.0). To support flight parameterisation and thermal index calculations, ambient air temperature and relative humidity were continuously recorded using an HMP155A probe (Vaisala, Vantaa, Finland) connected to a CR1000 datalogger (Campbell Scientific, Logan, UT, USA). Prior to data acquisition the thermal camera was powered on and allowed to stabilise for 15 minutes. Ambient temperature, relative humidity, and an emissivity value of 0.98 were entered into the DJI Pilot controller software (DJI, Shenzhen, China). Data were captured using the 4.9 mm lens (95°FOV) at an altitude of 80 m and a flight speed of 3 m/s, with 90% forward and side image overlap. The resulting imagery had an approximate ground resolution of 23 cm, with radiometric temperature values embedded in each pixel and stored in the proprietary R-JPEG format. Images were subsequently processed and orthomosaicked in Pix4D using the FLIR Vue Pro R camera profile.

To enable extraction of pure canopy pixels from thermal and hyperspectral data, LiDAR was collected and processed following Watt et al. (2024b). This facilitated individual tree segmentation and the generation of canopy polygons for guiding image-based data extraction. LiDAR was acquired on 8 September 2023 using the DJI Zenmuse L1 sensor on the Matrice 300 RTK platform (DJI, Shenzhen, China), flown at 55 m altitude and 3 m/s in a grid pattern with 10 m line spacing, achieving 85% overlap. The system operated in repetitive scan mode at 160 Hz, recording up to three returns per pulse. For accurate geolocation, the DJI Matrice 300 used a real-time kinematic (RTK) base station and 11 high-contrast ground control points (GCPs) with matte black and retro-reflective white panels, evenly distributed across the site. GCP centers were identified from the LiDAR point cloud and used to georeference the imagery. Post-capture, standard pre-processing steps were undertaken that included denoising, ground point classification, digital terrain model (DTM) generation and ground normalisation using LASTools (v2.0.2) following (Watt et al., 2024b).

Tree detection and crown delineation were performed in R (R Core Team, 2023). A 25 cm resolution canopy height model (CHM) was generated using the pit-free method of (Khosravipour et al., 2014) and then smoothed with a 3×3 pixel moving window. Tree peaks were identified using a local maxima algorithm with a 3 m window, reflecting initial planting spacing, and verified against high-resolution imagery and field-validated stem maps. The “mcws” function in the ForestTools R package (Plowright, 2025) was used to delineate crowns via the watershed algorithm, restricting delineation to the top 75% of the CHM to avoid low vegetation. The resulting crown polygons were visually checked and erroneous crowns removed.

Extraction of hyperspectral and thermal predictors

The hyperspectral and thermal data extracted for each canopy was used to calculate several plant traits and indices that could highlight tree sensitivity to infection. A range of narrowband hyperspectral indices (NBHI), designed to capture variations in plant biochemical and structural properties, were calculated using the formulations in Supplementary Table S1.

The PRO4SAIL2 radiative transfer model, which couples the leaf model PROSPECT-D (Féret et al., 2017) with the canopy model 4SAIL2 (Verhoef and Bach, 2007), was used to simulate canopy structural traits and leaf biochemical properties from pure vegetation pixels in the radiata pine canopies. By isolating pure canopy pixels, the influence of shade and soil was minimised reducing structural artifacts (Zarco-Tejada et al., 2021). This enabled accurate simulation of chlorophyll a+b (C_a+b), anthocyanin (Anth.), and carotenoid (C_x+c) concentrations as well as structural parameters such as the leaf inclination distribution function (LIDF_a) and leaf area index (LAI). A lookup table with 200,000 simulations was generated in forward mode by running the PRO4SAIL2 model, using the parameters and ranges in Supplementary Table S2 which were randomly varied based on continuous uniform distributions, with other model settings kept at default. Data resampling was undertaken to match the hyperspectral sensor’s spectral resolution using a Gaussian response function. Model inversion was performed using Support Vector Regression (SVR), a non-parametric approach grounded in statistical learning theory (Vapnik, 1999). SVR models were trained in parallel using MATLAB (Statistics and Machine Learning, Parallel Computing, and Deep Learning toolboxes; MathWorks, Natick, MA, USA), with resampled reflectance as input and predicted biochemical and structural traits as outputs. A radial basis function kernel was used, with hyperparameters optimised during training.

Solar-induced fluorescence (SIF), a proxy for photosynthetic efficiency, was estimated using hyperspectral radiance and solar irradiance data through the Fraunhofer Line Discriminator (FLD) approach (Plascyk, 1975). Specifically, the three-band FLD method (3FLD) proposed by Maier et al. (2004) was applied. This method extends the traditional FLD technique by incorporating three spectral bands — two situated outside the O₂-A absorption band and one within it — to improve the accuracy of SIF retrievals.

Thermal imagery data and field-measured ambient temperature were used to compute the normalised canopy temperature index (T_c – T_a), where T_c represents canopy-level temperatures and T_a reflects the ambient air temperature recorded during the flight period (Smigaj et al., 2023).

Univariate analyses of severity scores and remote sensing traits

Single phenotypes from visual assessment (severity scores) and remote sensing captures were analysed using the following linear mixed model implemented in the ASReml-R statistical package (Butler et al., 2009):

$$\mathit{\text{y}}=\mathit{\text{X}}\mathit{b}+{\mathit{\text{Z}}}_{\mathit{\text{a}}}{\mathit{\text{u}}}_{\mathit{\text{a}}}+{\mathit{\text{Z}}}_{\mathit{\text{g}}}{\mathit{\text{u}}}_{\mathit{\text{g}}}+{\mathit{\text{Z}}}_{\mathit{\text{d}}}{\mathit{\text{u}}}_{\mathit{\text{d}}}+\mathit{e}$$

where y is the vector of observed single-tree phenotypes, b is the vector of fixed effects with the associated design matrix X, ${\mathit{u}}_{\mathit{a}}$, ${\mathit{u}}_{\mathit{g}}$, and ${\mathit{u}}_{\mathit{d}}$ are the vectors of random additive genetic, nonadditive genetic, and design effects, respectively, with associated design matrices ${\mathit{Z}}_{\mathit{a}}$, ${\mathit{Z}}_{\mathit{g}}$, and ${\mathit{Z}}_{\mathit{d}}$, and e is the vector of residuals. Random effects are normally distributed following ${\mathit{\text{u}}}_{\mathit{\text{a}}}\sim N(0,\mathit{\text{H}}{\sigma }_a^2)$, ${\mathit{\text{u}}}_{\mathit{\text{g}}} \sim N(0,\mathit{\text{I}}{\sigma }_g^2)$, and ${\mathit{\text{u}}}_{\mathit{\text{d}}} \sim N(0,\mathit{\text{I}}{\sigma }_d^2)$, where I is the identity matrix, H is the average numerator relationship matrix between genotypes, and ${\sigma }_a^2$, ${\sigma }_g^2$, and ${\sigma }_d^2$ are variances for additive genetic, nonadditive genetic, and design effects, respectively. In this model, which is hereafter referred to as non-spatial, residuals are independent, following $\mathit{\text{e}} \sim N(0,\mathit{\text{I}}{\sigma }_e^2)$. To account for the spatial correlation between trees in the trials, we consider an alternative model which allocates a separable first-order autoregressive process component to residuals. In this case, $\mathit{\text{e}} \sim N(0,\mathit{\text{R}})$ where $\mathit{\text{R}}=[AR1({\rho }_{col})\otimes AR1({\rho }_{row})]{\sigma }_{\xi }^2+\mathit{\text{I}}{\sigma }_{\eta }^2$. $AR1(\rho )$ is a first-order autoregressive correlation matrix with autocorrelation parameter $\rho$, and ${\sigma }_{\xi }^2$ and ${\sigma }_{\eta }^2$ are spatial and independent residual variances, respectively.

In the first model instance, the average numerator relationship matrix H is the pedigree-based relationship matrix A calculated from the pedigree on all genotypes present in the trial. Alternatively, a blended matrix ${\mathit{H}}_{\mathit{\lambda }}$ was used for the variance structure of the additive genetic effects. ${\mathit{H}}_{\mathbf{\lambda }}$ was computed for all genotypes using A and the genomic marker-based relationship matrix G, calculated on genotyped individuals. A and G were calculated using the AGHmatrix R package (Amadeu et al., 2016), using the VanRaden method for G (VanRaden, 2008). The blended genetic relationship matrix ${\mathit{H}}_{\mathbf{\lambda }}$ was calculated using the ASRgenomics package (Gezan et al., 2022) and can be defined by its inverse as follows:

$${\mathbf{H}}_{\mathbf{\lambda }}-1={\mathbf{A}}^{-1}+\left[\begin{array}{cc}\mathbf{0}& \mathbf{0}\\ \mathbf{0}& \lambda ({\mathit{\text{G}}}^{-1}-{\mathit{\text{A}}}_{\mathbf{22}}-1)\end{array}\right]$$

where ${\mathit{A}}_{\mathbf{22}}$ is the genetic relationship matrix for genotyped individuals and λ is a scaling factor bound between 0 and 1 representing the weight associated to G relative to A for genotyped individuals.

Thorough model building and model selection were performed on severity scores. The Bayesian Information Criterion (BIC) as well as graphical checks of model diagnostics and variograms were used to guide model selection. One-sided likelihood ratio tests on random effect variances were used to compare models with nested random effects, and incremental addition or removal of fixed effects was guided by the Wald test. The model was constructed using the pedigree-based relationship matrix A. We started with a base model with an intercept, additive (${\mathit{u}}_{\mathit{a}}$) and nonadditive (${\mathit{u}}_{\mathit{g}}$) genetic effects, and block as a design effect (${\mathit{u}}_{\mathit{d}}$). Cubic spline transformations of trial rows and columns were added as fixed effects in b after visual examination of residuals against trial rows and columns, to account for environmental effects at the trial level. Two alternative spatial models for residuals were tested, one assuming no independent residual component (${\sigma }_{\eta }^2=0$), and one including the residual component. Once the most appropriate residual structure was selected, a reassessment of the relevance of block effects and row and column spline effects was conducted.

As part of the analysis of severity scores, after selecting the final set of fixed and random effects and the appropriate residual structure, we fitted the model using ${\mathit{H}}_{\mathit{\lambda }}$ and tested different values of λ in 0.1 increments across its possible range. The optimal λ value was selected based on model BIC. The final model including the optimal ${\mathit{H}}_{\mathit{\lambda }}$ matrix was implemented for the univariate genetic analyses of all remote sensing traits.

Narrow-sense heritability was calculated as follows:

$$\widehat{h^2}=\frac{\widehat{{\sigma }_a^2}}{\widehat{{\sigma }_a^2}+\widehat{{\sigma }_g^2}+\widehat{{\sigma }_r^2}}$$

where $\widehat{{\sigma }_r^2}=\widehat{{\sigma }_e^2}$ for non-spatial models and $\widehat{{\sigma }_r^2}=\widehat{{\sigma }_{\eta }^2}$ for spatial models.

Estimated breeding values were extracted from random genotype effects and their accuracy was obtained from univariate models using the following calculation:

$$accurac{y}_i= \sqrt{1-\frac{\widehat{S{E}_{i}2}}{\widehat{{\sigma }_a^2}}}$$

where $S{E}_i$ is the estimated standard error of the breeding value for genotype i.

Multi-trait analyses

Multi-trait analyses consist of simultaneously modelling several response traits and allow estimating the correlations between traits through structured variance matrices. In our case, the genetic correlation is estimated through appropriate use of the genetic relationship matrix. When the genetic correlation between two traits is non-zero, prediction performance of genotypes for one or more of the modelled traits can be improved. In our case, multi-trait analyses are advantageous when the trait in focus for selection has a lower heritability than other measured correlated traits. They can also be used when measurements for the focus trait are only available for a subset of observations in the trial: the trait can still be predicted from unmeasured units using the correlation with other traits that have been measured. The general multivariate model with d traits is as follows:

$$\mathit{\text{y}}=\mathit{\text{X}}\mathit{b}+{\mathit{\text{Z}}}_{\mathit{a}}{\mathit{\text{u}}}_{\mathit{a}}+{\mathit{\text{Z}}}_{\mathit{g}}{\mathit{\text{u}}}_{\mathit{\text{g}}}+{\mathit{\text{e}}}_{\mathit{\xi }}+{\mathit{\text{e}}}_{\mathit{\eta }}$$

where $\mathit{\text{y}}={[{\mathit{\text{y}}}_{\mathbf{1}},\dots ,{\mathit{\text{y}}}_{\mathit{\text{d}}}]}^T$ is the stacked column vector of phenotypic observations relating to each trait. Secondary traits ${\mathit{y}}_{\mathbf{2}},\dots ,{\mathit{y}}_{\mathit{d}}$ were rescaled to a similar scale as ${\mathit{y}}_{\mathbf{1}}$ if necessary, to improve convergence of the restricted maximum likelihood (REML) algorithm. $\mathit{\text{b}}={[{\mathit{\text{b}}}_{\mathbf{1}},\dots ,{\mathit{\text{b}}}_{\mathit{\text{d}}}]}^T$ is the column vector of fixed effects selected in univariate analyses (intercept and spline functions of trial rows and columns) with the associated design matrix $\mathit{\text{X}}=\left[\begin{array}{ccc}{\mathit{\text{X}}}_{\mathbf{1}}& \cdots & \mathbf{0}\\ ⋮& \ddots & ⋮\\ \mathbf{0}& \cdots & {\mathit{\text{X}}}_{\mathit{\text{d}}}\end{array}\right]$. ${\mathit{\text{u}}}_{\mathit{\text{a}}}={[{\mathit{\text{u}}}_{\mathit{\text{a}}\mathbf{1}},\dots ,{\mathit{\text{u}}}_{\mathit{\text{a}}\mathit{d}}]}^T$ and ${\mathit{\text{u}}}_{\mathit{\text{g}}}={[{\mathit{\text{u}}}_{\mathit{\text{g}}\mathbf{1}},\dots ,{\mathit{\text{u}}}_{\mathit{\text{g}}\mathit{d}}]}^T$ are the column vectors of random additive genetic and nonadditive genetic effects, respectively, with associated design matrices ${\mathit{\text{Z}}}_a=\left[\begin{array}{ccc}{\mathit{\text{Z}}}_{\mathit{\text{a}}\mathbf{1}}& \cdots & \mathbf{0}\\ ⋮& \ddots & ⋮\\ \mathbf{0}& \cdots & {\mathit{\text{Z}}}_{\mathit{\text{a}}\mathit{d}}\end{array}\right]$ and ${\mathit{\text{Z}}}_{\mathit{\text{g}}}=\left[\begin{array}{ccc}{\mathit{\text{Z}}}_{\mathit{\text{g}}\mathbf{1}}& \cdots & \mathbf{0}\\ ⋮& \ddots & ⋮\\ \mathbf{0}& \cdots & {\mathit{\text{Z}}}_{\mathit{\text{g}}\mathit{d}}\end{array}\right]$. In this model notation we separated the residuals into a spatially structured component ${\mathit{e}}_{\mathit{\xi }}$ and an independent component ${\mathit{e}}_{\mathit{\eta }}$, where ${\mathit{\text{e}}}_{\xi }={[{\mathit{\text{e}}}_{\xi \mathbf{1}},\dots ,{\mathit{\text{e}}}_{\xi \mathit{\text{d}}}]}^T$ and ${\mathit{\text{e}}}_{\eta }={[{\mathit{\text{e}}}_{\eta \mathbf{1}},\dots ,{\mathit{\text{e}}}_{\eta \mathit{\text{d}}}]}^T$. Random effects and residuals are modelled as follows: ${\mathit{u}}_{\mathit{g}} \sim N(0,{\mathbf{\sum }}_g \otimes {\mathbf{I}}_{\mathbf{p}})$, ${\mathbf{u}}_{\mathbf{a}} \sim N(0,{\Sigma }_a \otimes {\mathit{\text{H}}}_{\mathit{\text{λ}}})$, ${\mathit{\text{e}}}_{\xi } \sim N(0,\mathbf{ }{\Sigma }_{\xi }\otimes {\mathit{\text{R}}}_{\mathit{\text{ξ}}})$, and ${\mathit{e}}_{\mathit{\eta }} \sim N(0,{\mathit{\text{Σ}}}_{\mathit{\text{η}}}\otimes {\mathbf{I}}_{\mathbf{n}})$, where ${\Sigma }_{\mathit{g}}$, ${\Sigma }_{\mathit{a}}$, ${\Sigma }_{\mathit{\xi }}$, and ${\mathbf{\text{Σ}}}_{\mathit{\text{η}}}$ are general heterogenous variance structures of the form

$${\mathit{\sum }}_{\mathit{\text{x}}}=\left[\begin{array}{ccc}{\sigma }_x{{\displaystyle }}_1^2& \cdots & {\rho }_{x_{1d}}{\sigma }_{x_1}{\sigma }_{x_d}\\ ⋮& \ddots & ⋮\\ {\rho }_{x_{d1}}{\sigma }_{x_1}{\sigma }_{x_d}& \cdots & {\sigma }_x{{\displaystyle }}_d^2\end{array}\right]$$

${\mathit{I}}_{\mathit{n}}$ and ${\mathit{I}}_{\mathit{p}}$ are identity matrices of sizes n, the number of trees in the trial, and p, the number of genotypes represented in the trial. ${\mathit{H}}_{\mathit{\lambda }}$ is the blended genetic relationship matrix. Finally, ${\mathit{\text{R}}}_{\mathit{\text{ξ}}}=[AR1({\rho }_{col})\otimes AR1({\rho }_{row})]$ is the separable autoregressive variance structure of spatially correlated residuals. We note that this multivariate model has the same fixed effects, random effects, genetic relationship matrix and type of residual structure as the final univariate model selected for all studied traits. We also note that this model assumes the same spatial row and column autocorrelation values ${\rho }_{col}$ and ${\rho }_{row}$ across jointly modelled traits.

We first identified the remote sensing index with the highest observed correlation with severity scores. To estimate genetic correlations between the index and severity scores, and to assess whether joint modelling enhances the accuracy of estimated breeding values (EBVs) of severity scores, we performed a bivariate analysis with severity scores as the primary trait and the selected index as the secondary trait.

To test whether visual phenotyping can be reduced to a sample of trees in the trial and complemented by remote sensing data, we repeated and modified the previous bivariate analysis. Instead of using the complete set of severity scores as the primary trait, we randomly selected from 5% to 100% of severity scores and used the selected subset of severity scores as the primary trait. We repeated the procedure 10 times for each sampled proportion.

We then explored the use of remote sensing indices only in DNB severity phenotyping, by incrementally building a multi-trait model with a combination of high-heritability indices that are correlated with severity scores, starting with the most correlated index and adding additional indices in decreasing order of heritability values.

Results

Distribution and genetic analysis of severity scores

The distribution of DNB severity over the 2054 scored trees covered the full range of possible values, from 5% to 100%, with a median score of 55% (Figure 3). Tree-level and trial-level spatial dependence is apparent, (Figure 3), suggesting that tree-to-tree transmission of the disease and environmental heterogeneity in the trial had an influence on disease severity.

Figure 3

Figure 3. Histogram (left) of severity scores and their spatial distribution (right) across rows and columns in the trial. The coloured histogram (left) also acts as a colour legend for the spatial distribution of severity scores (right).

Selected genetic models of severity scores were in accordance with this observation. Model selection was performed using the pedigree-based relationship matrix A for the estimation of additive genetic effects. Including a spatial autocorrelation variance structure of residuals led to a significantly better fit than modelling independent residuals (Table 1), and the addition of splines for trial rows and columns in the fixed effects also improved the model. These results suggest both a landscape effect (variation across the trial site) and a neighbour effect (proximity related disease transmission) on the severity of disease. The addition of a spatially independent residual term gave a likelihood ratio test p-value of 0.05. Row and column correlations were around 0.5 and the spatially dependent error variance was five times higher than the independent error variance, suggesting that most of the residual structure could be accounted for by spatial factors. However, as the best model overall (lowest BIC) included both separable autoregressive and spatially independent error terms, both were included in the final model. As block effects became insignificant after accounting for spatial structure in the residuals, we omitted the block term from the final model, effectively pooling any block effect variance with the residual variance. Table 1 shows variance components and estimated genetic parameters for the best spatial and non-spatial models when using the A matrix.

Table 1

Table 1. Parameter estimates, narrow-sense heritability, Bayesian information criterion (BIC) and estimated breeding value (EBV) accuracy for four alternative models of Dothistroma needle blight severity scores.

The test on proportion of G vs. A matrices in the blended ${\mathit{H}}_{\mathit{\lambda }}$ matrix was performed using the model selected above. We found that values of λ between 0.4 and 0.6 were optimal. We set $\lambda =0.5$ as it provided both the lowest BIC and highest heritability for severity scores. The final model selected to calculate and compare breeding values across traditional and remote sensing traits was therefore the spatial model with the ${\mathit{H}}_{\mathbf{0.5}}$ matrix, and its variance components and estimated genetic parameters are presented in Table 1, together with results for the best non-spatial model for comparison purposes. Narrow-sense heritability of severity scores using ${\mathit{H}}_{\mathbf{0.5}}$ was 0.38 with a standard error (SE) of 0.06 when estimated from the non-spatial model and 0.57 (SE=0.08) when estimated from the final spatial model.

Univariate analysis of remote sensing indices

The same final univariate spatial model and ${\mathit{H}}_{\mathbf{0.5}}$ matrix that were used in the analysis of visual severity scores were also applied to model the genetic parameters of each remote sensing metric. Supplementary Table S1 reports narrow-sense heritability for each remote sensing index, estimated from variance components from non-spatial models with design effects. Calculating heritability on the non-spatial model with design effects makes values more comparable across indices as design effects are the same across indices and all the residual variance is included in the calculation. Narrow-sense heritability of remote sensing indices ranged from 0 to 0.37 when estimated from non-spatial models, reaching values similar to severity scores (0.38). Absolute observed correlations with severity scores reached 0.79 (Supplementary Table S1). The relationship between heritability (h²) and absolute Pearson’s correlation with severity scores (r_P) pictured in Figure 4A allows the identification of indices with high heritability and high absolute correlation with severity scores. These traits are potential candidates to assist or replace visual scoring in genetic selection for DNB resistance.

Figure 4

Figure 4. (A) Relationship between narrow-sense heritability estimated from univariate models of all remote sensing indices, and their absolute observed correlation with severity scores, by sensor type and derivation type. Note that the multispectral category corresponds to hyperspectral indices whose calculation is compatible with multispectral capture. The dashed line indicates the estimated narrow-sense heritability of severity scores. The box delineates the area of (B). RTM, Radiative transfer model; SIF, Solar-induced fluorescence. (B) Expansion of (A) for indices with heritability above 0.3, with point and index labels coloured by sensitivity to different physiological characteristics.

The most promising indices were mostly NBHIs (Figure 4B), some of which included spectral bands that can be calculated from multispectral instruments (Supplementary Table S1). The index RNIR_CRI700, RNIR_CRI550 and PSSR_c were the best performing indices (Figure 4B, Supplementary Table S1), which are proxies for variation in carotenoid content: RNIR_CRI550 and RNIR_CRI700 (h² = 0.37), are carotenoid reflectance indices (Gitelson et al., 2003, 2006). RNIR_CRI700 is also the NBHI with the highest absolute correlation with severity scores (r_P=-0.79). Many other important indices (CI, GM1, GM2, GM4, PSSR_a, PSSR_b, VOG3, BGI2, LIC3) characterise chlorophyll content and three were sensitive to both chlorophyll content and plant structure (MCARI3, MSR, SR). Notably, none of the water indices had either high heritability (h² = 0.004 - 0.04) or high correlation with disease (0.065 - 0.375). Mean and median normalised canopy temperature (TcTa_InMn and TcTa_InMed), derived from thermal captures, were highly correlated with severity scores but had a lower heritability. Solar-induced fluorescence only had a low to moderate correlation with severity scores and low heritability. Traits inverted from radiative transfer modelling were not as useful as the best NBHIs and mainly had heritability values between 0.15 and 0.28.

We selected the 21 indices with heritability values above 0.30 (Figure 4B) and computed their pairwise correlations using Pearson’s correlation coefficient (Figure 5). Pairwise correlations ranged from 0.71 to 1. Many of these NBHIs were highly correlated to each other as they are derived from similar wavelengths (Supplementary Table S1). For instance, RNIR_CRI550 and RNIR_CRI700 (h² = 0.37), are carotenoid reflectance indices (Gitelson et al., 2003, 2006) and were strongly correlated with each other (r_P > 0.99). RNIR_CRI700 is also the NBHI with the highest absolute correlation with severity scores (r_P=-0.79). Most other indices with heritability values > 0.30 were strongly correlated with RNIR_CRI700, with Pearson’s correlation coefficients ranging from 0.83 to 0.99.

Figure 5

Figure 5. Pairwise correlation between indices with heritability values above 0.30. The left-to-right and top-to-bottom order of indices reflects their ranking in decreasing heritability values.

Figure 6 shows the relationship between severity scores and RNIR_CRI700 both in terms of phenotypic observations (A) and EBVs from final spatial univariate models (B). Other promising indices included three Red/Green/Blue indices: the blue/green index (BGI2), the ratio analysis of reflectance spectra (RARS) and the Lichtenthaler index (LIC3). These indices use spectral bands compatible with multispectral capture. One high-heritability NBHI, the Sugar Beet Rust Index (SBRI), was developed for disease assessment (Mahlein et al., 2013). All other high-ranking indices are NBHI relating to pigments or structural characteristics, and their description can be found in Supplementary Table S1.

Figure 6

Figure 6. Relationship between severity scores and RNIR_CRI700. (A). Relationship between RNIR_CRI700 and observed severity. (B) Relationship between RNIR_CRI700 estimated breeding values (EBVs) and EBVs from the univariate model. (C) Relationship between RNIR_CRI700 EBVs from the bivariate model with SIPI and severity EBVs from the univariate model. A smoothing spline (blue line) has been fitted to aid visualisation and Pearson’s correlation coefficient is shown for all relationships.

Multi-trait analysis of severity scores

We performed a bivariate spatial analysis with the full set of severity scores as the primary trait and RNIR_CRI700 as the secondary trait, as its absolute correlation with severity scores was the highest (r_P = -0.79), and its relationship with severity scores appeared to be relatively linear (Figure 6A). The bivariate model estimated a genetic correlation of -0.80 (SE 0.09) between severity scores and RNIR_CRI700. The accuracy of severity score EBVs was on average 0.75. This represents little improvement from the mean accuracy of 0.74 obtained from the univariate analysis of severity scores, which is not surprising as severity scores have a higher estimated heritability than remote sensing indices.

Results from all instances of the bivariate model with various proportions of scored trees as the primary trait and RNIR_CRI700 as the secondary trait are shown in Figure 7. When 80% or more trees are scored, mean EBV accuracy of severity scores were similar or lower than the EBV accuracy of severity scores from the univariate model, and steadily decreased to around 0.5 with 5% trees scored (Figure 7A). The precision and accuracy of genetic correlation estimates decreases as fewer trees are scored, with estimated values becoming unreliable below 40% of sampled trees (Figure 7B). Severity score EBVs are highly correlated to EBVs from the univariate model with full scoring, even at low scoring proportions (Figure 7C). However, a sharper decrease in correlation is observed when fewer than 25% trees are scored, with correlation values dropping to below 0.85.

Figure 7

Figure 7. (A) Accuracy of estimated breeding values (EBVs) for severity scores for each proportion of severity scores sampled. Grey points represent each iteration, black points and error bars are the average and standard deviation for all 10 iterations of a given scoring proportion. The dashed line is the accuracy of severity scores EBVs from the univariate model using the full set of severity scores. (B) Estimated additive genetic correlation between severity scores and RNIR_CRI700 for each iteration of the simulation. (C) Pearson’s correlation coefficient between EBVs of sampled severity scores in each iteration of the bivariate model and EBVs from the univariate model using the full set of severity scores.

Multi-trait analysis of remote sensing indices

To assess the ability of remote sensing indices to fully replace severity scores in DNB resistance phenotyping, we selected the index most correlated with severity scores, RNIR_CRI700, as the primary trait in the spatial multi-trait analysis of remote sensing metrics. We selected remote sensing indices with estimated heritability higher than 0.2 and an absolute correlation with RNIR_CRI700 lower or equal to 0.90 as potential secondary traits. We avoided selecting indices with a higher absolute correlation as they would add little information compared to the univariate model of the primary trait and lead to multi-trait model convergence issues. Resulting EBV accuracies and EBV correlation with severity scores from a univariate model are presented in Table 2.

Table 2

Table 2. Observed correlation, estimated genetic correlation, mean estimated breeding value accuracy and correlation with severity scores resulting from bivariate models with RNIR_CRI700 as the primary trait.

Modelling RNIR_CRI700 with the structure-intensive pigment index SIPI as the secondary trait increased its EBV accuracy by 6.8%, from 0.71 to 0.76. Modelling RNIR_CRI700 with the blue/green index BGI2 as the secondary trait increased its EBV accuracy by 4.3%, from 0.71 to 0.74. Using any other index as a secondary trait led to a decrease in EBV accuracy from the univariate analysis of RNIR_CRI700. The absolute correlation between RNIR_CRI700 EBVs and severity score EBVs ranged from 0.55 to 0.79. Figure 6C shows the relationship between RNIR_CRI700 EBVs from the bivariate analysis with SIPI, which has a Pearson correlation coefficient of -0.76 with severity scores. This correlation is less strong than the correlation calculated with RNIR_CRI700 EBVs from the univariate analysis (r_P = -0.79, Figure 6B), which suggests that jointly modeling traits adds complexity without necessarily improving the alignment of RNIR_CRI700 with severity scores. Using RNIR_CRI700 alone may therefore be more straightforward and equally effective.

Discussion

This study leverages a Dothistroma needle blight outbreak displaying the full range of severity levels within a radiata pine clonal breeding trial to extend phenotyping methods for the assessment of disease susceptibility using remote sensing techniques. We quantified the extent of genetic control of novel thermal and spectral disease susceptibility phenotypes through a large set of replicated genotypes in a clonal breeding trial. Our research focused on identifying the most appropriate phenotypes and statistical methodologies for estimating genetic parameters and tested several operational scenarios for phenotyping. These scenarios included approaches solely based on remote sensing and those integrating remote sensing with traditional phenotyping techniques. Although a few studies have calculated narrow-sense heritability of spectral indices in tree species (Čepl et al., 2018; Tao et al., 2021), this is the first time that such indices have been used to characterise disease susceptibility in an operational breeding trial, and the first time field-based thermal measurements on forest stands have been used for tree genetic assessment. The successful combination of these measurement methods with UAV technology paves the way for a future implementation of high-throughput phenotyping of plant stress in forestry breeding trials.

The narrow-sense heritability values for visual severity scores obtained in this study align with previous estimates from breeding trials in New Zealand and Australia (Li et al., 2018; Ismael et al., 2020; Klápště et al., 2020). We confirmed that foliar disease modelling requires the addition of spatial structure through both trial-level geographical information (here trial rows and columns) and neighbourhood effects to take into account the landscape variability and mode of transmission of the disease. The application of single-step genomic evaluation, integrating pedigree and genomic marker information in the genetic relationship matrix, improved model performance. The optimal proportion of the G matrix, which reflects the weighting between pedigree and genomic marker information, was found to be between 0.4-0.6, which is consistent with previous studies within the same breeding programme (Klápště et al., 2020, 2022).

UAV hyperspectral and thermal data acquisition at the time of the outbreak enabled the remote phenotyping of disease severity through the calculation of spectral, thermal, and fluorescence indices. Pigment- and structure-sensitive narrow-band hyperspectral indices had a strong observed correlation with severity scores and a narrow-sense heritability similar to severity scores. The importance of these indices is consistent with previous research describing the impacts of DNB and disease on changes in tree structure and fundamental leaf-level physiological processes. Previous research has demonstrated the importance of carotenoids as a symptom of DNB severity (Watt et al., 2023b) and the carotenoid:chlorophyll ratio has been found to increase as a result of increasing disease severity (Watt et al., 2023b). This is consistent with our results of carotenoid-sensitive indices such as RNIR_CRI700, RNIR_CRI550, and PSSR_c showing highest correlation with severity scores. It was interesting to note that the RTM plant trait that performed best was carotenoid content (CAR). Although CAR did not perform as well as the pigment-related NBHIs discussed above, the relative importance of CAR amongst other RTM traits is consistent with rankings between indices. The utility of the important PRI indices (PRI_CI, PRI515, PRI550, PRIm1), which can account for xanthophyll dynamics, are consistent with a large body of previous research demonstrating their utility for detecting a range of stresses (Peguero-Pina et al., 2008; Zarco-Tejada et al., 2009; Watt et al., 2024a). Photochemical reflectance index is an effective proxy for the reductions in stomatal conductance and assimilation rate that are associated with pre-visual and early disease (Hernández-Clemente et al., 2019) and was found to be one of the most important classifiers of DNB severity within a different radiata pine field trial (Watt et al., 2023b). It is interesting to note that non-standard formulations of PRI appeared to be more important than the standard PRI which uses R₅₃₁ as a xanthophyll-sensitive spectral band. Compared to the standard formulation, the PRI variants using bands from the 500–515 nm spectral region (i.e., PRI_m1 and PRI₅₁₅) have been shown to be less sensitive to leaf area index, tree densities, and structural effects in conifer canopies (Hernández-Clemente et al., 2011), which highlights again the importance of accounting for stand structure in identifying changes in pigments.

Canopy temperature indices from thermal captures showed promising performance in this study. Following energy balance theory, increases in leaf temperature result from reductions in stomatal conductance and lower rates of transpiration (Smith et al., 1986; Still et al., 2019), and our results show that these modifications are detectable by UAV thermal imagery at the canopy level.

Through an extensive simulation exercise, we showed that the joint modelling of severity scores and the best-performing NBHI, RNIR_CRI700, allows the estimation of severity scores breeding values with very high fidelity (absolute EBV correlation above 0.9) when at least 50% of the trial is manually scored. Even when only 20% of the trial is scored, the fidelity remains relatively high (absolute EBV correlation above 0.8). The associated loss in accuracy is minimal, only decreasing from 0.74 to 0.68 at 50% scoring and to 0.59 at 20% scoring. We found that using the remote sensing trait RNIR_CRI700 as the sole disease severity phenotype in a univariate analysis led to an EBV accuracy of 0.71 with an absolute correlation of 0.79 with severity scores. These results suggest that partial visual scoring could be entirely replaced by hyperspectral remote sensing, with little loss in EBV correlation and no reduction in EBV accuracy compared to a bivariate analysis that incorporates partial visual scoring. We also assessed the potential benefits of combining several indices to improve phenotyping from only one index. However, using a combination of remote sensing indices in a multi-trait setup did not improve the quality of genetic results over univariate results, as it only minorly increased EBV accuracy and considerably reduced the correlation of the selected remote sensing index EBVs with severity score EBVs. This can be explained by the redundancy of many of the pigment indices. Indeed, although we identified 21 promising pigment-related indices, we found a high level of correlation among them, which can partially be explained by the similarity in their calculation from spectral bands (Supplementary Table S1). We therefore suggest that using the remote sensing trait most correlated with traditional severity scores in a univariate analysis can be sufficient and leads to EBV estimations more in line with severity score EBVs than those estimated from the joint analysis of multiple remote sensing traits.

The focus of this study was restricted to exploring remote sensing options to supplement or replace currently time-consuming disease phenotyping methods, with the objective to maximise high fidelity of alternative methods to traditional ones. Although the technical requirements of high-quality spectral and thermal field captures presented here might currently be prohibitive in the forestry context, we found that many of the high-heritability indices use spectral bands that can be captured with multispectral instruments, which are more accessible than hyperspectral ones. As UAV remote-sensing assessments become more routine and gain trust among foresters and tree breeders, we believe the advantages of remote sensing for tree disease phenotyping will expand and outweigh traditional methods. First, the precision, accuracy, and consistency of remote-sensing measurements may lead to more effective breeding outcomes than with traditional phenotyping. Also, structural traits captured in years following disease using LiDAR or spectral methods or derived through radiative transfer modelling such as leaf area index can complement hyperspectral and thermal data to assess longer-term outcomes of disease on tree health such as defoliation. As a number of stress responses detected by hyperspectral and thermal metrics are likely common to several different diseases, remote sensing phenotyping of stress responses might be useful in simultaneously selecting for resistance to multiple foliar diseases. Ismael et al. (2020) showed evidence for strong genetic correlation of resistance to multiple foliar diseases such as Dothistroma needle blight and Cyclaneusma needle cast in radiata pine. More generally, hyperspectral and thermal indices can be used not only for disease-related phenotypes but also for other biotic and non-biotic stress responses. For instance, canopy temperature can specifically be used to characterise drought stress phenotypes, as it is a proxy for stomatal conductance, transpiration and assimilation rate.

Conclusions

Our study illustrates the potential for the remote sensing capture of hyperspectral and thermal traits to be used in disease resistance phenotyping for tree breeding. Genotype performance for traditional and remote sensing phenotypes was enhanced by the use of genomic and pedigree data through spatially explicit single-step predictive modelling and led to high breeding value accuracy. Single narrow-band hyperspectral indices calculated on all trees in the breeding trial can be used to complement partial visual disease scoring or replace it altogether with relatively high fidelity, with close to no loss in breeding value accuracy. However, the success of high-throughput single-tree characterisation of hyperspectral and thermal traits in a field trial setup achieved in this study is highly innovative and required a large amount of expertise and resources. Additional advanced experimentation in field methodologies will help streamline and reduce the cost of data capture processes. This will enable the extension of this disease phenotyping approach to other trials and other phenotyping applications and solidify the adequate identification of remote sensing metrics as phenotypes for physiological traits of interest.

Data availability statement

The datasets generated for this study can be found in the Dryad data repository: DOI: 10.5061/dryad.5hqbzkhhk.

Author contributions

JE: Conceptualization, Data curation, Formal analysis, Investigation, Methodology, Visualization, Writing – original draft, Writing – review & editing. RM: Data curation, Investigation, Methodology, Visualization, Writing – original draft, Writing – review & editing. RH: Data curation, Methodology, Resources, Writing – review & editing. MW: Conceptualization, Investigation, Project administration, Resources, Writing – original draft, Writing – review & editing.

Funding

The author(s) declare that financial support was received for the research and/or publication of this article. This research was funded by the Ministry of Business, Innovation and Employment (MBIE), as part of the Endeavour programme entitled “Seeing the forest for the trees: transforming tree phenotyping for future forests” (grant number C04X2101).

Acknowledgments

We thank forest managers for Kinleith Forest, Manulife Investment Management Forest Management NZ Ltd, for permission to access the site. We thank the Radiata Pine Breeding Company (RPBC) for providing access to breeding trials to capture UAV data, DNB severity scores and genomic data (SNPchip), which has been instrumental to perform this study. We thank the Scion Autonomous Systems team for their great effort conducting the UAV data captures, especially Peter Massam, Honey Jane Estarija, Warren Yorston, and David Cajes. Thanks also to Andrew MacMillan (Manaaki Whenua - Landcare Research) for capturing and processing the hyperspectral UAV data. We are also grateful to Tomas Poblete, Na Wang, and Pablo J. Zarco-Tejada (School of Agriculture, Food and Ecosystem Sciences (SAFES), University of Melbourne) who calculated the radiative transfer modelling-based plant traits used in this study. We are deeply grateful to Sadeepa Jayathunga for her processing of the LiDAR data captured at this site, which has been instrumental to this and other studies. We thank Dalila Pasquini, Jack Moss and Andrew Macalister for their assistance in the field. Natalie Graham and Sai Arojju (RPBC) kindly provided assistance with access to genomic data. Sai Arojju provided useful advice on the analyses conducted in this study.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Generative AI statement

The author(s) declare that no Generative AI was used in the creation of this manuscript.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fpls.2025.1574720/full#supplementary-material

References

Amadeu, R. R., Cellon, C., Olmstead, J. W., Garcia, A. A., Resende, M. F., and Munoz, P. R. (2016). AGHmatrix: R Package to Construct Relationship Matrices for Autotetraploid and Diploid Species: A Blueberry Example. Plant Genome 9 (3). doi: 10.3835/plantgenome2016.01.0009

PubMed Abstract | Crossref Full Text | Google Scholar

Barnes, I., Crous, P. W., Wingfield, B. D., and Wingfield, M. J. (2004). Multigene phylogenies reveal that red band needle blight of Pinus is caused by two distinct species of Dothistroma, D. septosporum and D. pini. Stud. Mycol. 50, 551–565.

Google Scholar

Bassett, C., Buchanan, M., Gallagher, R., and Hodges, R. (1970). A toxic difuroanthraquinone from Dothistroma pini. Chemistry & Industry 52, 1659–1660

Google Scholar

Bradshaw, R. (2004). Dothistroma (red-band) needle blight of pines and the dothistromin toxin: a review. For. Pathol. 34, 163–185. doi: 10.1111/j.1439-0329.2004.00356.x

Crossref Full Text | Google Scholar

Bulman, L. S., Dick, M. A., Ganley, R. J., McDougal, R. L., Schwelm, A., and Bradshaw, R. E. (2013). “Dothistroma needle blight,” in Infectious forest diseases (Wallingford UK: CABI), 436–457.

Google Scholar

Bulman, L., Gadgil, P., Kershaw, D., and Ray, J. (2004). Assessment and control of Dothistroma needle-blight. Forest Research Bulletin 229, 1–48.

Google Scholar

Butler, D. (2013). On the optimal design of experiments under the linear mixed model. PhD Thesis (School of Mathematics and Physics, The University of Queensland).

Google Scholar

Butler, D., Gilmour, A., and Gogel, B. (2009). ASReml-R reference manual (Brisbane: The State of Queensland: Department of Primary Industries and Fisheries).

Google Scholar

Carson, S. (1989). Selecting Pinus radiata for resistance to Dothistroma needle blight. New Z. J. Forest. Sci. 19, 3–21.

Google Scholar

Čepl, J., Stejskal, J., Lhotáková, Z., Holá, D., Korecký, J., Lstibůrek, M., et al. (2018). Heritable variation in needle spectral reflectance of Scots pine (Pinus sylvestris L.) peaks in red edge. Remote Sens. Environ. 219, 89–98. doi: 10.1016/j.rse.2018.10.001

Crossref Full Text | Google Scholar

Corbin, J. P., Best, R. J., Garthwaite, I. J., Cooper, H. F., Doughty, C. E., Gehring, C. A., et al. (2024). Hyperspectral leaf reflectance detects interactive genetic and environmental effects on tree phenotypes, enabling large-scale monitoring and restoration planning under climate change. Plant Cell Environ 48 (3), 1842–1857. doi: 10.1111/pce.15263

PubMed Abstract | Crossref Full Text | Google Scholar

Dungey, H. S., Dash, J. P., Pont, D., Clinton, P. W., Watt, M. S., and Telfer, E. J. (2018). Phenotyping whole forests will help to track genetic performance. Trends Plant Sci. 23, 854–864. doi: 10.1016/j.tplants.2018.08.005

PubMed Abstract | Crossref Full Text | Google Scholar

Féret, J.-B., Gitelson, A., Noble, S., and Jacquemoud, S. (2017). PROSPECT-D: Towards modeling leaf optical properties through a complete lifecycle. Remote Sens. Environ. 193, 204–215. doi: 10.1016/j.rse.2017.03.004

Crossref Full Text | Google Scholar

Gezan, S., Avelar, A., Galli, G., and Murray, D. (2022). ASRgenomics: An R package with complementary genomic functions. Ed. International, V. (United Kingdom: Hemel Hempstead).

Google Scholar

Gitelson, A. A., Gritz, Y., and Merzlyak, M. N. (2003). Relationships between leaf chlorophyll content and spectral reflectance and algorithms for non-destructive chlorophyll assessment in higher plant leaves. J. Plant Physiol. 160, 271–282. doi: 10.1078/0176-1617-00887

PubMed Abstract | Crossref Full Text | Google Scholar

Gitelson, A. A., Keydan, G. P., and Merzlyak, M. N. (2006). Three-band model for noninvasive estimation of chlorophyll, carotenoids, and anthocyanin contents in higher plant leaves. Geophys. Res. Lett. 33, 1–5. doi: 10.1029/2006GL026457

Crossref Full Text | Google Scholar

Graham, N., Telfer, E., Frickey, T., Slavov, G., Ismael, A., Klápště, J., et al. (2022). Development and validation of a 36K SNP array for radiata pine (Pinus radiata D.Don). Forests 13, 1–14. doi: 10.3390/f13020176

Crossref Full Text | Google Scholar

Hernández-Clemente, R., Hornero, A., Mottus, M., Penuelas, J., González-Dugo, V., Jiménez, J. C., et al. (2019). Early diagnosis of vegetation health from high-resolution hyperspectral and thermal imagery: lessons learned from empirical relationships and radiative transfer modelling. Curr. Forest. Rep. 5, 169–183. doi: 10.1007/s40725-019-00096-1

Crossref Full Text | Google Scholar

Hernández-Clemente, R., Navarro-Cerrillo, R. M., Suárez, L., Morales, F., and Zarco-Tejada, P. J. (2011). Assessing structural effects on PRI for stress detection in conifer forests. Remote Sens. Environ. 115, 2360–2375. doi: 10.1016/j.rse.2011.04.036

Crossref Full Text | Google Scholar

Isik, F., Holland, J., and Maltecca, C. (2017). Genetic data analysis for plant and animal breeding (Cham, Switzerland: Springer).

Google Scholar

Ismael, A., Suontama, M., Klapste, J., Kennedy, S., Graham, N., Telfer, E., et al. (2020). Indication of quantitative multiple disease resistance to foliar pathogens in pinus radiata D.Don in New Zealand. Front. Plant Sci. 11. doi: 10.3389/fpls.2020.01044

PubMed Abstract | Crossref Full Text | Google Scholar

Kabir, M., Ganley, R., and Bradshaw, R. (2015). Dothistromin toxin is a virulence factor in dothistroma needle blight of pines. Plant Pathol. 64, 225–234. doi: 10.1111/ppa.2015.64.issue-1

Crossref Full Text | Google Scholar

Khosravipour, A., Skidmore, A. K., Isenburg, M., Wang, T., and Hussin, Y. A. (2014). Generating pit-free canopy height models from airborne lidar. Photogram. Eng. Remote Sens. 80, 863–872. doi: 10.14358/PERS.80.9.863

Crossref Full Text | Google Scholar

Klápště, J., Dungey, H. S., Graham, N. J., and Telfer, E. J. (2020). Effect of trait’s expression level on single-step genomic evaluation of resistance to Dothistroma needle blight. BMC Plant Biol. 20, 1–13. doi: 10.1186/s12870-020-02403-6

PubMed Abstract | Crossref Full Text | Google Scholar

Klápště, J., Ismael, A., Paget, M., Graham, N. J., Stovold, G. T., Dungey, H. S., et al. (2022). Genomics-enabled management of genetic resources in radiata pine. Forests 13, 1–19. doi: 10.3390/f13020282

Crossref Full Text | Google Scholar

Li, Y., Dungey, H. S., Carson, M., and Carson, S. (2018). Genotype by environment interaction for growth and Dothistroma resistance and clonal connectivity between environments in radiata pine in New Zealand and Australia. PloS One 13, e0205402. doi: 10.1371/journal.pone.0205402

PubMed Abstract | Crossref Full Text | Google Scholar

Liu, Y., Feng, H., Fan, Y., Yue, J., Yang, F., Fan, J., et al. (2025). Utilizing UAV-based hyperspectral remote sensing combined with various agronomic traits to monitor potato growth and estimate yield. Comput. Electron. Agric. 231, 109984. doi: 10.1016/j.compag.2025.109984

Crossref Full Text | Google Scholar

Mahlein, A.-K., Rumpf, T., Welke, P., Dehne, H.-W., Plümer, L., Steiner, U., et al. (2013). Development of spectral indices for detecting and identifying plant diseases. Remote Sens. Environ. 128, 21–30. doi: 10.1016/j.rse.2012.09.019

Crossref Full Text | Google Scholar

Maier, S. W., Günther, K. P., and Stellmes, M. (2004). “Sun-induced fluorescence: A new tool for precision farming,” in Digital imaging and spectral techniques: Applications to precision agriculture and crop physiology, vol. 66. (American Society of Agronomy), 207–222.

Google Scholar

Oerke, E.-C. (2020). Remote sensing of diseases. Annu. Rev. Phytopathol. 58, 225–252. doi: 10.1146/annurev-phyto-010820-012832

PubMed Abstract | Crossref Full Text | Google Scholar

Pan, J., Lin, J., and Xie, T. (2023). Exploring the potential of UAV-based hyperspectral imagery on pine wilt disease detection: influence of spatio-temporal scales. Remote Sens. 15, 2281. doi: 10.3390/rs15092281

Crossref Full Text | Google Scholar

Peguero-Pina, J. J., Morales, F., Flexas, J., Gil-Pelegrín, E., and Moya, I. (2008). Photochemistry, remotely sensed physiological reflectance index and de-epoxidation state of the xanthophyll cycle in Quercus coccifera under intense drought. Oecologia 156, 1–11. doi: 10.1007/s00442-007-0957-y

PubMed Abstract | Crossref Full Text | Google Scholar

Plascyk, J. A. (1975). The MK II fraunhofer line discriminator (FLD-II) for airborne and orbital remote sensing of solar-stimulated luminescence. Opt. Eng. 14, 339–330. doi: 10.1117/12.7971842

Crossref Full Text | Google Scholar

Plowright, A. (2025). ForestTools: Tools for Analyzing Remote Sensing Forest Data. R package version 1.0. 3.

Google Scholar

Poland, J. A. and Nelson, R. J. (2011). In the eye of the beholder: the effect of rater variability and different rating scales on QTL mapping. Phytopathology 101, 290–298. doi: 10.1094/PHYTO-03-10-0087

PubMed Abstract | Crossref Full Text | Google Scholar

R Core Team (2023). R: A Language and Environment for Statistical Computing (Vienna, Austria: Computing).

Google Scholar

Smigaj, M., Agarwal, A., Bartholomeus, H., Decuyper, M., Elsherif, A., de Jonge, A., et al. (2024). Thermal infrared remote sensing of stress responses in forest environments: a review of developments, challenges, and opportunities. Curr. Forest. Rep 10 (1), 56–76. doi: 10.1007/s40725-023-00207-z

Crossref Full Text | Google Scholar

Smigaj, M., Gaulton, R., Suárez, J. C., and Barr, S. L. (2019). Combined use of spectral and structural characteristics for improved red band needle blight detection in pine plantation stands. For. Ecol. Manage. 434, 213–223. doi: 10.1016/j.foreco.2018.12.005

Crossref Full Text | Google Scholar

Smith, R., Heritage, A., Stapper, M., and Barrs, H. (1986). Effect of stripe rust (Puccinia striiformis West.) and irrigation on the yield and foliage temperature of wheat. Field Crops Res. 14, 39–51. doi: 10.1016/0378-4290(86)90045-6

Crossref Full Text | Google Scholar

Song, P., Wang, J., Guo, X., Yang, W., and Zhao, C. (2021). High-throughput phenotyping: Breaking through the bottleneck in future crop breeding. Crop J. 9, 633–645. doi: 10.1016/j.cj.2021.03.015

Crossref Full Text | Google Scholar

Still, C., Powell, R., Aubrecht, D., Kim, Y., Helliker, B., Roberts, D., et al. (2019). Thermal imaging in plant and ecosystem ecology: applications and challenges. Ecosphere 10, e02768. doi: 10.1002/ecs2.2019.10.issue-6

Crossref Full Text | Google Scholar

Tao, X., Li, Y., Yan, W., Wang, M., Tan, Z., Jiang, J., et al. (2021). Heritable variation in tree growth and needle vegetation indices of slash pine (Pinus elliottii) using unmanned aerial vehicles (UAVs). Industr. Crops Prod. 173, 114073. doi: 10.1016/j.indcrop.2021.114073

Crossref Full Text | Google Scholar

Telfer, E., Graham, N., Macdonald, L., Li, Y., Klapste, J., Resende, M., Jr., et al. (2019). A high-density exome capture genotype-by-sequencing panel for forestry breeding in Pinus radiata. PloS One 14, e0222640. doi: 10.1371/journal.pone.0222640

PubMed Abstract | Crossref Full Text | Google Scholar

VanRaden, P. M. (2008). Efficient methods to compute genomic predictions. J. dairy Sci. 91, 4414–4423. doi: 10.3168/jds.2007-0980

PubMed Abstract | Crossref Full Text | Google Scholar

Vapnik, V. (1999). The nature of statistical learning theory (New York: Springer science & business media).

Google Scholar

Verhoef, W. and Bach, H. (2007). Coupled soil–leaf-canopy and atmosphere radiative transfer modeling to simulate hyperspectral multi-angular surface reflectance and TOA radiance data. Remote Sens. Environ. 109, 166–182. doi: 10.1016/j.rse.2006.12.013

Crossref Full Text | Google Scholar

Virlet, N., Costes, E., Martinez, S., Kelner, J.-J., and Regnard, J.-L. (2015). Multispectral airborne imagery in the field reveals genetic determinisms of morphological and transpiration traits of an apple tree hybrid population in response to water deficit. J. Exp. Bot. 66, 5453–5465. doi: 10.1093/jxb/erv355

PubMed Abstract | Crossref Full Text | Google Scholar

Watt, M. S., Bartlett, M., Soewarto, J., de Silva, D., Estarija, H. J. C., Massam, P., et al. (2023a). Previsual and early detection of myrtle rust on rose apple using indices derived from thermal imagery and visible-to-short-infrared spectroscopy. Phytopathology® 113, 1405–1416. doi: 10.1094/PHYTO-02-23-0078-R

PubMed Abstract | Crossref Full Text | Google Scholar

Watt, M. S., Estarija, H. J. C., Bartlett, M., Main, R., Pasquini, D., Yorston, W., et al. (2024a). Early detection of myrtle rust on pōhutukawa using indices derived from hyperspectral and thermal imagery. Remote Sens. 16, 1050. doi: 10.3390/rs16061050

Crossref Full Text | Google Scholar

Watt, M. S., Jayathunga, S., Hartley, R. J. L., Pearse, G. D., Massam, P. D., Cajes, D., et al. (2024b). Use of a consumer-grade UAV laser scanner to identify trees and estimate key tree attributes across a point density range. Forests 15, 899. doi: 10.3390/f15060899

Crossref Full Text | Google Scholar

Watt, M. S., Poblete, T., De Silva, D., Estarija, H. J. C., Hartley, R. J. L., Leonardo, E. M. C., et al. (2023b). Prediction of the severity of Dothistroma needle blight in radiata pine using plant based traits and narrow band indices derived from UAV hyperspectral imagery. Agric. For. Meteorol. 330, 109294. doi: 10.1016/j.agrformet.2022.109294

Crossref Full Text | Google Scholar

Zarco-Tejada, P. J., Berni, J. A., Suárez, L., Sepulcre-Cantó, G., Morales, F., and Miller, J. R. (2009). Imaging chlorophyll fluorescence with an airborne narrow-band multispectral camera for vegetation stress detection. Remote Sens. Environ. 113, 1262–1275. doi: 10.1016/j.rse.2009.02.016

Crossref Full Text | Google Scholar

Zarco-Tejada, P. J., Poblete, T., Camino, C., González-Dugo, V., Calderon, R., Hornero, A., et al. (2021). Divergent abiotic spectral pathways unravel pathogen stress signals across species. Nat. Commun. 12, 6088. doi: 10.1038/s41467-021-26335-3

PubMed Abstract | Crossref Full Text | Google Scholar

Zhang, H., Li, W., Wang, G., Song, F., Wen, Z., Zhang, H., et al. (2025). Predicting stomatal conductance of chili peppers using TPE-optimized LightGBM and SHAP feature analysis based on UAVs’ hyperspectral, thermal infrared imagery, and meteorological data. Comput. Electron. Agric. 231, 110036. doi: 10.1016/j.compag.2025.110036

Crossref Full Text | Google Scholar