Genotype, Environment, Year, and Harvest Effects on Fruit Quality Traits of Five Blueberry (Vaccinium corymbosum L.) Cultivars

Abstract

Blueberries (Vaccinium spp.) comprise a broad range of perennial woody species. Introgression of native species into cultivated germplasm has adapted Vaccinium germplasm to a range of climates and growing conditions for cultivated blueberry. Genetic differences signify phenotypic variance that is observed among blueberry accessions. In addition, variability in geographic and climatic growing conditions between environments or within the same environment across different years may further affect fruit and plant phenotypic expression. As a result, a phenotype is a function of genetic background (G), environment (E), and their interaction (G × E). In addition, other temporally regulated factors such as year (Y) and harvest time (H) impact plant and fruit quality phenotypic variation. Our research aimed to assess the genotypic performance of five blueberry cultivars, including ‘Echota’, ‘O’Neal’, ‘Reveille’, ‘Summit’, and ‘Sunrise’. The selected cultivars were phenotyped for various fruit quality-related traits over two sequential harvests in two years and two locations. Our results indicated that genotype was a significant source of variation for most phenotypic characteristics. Further, the effect of Y × H and G × Y × H significantly affected the majority of studied phenotypic traits. Within the studied genotypes, ‘Reveille’ and ‘O’Neal’ phenotypic stability were consistent across locations and years; additionally, ‘Summit’ phenotypic characteristics were stable across years, environments, and harvests. Clonal plant replicates within a genotype, harvest, and environment, in addition to individual fruit measures, were the most significant sources of variability.

1. Introduction

Within North America, the geographic range of native Vaccinium species ranges from Canada through Mexico, encompassing a broad spectrum of different climatic growing conditions. In the southeastern United States, select native Vaccinium species including V. angustifolium Aiton, V. arboreum Marsh, V. constablaei Gray, V. corymbosum L., V. darrowii Camp, V. elliottii Chapm., V. myrtilloides Michx., V. pallidum Aiton, V. tenellum Aiton, V. simulatum Small, V. stamineum L., and V. virgatum Aiton have been utilized in breeding populations for specific adaptations and trait introgression. Introgression of native species into cultivated germplasm has adapted Vaccinium germplasm to a range of climates and growing conditions for cultivated blueberry, resulting in both low-chill and high-chill adapted cultivars [1,2]). Breeding for low-chill and high-chill cultivation has enabled blueberry production expansion, increasing the geographic range and facilitating fresh market availability year-round due to longer production windows across regions. Furthermore, blueberry production has gained traction due to increased consumer and commercial attention due to recent reports regarding the health benefits of blueberry total anthocyanin content and anthocyanin composition, including aglycone, sugar moiety, and acyl group [3].

Genetic differences signify phenotypic variance and are observed as fruit flavor or size differences between different blueberry accessions. Variability in geographic and climatic growing conditions between environments may further affect fruit and plant phenotypic expression. A comparative study of geographically diverse locations found that blueberry cv. ‘Brightwell’ (V. virgatum) grown in higher altitude environments reported increased total soluble solids (TSS), flavonoids, phenols, proanthocyanidins, and anthocyanins in China [4]. In addition to altitude, environmental differences encompass climatic factors, including temperature and precipitation. Increased TSS can be attributed to decreased fruit volume and moisture loss resulting from reduced rainfall events [5]. Additionally, decreased precipitation during blueberry fruit ripening resulted in small, firm fruit, which correlated with increased soluble solids and acidity [6]. Water deficit treatments also decreased fruit antioxidant levels in blueberry [7]. Alternatively, increased TSS and decreased fruit weight resulted from warmer temperatures in tomato [8]. In addition, temperature also plays a role in anthocyanin content. Genetic control of apple anthocyanin biosynthesis was found to be regulated by temperature [9]. Hot climates decreased the expression of anthocyanin biosynthesis genes, which failed to accumulate red fruit pigmentation in apples grown in hot climates compared to those grown under temperate climate conditions [9].

Phenotypic expression and observed variation of plant growth and development are functions of both genetic background (G) and environment (E) and their interaction (G × E) [10]. Temporal variation can also impact plant and fruit quality phenotypic variation. Temporal regulation includes the effects of the year (Y) and harvest time (H), wherein biotic and abiotic events between years or maturation time between harvests are involved in fruit production and quality differences. Seasonal temperature differences over the years can significantly affect the fruit production of blueberries [11,12,13]. Winter temperatures control dormancy release. Thus, the earlier onset of warmer temperatures exemplified by the three-week earlier snowmelt in Vaccinium myrtillus in Norway in 2011 resulted in two-week earlier blooming compared with the 2010 production season [12]. Furthermore, severe weather events such as spring cold spells can cause significant cold injury to low-chill and precocious blueberry cultivars, impacting yield [14]. In addition to cold-related events, yearly differences in temperature fluctuations and precipitation events, including heatwaves, drought, and flooding, affect developing fruits’ quality [15,16]. Within a harvest season, weekly climate changes and increased growing degree days (GDD) can impact blueberry fruit maturation rates between sequential harvests within the same plant; exemplified by the and cause observed variation in phenolic, acid, and total soluble solid composition between harvests [17]. Within a harvest season, weekly climate changes and increased growing degree days (GDD) can impact blueberry fruit maturation rates between sequential harvests within the same plant; exemplified by the and cause observed variation in phenolic, acid, and total soluble solid composition between harvests [17]. By studying the impact of these interactions on phenotypic traits, we can evaluate factors affecting consistency in genotype performance.

In studying the effects of environment, years, harvests, and genotype on phenotypic traits, our research aimed to assess genotypic performance in select environments and evaluate phenotypic consistency. The cultivars used in this study are highbush blueberry ‘Echota’, ‘O’Neal’, ‘Reveille’, ‘Summit’, and ‘Sunrise’ (Supplementary Figure S1). The southern highbush (SHB) cultivars ‘O’Neal’ and ‘Summit’ have V. darrowii in their pedigrees; ‘Summit’ and ‘Reveille’ also have low-chill V. virgatum introgression; ‘Reveille’ has V. tenellum in its pedigree history. ‘Echota’ is a high-chill northern highbush (NHB) derived from a native low-chill North Carolina V. corymbosum accession (NC102). ‘Sunrise’ is a modern northern highbush cultivar developed from high-chill V. corymbosum parents with pedigrees containing V. angustifolium; as such, ‘Sunrise’ does not have low-chill species in its pedigree history. Using these selected cultivars in both North Carolina and Oregon environments, our research informs on the stability of fourteen phenotypic traits. The phenotypic stability of a trait across environments, years, and harvests indicates a greater emphasis on the genetic control of the trait, which is essential for blueberry genotype uniformity for commercial production. The objectives of this study were to evaluate the genetic consistency of five blueberry cultivars and assess them for phenotypic variation between two environments, two years of research, two sequential harvests per year, and their interactions. Further, we aimed to evaluate the correlation between phenotypic traits. The results of this study may impact blueberry breeding program decision-making for target environments and further identify blueberry production benchmarks for specified traits.

2. Materials and Methods

2.1. Plant Material

Fresh fruits were collected from selected cultivars from two field locations, NC State University (NCSU) Sandhills Research Station located in Jackson Springs, NC (35°11′ N, 79°40′ W and 176 m above sea level) and USDA-ARS National Clonal Germplasm Repository (NCGR) in Corvallis, OR (44°33′ N, 123°13′ W and 74 m above sea level), weekly from May through August in 2019 and 2020. The Jackson Springs soil is Fuquay sand and is amended by pine bark for acidification. The Jackson Springs location was irrigated with a Valley Linear Precision Irrigation system every 7–10 days. The Corvallis site soil is predominantly Malabon silty clay loam with a pH of 4.5–5.5 that is regularly amended with organic sawdust mulching. The Corvallis NCGR farm was irrigated via an underground sprinkler system that was run twice a week.

The studied genotypes consisted of Vaccinium corymbosum SHB, ‘O’Neal’, ‘Reveille’, and ‘Summit’, as well as the NHB ‘Echota’ and ‘Sunrise’. Each cultivar was represented by two clonal replicate plants planted 1 m apart in the field at each fixed location. The fruit ripening window was based on fruit maturity on a per genotype and location basis. Ripe fruits from each plant were harvested for two sequential weeks once the plant achieved a ripeness status of >30% blue fruit. A minimum of ten individual fruits were hand-harvested per clonal replicate of each genotype for each harvest. Fruit collected from Jackson Springs were stored at 4 °C for next-day analysis. Fruit collected at the Corvallis location were harvested and shipped overnight on ice and in refrigerated shipping conditions (4 °C) to the NC State University, Raleigh Campus, Blueberry Genetics and Genomics Laboratory for next-day analysis.

Ten fruit from each clonal replicate were randomly sampled per harvest at each location. The same random sample of fruit was used for firmness, weight, diameter, and puncture measurements; these fruits were later homogenized into a fruit purée for quantitative measurement of the phenotypic traits of TSS and TA. Pooling the fruit samples for colorimetric, TSS, and TA quantitative measures was conducted to ensure sufficient material for their respective procedures.

2.2. Phenotypic Trait Measurements

2.2.1. Colorimetric Analysis

Color indices were quantified for 40 mL volumetric samples of each genotype, replicate, and harvest, encased in a light-impermeable Zero Calibration Box (CM-A124 Konica Minolta). CIELAB values, L*, a*, and b*, were quantified using a Konica Minolta CR-5 Chroma Meter (Konica Minolta, Inc., Tokyo, Japan). Three automatic readings were taken per sample, and the instrument recorded the average of the three readings. Quantified values indicated spectrum position for color identification. Specifically, higher L* values indicate increased sample luminescence represented by light blue fruit; a* values are a red-green color spectrum where positive values are red and negative values are green; b* values are a blue-yellow color spectrum where positive values are yellow and negative values are blue. Luminescence (L*) values depreciate with fruit handling and bloom removal, observed by darker fruit skin and elevated a* and b* values. Therefore, CIELAB values were measured first in the phenotypic workflow prior to excessive fruit handling. The a* and b* values were converted to chroma (CRM) and hue (HUE) values according to the following equations [18].

$$\begin{array}{c}Chroma=\sqrt{\left(a^{*}{}^2+b^{*}{}^2\right)}\\ \theta =\frac{{\tan }^{-1}\frac{b^{*}}{a^{*}} }{6.2832}\ast 360\\ hue angle=\theta if a^{*}>0 and b^{*}\ge 0 \\ hue angle=\theta +180 if a^{*}<0\\ hue angle=\theta +360 if a^{*}>0 and b^{*}<0\end{array}$$

2.2.2. Fruit Firmness Measurements

Fresh fruit firmness (FRM) was quantified via a Firmtech II (BioWorks, KS, USA). The instrument load cell depression speed was set at 12 mm·s⁻¹ with a table speed of 0.79 mm·s⁻¹. Ten fruits were selected and measured per clonal genotype replicate and harvest. Individual fruits were positioned for equilateral compression by the Firmtech II instrument, which recorded fruit deflection force. Deflection force was quantified as grams necessary to compress the fruit 1 mm (g∙mm⁻¹). Deflection force and compression arm height were calibrated with pure gum rubber balls, 15.4 mm in diameter, with a tensile strength of 2700 psi (MCMasterCarr part 96385K61) [19].

2.2.3. Fruit Weight and Size

Fruit weight (FW) was taken from these randomly selected individual fresh fruit (n =10) using a standard laboratory scale (Mettler Toledo, Columbus, OH). The square root of FW was taken and used in statistical analysis to stabilize the variance.

The diameter of the ten randomly selected and previously weighed fruits was measured via caliper from calyx to stem scar, gauging equatorial diameter (ED) and polar diameter (PD). The roundness index (RI) was determined as the ratio of polar diameter to the equatorial diameter for each fruit [20].

$$RI=\frac{PD}{ED}$$

2.2.4. Texture Analysis

The mechanical force required to puncture the fruit skin of the same ten individual fruits used in diameter and weight measurements was determined using a TA.XTplus texture analyzer (Stable MicroSystem Ltd., Godalimng, UK). The texture analyzer was fitted with a 4 mm flathead probe as previously described [21]. The instrument’s decompression speed was set to 2.5 mm·s⁻¹ with a retraction speed of 5 mm·s⁻¹. The flathead probe applied a maximum force of 5 g, indicating completion of fruit puncture and platform contact, after which probe retraction was initiated. Consistent with firmness measurements, individual fruits were positioned equatorially for puncture. Exponent Connect software (Stable MicroSystem Ltd., Godalimng, UK) was used to produce graphical profiles for each fruit calculating absolute positive force (g) (APF), force at target (g) (FT), and distance at absolute positive force (mm) (DPF). APF is the pressure exerted upon the fruit at the time of puncture (Supplementary Figure S2). Likewise, the FT is the innate resistance to depression of the fruit upon initial contact with the flathead probe (Supplementary Figure S2). The DPF indicates the fruit’s elasticity (FE), measuring the distance from initial target contact to APF and fruit puncture. Using DPF divided by the ED of the blueberry, we established a ratio of FE.

$$FE=\frac{DPF}{ED}$$

Total soluble solids (TSS) were measured from the supernatant of the fruit purée. The supernatant was extracted from homogenized (Fisher Scientific Homogenizer 150, Fisher Scientific, Pittsburg, PA) and centrifuged (10 min 4200 RPM/10,120× g with HIGHConic II Rotor ofSorvall Legend X1R Centrifuge, Fisher Scientific, Pittsburg, PA) fruit pulp of fresh fruit purée at room temperature (n =10) using a handheld ATAGO PAL-BX|ACID F5 refractometer (ATAGO, USA, Inc., Bellevue, WA, USA). The fruit was contained in a 50 mL centrifuge tube held in ice during homogenization to prevent purée heating and soluble solid decomposition [22]. A two-hundred microliter volume of liquid supernatant was extracted and placed on the refractometer sensor for TSS measurement, and the TSS values were recorded.

2.2.5. Titratable Acidity

Titratable acidity (TA) was measured using a Mettler Toledo G20S Compact Titrator (Mettler Toledo, Columbus, OH). The fruit was homogenized and vortexed, creating a blueberry fruit purée. Subsequently, 2 ± 0.01 g of the purée was mixed with 60 mL of water and subsequently titrated with 0.1 N NaOH to an endpoint of pH 8.2. Titratable acidity was converted to citric acid percentage using the following equation.

$$Titratable Acidity \%=\frac{mL NaOH\times 0.1 N NaOH\times milliequivalent factor\times 100}{sample \left(g\right)}$$

2.2.6. Climatic Data

Monthly averages were extrapolated from maximum and minimum daily air and soil temperatures, as well as precipitation and evapotranspiration measures. These data were obtained during blueberry flower and fruit development, from March to August, using State Climate Office of North Carolina, Sandhill Research Station weather data, and the Hyslop weather station data maintained by Oregon State University (Corvallis, OR, USA).

2.2.7. Statistical Analysis

Statistical analyses for phenotypic data were performed using the REML method of Proc GLIMMIX in SAS v9.4 (SAS, Cary, NC, USA). All phenotypic traits were checked for variance stability and heteroskedasticity. Subsequently, stabilizing transformation was performed only on FW data, which was transformed by taking the square root of each data point. Environment (E) was analyzed as a blocking factor for five identical genotypes (G) between the Corvallis and Jackson Springs to evaluate statistical significance between sites for each phenotypic variable. Model 1 included environment (E), year (Y), genotype (G), harvest (H), random effects of replication (R), and their interactions. Replication refers to duplicate plants of each genotype and measures the variability among plants of the same genotype grown 1 m apart. Individual fruit replication within each clonal replicate plant for each measured phenotypic trait was used to estimate within-plant variability for these analyses:

Model 1

$$\begin{array}{l}{y}_{ijklmn}=\text{µ}+E_i+Y_j+{\epsilon }_{1_{ij}}+G_k+E{G}_{ik}+Y{G}_{jk}+R{\left(GE\right)}_{m\left(ki\right)}+EY{G}_{ijk}+R{\left(GEY\right)}_{m\left(kij\right)}+H_l+E{H}_{il}+Y{H}_{jl}\\ +H{G}_{lk}+YG{H}_{jkl}+R{\left(GHE\right)}_{m\left(lki\right)}+R{\left(GHY\right)}_{m\left(lkj\right)}+EG{H}_{ikl}+EY{H}_{ijl}+EYG{H}_{ijkl}\\ +R{\left(GEHY\right)}_{m\left(lkij\right)}+{\epsilon }_{2ijklmn}\end{array}$$

where $y_{ijklmn}$ is the measured phenotypic trait of the ith environment in the jth year and lth sequential harvest for the kth genotype and the mth clonal plant replicate per genotype and the nth individual fruit. In this model, µ is the grand mean, and $E_i, Y_j, G_k, \mathrm{and} H_l$ are the main effects of environment, year, genotype, and harvest, respectively, where i = 1, 2 environments; j = 1, 2 years; k = 1, 2, 3, 4, 5 genotypes; l = 1, 2 harvests; and m = 1, 2 clonal genotype replicates. Two-factor interactions included $E{G}_{ik}, Y{G}_{jk}, E{H}_{il}, Y{H}_{jl},\mathrm{and}H{G}_{lk}$, representing the interactions between environment and genotype, year and genotype, environment and harvest, year and harvest, and harvest and genotype, respectively. Three-factor interactions included $EY{G}_{ijk},EG{H}_{ikl}, \mathrm{and} EY{H}_{ijl}$ for environment-by-year-by-genotype, environment-by-genotype-by-harvest, and environment-by-year-by-harvest. The four-factor interactions included only $EYG{H}_{ijkl}$, representing the interaction for environment-by-year-by-genotype-by-harvest. Clonal genotype replicates were present in both environments and were treated as a random effect. The random effects included whole plot error, ${\epsilon }_{1_{ij}}$, of the E × Y interaction, random between-plant effects of the clonal replicate plants (all terms involving R_m), and the error between replicate measures of individual fruit for the phenotypic trait.

A second model (model 2) was developed for statistical analysis of pooled fruit. The phenotypic measurements on fruit purée comprised of the ten individual fruit per genotype, replicate, harvest, location, and year. This model included environment (E), year (Y), genotype (G), harvest (H), and their interactions while using the random effect of replication (R) for accounting for variability between plants:

Model 2

$$\begin{array}{c}{y}_{ijklm}=\text{µ}+E_i+Y_j+{\epsilon }_{1_{ij}}+G_k+L{G}_{ik}+Y{G}_{jk}+R{\left(GE\right)}_{m\left(ki\right)}+GE{Y}_{ijk}+R{\left(GEY\right)}_{m\left(kij\right)}+H_l+E{H}_{il}+Y{H}_{jl}+H{G}_{lk}+GY{H}_{jkl}+R{\left(GEH\right)}_{m\left(lki\right)}\\ +R{\left(GHY\right)}_{m\left(lkj\right)}+GE{H}_{ikl}+EY{H}_{ijl}+GEY{H}_{ijkl}+{\epsilon }_{2ijlkm}\end{array}$$

where $y_{ijklm}$ is the measured phenotypic trait of the ith environment in the jth year and lth sequential harvest for the kth genotype and the mth clonal plant replicate per genotype. In this model, µ is the grand mean and $E_i, Y_j, G_k, \mathrm{and} H_l$ are the main effects of environment, year, genotype, and harvest, respectively, where i = 1, 2 environments; j = 1, 2 years; k = 1, 2, 3, 4, 5 genotypes; l = 1, 2 harvests; and m = 1, 2 clonal genotype replicates. Two-factor interactions included $E{G}_{ik}, Y{G}_{jk}, E{H}_{il}, Y{H}_{jl},\mathrm{and}H{G}_{lk},$ representing the interactions between environment and genotype, year and genotype, environment and harvest, year and harvest, and harvest and genotype, respectively. Three-factor interactions included $EY{G}_{ijk},EG{H}_{ikl}, \mathrm{and} EY{H}_{ijl}$ for environment-by-year-by-genotype, environment-by-genotype-by-harvest, and environment-by-year-by-harvest. The four-factor interaction included only $EYG{H}_{ijkl}$, representing the interaction for environment-by-year-by-genotype-by-harvest. Clonal genotype replicates were present in both environments and were treated as a random effect. Similar to Model 1, the random effects included whole plot error, ${\epsilon }_{1_{ij}}$, of the E × Y interaction and the error attributed to the interaction of the clonal genotypic replicate within genotype and environment with year and harvest.

Correlation analyses were performed for genotypes within each location using the R “psych” package [23] (Revelle 2021), were analyzed by the Pearson’s method [24], and were considered significant at a threshold of p < 0.05.

3. Results

3.1. Variance across Environments, Years, and Harvests

Analysis of variance for individual fruit-measured phenotypic traits indicated that genotype (G), environment (E), year (Y), harvest (H), and their interactions significantly influenced (p < 0.05) the majority of the evaluated traits. Model 1 was used for nine phenotypic traits that used individual fruit measurements, including FW, PD, ED, RI, FRM, APF, DPF, FT, and FE. After taking measures for FW, PD, ED, RI, FRM, APF, DPF, FT, and FE, the ten fruit were homogenized into a fruit purée for quantitative measurement of the five phenotypic traits analyzed using model 2. Collectively homogenizing ten fruit was performed so that measures including L*, CRM, HUE, TSS, and TA had sufficient material. Within model 1, environment (E) significantly impacted phenotypic traits, including RI and DPF (

Table 1. Analysis of variance (ANOVA) p-value summary of nine phenotypic traits evaluated using model 1 in PROC GLIMMIX (SAS v9.4) with five blueberry cultivars from two environments collectedover sequential harvests in 2019 and 2020.

Factor	FW †	PD	ED	RI	FRM	APF	DPF	FT	FE
Environment (E)	0.214	0.722	0.077	0.030‡	0.222	0.104	0.025	0.368	0.069
Year (Y)	0.161	0.542	0.684	0.651	0.001	0.976	0.107	0.679	0.471
Genotype (G)	0.0008	0.016	0.002	0.003	0.0002	0.006	0.468	0.010	0.005
G × E	0.360	0.083	0.673	0.543	0.031	0.083	0.005	0.039	0.209
G × Y	0.143	0.015	0.003	0.070	0.366	0.434	0.023	0.767	0.688
G × E × Y	0.038	0.526	0.029	0.070	0.761	0.684	0.074	0.416	0.616
Harvest (H)	0.002	0.005	0.0004	0.001	0.929	0.299	0.881	0.849	0.061
E × H	0.318	0.085	0.229	0.378	0.879	0.388	0.137	0.690	0.172
Y × H	0.007	0.023	0.471	0.046	0.175	0.556	0.343	0.491	0.108
G × H	0.092	0.088	0.185	0.521	0.845	0.827	0.601	0.754	0.741
E × Y × H	0.0001	0.010	0.043	0.406	0.105	0.322	0.158	0.267	0.197
G × E × H	0.167	0.415	0.273	0.295	0.859	0.353	0.685	0.611	0.403
G × Y × H	0.094	0.0002	<0.0001	0.011	0.790	0.120	0.990	0.493	0.007
G × E × Y × H	0.042	0.054	0.057	0.122	0.108	0.458	0.570	0.146	0.464

^† FW = fruit weight; PD = polar diameter; ED = equatorial diameter; RI = roundness index; FRM = firmness; APF = absolute positive force; DPF = distance at absolute positive force; FT = force at target; FE = fruit elasticity. ^‡ Treatment effects in bold are significant at p ≤ 0.05.

), where values were 4.9% greater and 4.8% lower, respectively, in the Corvallis environment (

Table 2. Means of fourteen phenotypic traits across environments, years, harvests, and genotypes.

	FW †	PD	ED	RI	FRM	APF	DPF	FT	FE	L*	HUE	CRM	TSS	TA
	Environment
Jackson Springs	1.19 a	10.49 a	13.98 a	0.76 b‡	180.90 a	321.33 a	5.65 a	265.04 a	28.36 a	29.55 a	265.40 a	4.80 a	13.02 b	0.96 a
Corvallis	1.16 a	10.54 a	13.33 a	0.80 a	192.27 a	371.66 a	5.38 b	292.74 a	36.28 a	27.79 b	266.16 a	4.15 b	14.30 a	0.60 b
	Year
Year 1	1.17 a	10.47 a	13.56 a	0.78 a	174.73 b	348.31 a	5.63 a	283.89 a	32.27 a	29.92 a	264.72 b	4.43 a	14.23 a	0.73 b
Year 2	1.18 a	10.56 a	13.81 a	0.77 a	198.04 a	340.82 a	5.42 a	271.60 a	31.76 a	27.51 b	266.82 a	4.57 a	12.97 b	0.87 a
	Harvest
Harvest 1	1.21 a	10.76 a	14.21 a	0.76 b	186.74 a	336.82 a	5.51 a	280.03 a	29.48 a	28.90 a	265.87 a	4.59 a	13.17 b	0.80 a
Harvest 2	1.13 b	10.25 b	13.11 b	0.79 a	185.57 a	353.23 a	5.54 a	275.65 a	34.82 a	28.55 b	265.62 a	4.41 a	14.09 a	0.79 a
	Genotype
‘Echota’	1.27	10.55	14.39	0.73	167.42	293.38	5.63	233.16	20.75	30.49	265.28	5.19	12.34	1.05
‘O’Neal’	1.20	11.03	13.60	0.82	189.13	337.53	5.29	269.96	28.20	26.21	265.58	3.83	13.87	0.54
‘Reveille’	0.95	9.49	11.54	0.83	224.26	469.91	5.42	374.81	52.18	28.61	266.95	4.49	15.50	0.73
‘Summit’	1.40	11.40	15.62	0.73	166.92	288.91	5.74	231.03	25.51	29.84	266.13	4.79	12.88	0.88
‘Sunrise’	1.13	10.33	13.72	0.76	178.84	320.2	5.58	268.76	31.86	28.66	264.91	4.22	13.31	0.8

^† FW = fruit weight (g); PD = polar diameter (mm); ED = equatorial diameter (mm); RI = roundness index; FRM = firmness (g·mm⁻¹); APF = absolute positive force (g); DPF = distance at absolute positive force (mm); FT = force at target (g); FE = fruit elasticity; L* = luminescence; HUE = hue angle; CRM = chroma; TSS = total soluble solids (%); TA = titratable acidity (%). ^‡ Means followed by the same letter(s) within a column under subheading of environment, year, harvest, and genotype are not significantly different using the least squares means (LSMEANS) Tukey HSD multiple comparisons procedure (p < 0.05).

). In model 2, E was a significant source of variance for phenotypic traits including L*, CRM, TSS, and TA (

Table 3. Analysis of variance (ANOVA) p-value summary of five phenotypic traits evaluated on fruit purée using model 2 in PROC GLIMMIX (SAS v9.4) with five blueberry cultivars from two environments collected over sequential harvests in 2019 and 2020.

Factor	L* †	HUE	CRM	TSS	TA
Environment (E)	0.0002 ‡	0.117	0.007	0.001	<0.0001
Year (Y)	<0.0001	<0.0001	0.745	<0.0001	0.006
Genotype (G)	0.0001	0.133	0.002	<0.0001	0.003
G × E	0.391	0.061	0.678	0.012	0.137
G × Y	0.001	0.013	0.146	0.069	0.005
G × E × Y	0.003	0.495	0.305	0.087	0.001
Harvest (H)	0.049	0.940	0.146	0.001	0.777
E × H	0.649	0.002	0.135	0.402	0.410
Y × H	<0.0001	<0.0001	<0.0001	<0.0001	<0.0001
G × H	0.945	0.576	0.202	0.513	0.710
E × Y × H	<0.0001	<0.0001	<0.0001	<0.0001	<0.0001
G × E × H	0.039	0.110	0.389	0.619	0.198
G × Y × H	<0.0001	<0.0001	<0.0001	0.003	<0.0001
G × E × Y × H	<0.0001	<.0001	<0.0001	<0.0001	<0.0001

^† L* = luminescence; HUE = hue angle; CRM = chroma; TSS = total soluble solids; TA = titratable acidity. ^‡ Treatment effects in bold are significant at p ≤ 0.05.

). The fruit had significantly higher L*, CRM, and TA values at the Jackson Springs site, while the fruit from Corvallis had a higher TSS percentage (Table 2).

Only FRM in model 1 and L*, HUE, TSS, and TA in model 2 were impacted by the effect of year (Y) (Table 1 and Table 3). Overall, the first-year, L* and TSS measurements were >8% higher than the second year. In the second year, the measured FRM and TA were more than 13% of the first year. Compared to year 1, the HUE angle increased from 264.7 to 266.8 in year 2 (Table 2) (Supplementary Figure S3).

The difference between genotypes (G) was significant for all measured traits except for DPF and HUE (Table 1 and Table 3). For example, ‘Reveille’ had the lowest FW, PD, and ED, but the highest RI, FRM, APF, FT, FE, and TSS (Table 2, Figure 1 and Figure 2). Conversely, ‘Summit’ had the highest measured FW, PD, and ED and the lowest measured RI, FRM, and FT. Correspondingly, FW and ED have a moderate negative correlation with RI, FRM, APF, and FT (R² > |0.40|) (Figure 3). On the other hand, ‘O’Neal’ had the lowest L* and CRM values, while ‘Echota’ had the highest L* and CRM values of the evaluated cultivars (Table 2). As expected, the L* and CRM had a strong positive correlation (R = 0.69) (Figure 3).

Compared with the G × E (location) interaction, G × Y (time) affected more phenotypic traits, including PD, ED, RI, and DPF in model 1 and L* in model 2. In contrast, the G × E interaction was only significant for FRM, DPF, FT, and TSS. Combining the sources of variation in the three-factor interaction of G × E × Y, FW, RI, and DPF were significant in model 1, and TSS and TA in model 2 (Table 1 and Table 3).

More so than either E or Y, Harvest (H) was a significant source of variation for both models 1 and 2. In total, H was significant for nine phenotypic traits, including FW, PD, ED, and RI in model 1, and L*, HUE, CRM, TSS, and TA in model 2 (Table 1 and Table 3). FW, PD, ED, L*, CRM, TA, and HUE were higher in the first harvest compared with the second harvest. In contrast, the second harvest increased in RI and TSS (see Section 3.2).

The effect of E on H was less significant as a source of variation in the model compared to the impact of Y on H. The interaction of E × H was only a significant source of variation in HUE in model 2 (Table 3). Overall, HUE decreased between harvests in Jackson Springs, while Corvallis increased between harvests ( see Section 3.2). Y × H had a more significant impact on phenotypic traits used in both models compared with E × H. Y × H was significant for phenotypic traits including FW, PD, and RI in model 1 and L*, HUE CRM, TSS, and TA in model 2 (Table 1 and Table 3). Phenotypic traits including FW, PD, and CRM decreased between sequential harvests in both years of study; conversely, RI and TSS increased between sequential harvests in both years of study (Tables 6 and 7). Both L* and HUE increased between harvests in the first year and decreased between harvests in the second year, whereas TA decreased between harvests in the first year, while increasing in the second year. The effect of G on H (G × H) was not a significant source of variance to any of the studied phenotypic traits (Table 1).

Complex three-way interactions included E × Y × H, G × E × H, G × Y × H, and G × E × Y × H. Among these, E × Y × H was significant for phenotypic traits including FW, PD, ED, L*, HUE, CRM, TSS, and TA; G × E × H was significant for L*; G × E × Y × H was significant for phenotypic traits including FW, L*, HUE, CRM, TSS, and TA. Though G × H was not significant in the models for any phenotypic trait, G × Y × H was a significant three-way interaction for twelve of the fourteen measured traits, including PD, ED, RI, FE, L*, HUE, CRM, TSS, and TA (Table 1 and Table 3).

More scrutiny of the variance components contributing to the total variance in the two models revealed that the most significant variation in model 1 was the within-plant error corresponding to individual fruit variation, which accounted for >42% for all phenotypic traits (

Table 4. Variance component percentage of total variance of nine phenotypic traits evaluated on fresh fruits using model 1 in PROC GLIMMIX (SAS v9.4) with five blueberry cultivars from two environments collected over sequential harvests in 2019 and 2020.

Covariance Factor	FW †	PD	ED	RI	FRM	APF	DPF	FT	FE
Whole Plot Error: Environment (E) × Year (Y)	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00
Genotype (G) × E × Replicate (R)	6.93	7.89	19.25	20.72	0.00	16.81	0.00	9.47	14.29
G × E × Y × R	8.95	12.61	0.00	7.51	11.32	0.00	2.31	9.78	0.00
Harvest (H) × G × E × R	15.18	16.12	10.69	1.94	6.23	19.56	0.00	4.70	20.45
G × E × Y × H × R	22.91	15.64	26.37	10.81	10.98	12.16	3.79	6.37	23.05
Error among individual fruit	46.03	47.74	43.69	59.02	71.47	51.47	93.89	69.68	42.21

^† FW = fruit weight; PD = polar diameter; ED = equatorial diameter; RI = roundness index; FRM = firmness; APF = absolute positive force; DPF = distance at absolute positive force; FT = force at target; FE = fruit elasticity.

). The dominant variance components were dependent upon the phenotypic trait of study, often with greater emphasis on G × E × replicate (R) and G × E × Y × H × R (Table 4). In model 2, G × E × H × R accounted >45% of the variation in all phenotypic traits (

Table 5. Variance component percentage of the total variance of five phenotypic traits evaluated on fruit purée using model 2 in PROC GLIMMIX (SAS v9.4) with five blueberry cultivars from two environments collected over sequential harvests in 2019 and 2020.

Covariance Factor	L* †	HUE	CRM	TSS	TA
Whole Plot Error: Environment (E) × Year (Y)	1.69	0.00	3.70	0.00	3.13
Genotype (G) × E × Replicate (R)	0.00	0.00	0.00	0.00	11.86
G × E × Y × R	33.09	27.12	39.38	27.21	25.61
Harvest (H) × G × E × R	54.33	51.80	46.45	67.03	46.23
G × E × Y × H × R	10.88	21.08	10.48	5.77	13.17

^† L* = luminescence; HUE = hue angle; CRM = chroma; TSS = total soluble solids; TA = titratable acidity.

). The variance component G × E × Y × R accounted for 25–40% of model 2 variation.

3.2. Genotypic Means within Environment and Year

Among the genotypes, ‘Reveille’ had the most significant variability for nine measured traits across environmental locations (

Table 6. Harvest means for each of the five cultivars within a year, environment, and genotype for fourteen phenotypic traits using model 1 in PROC GLIMMIX (SAS v9.4).

Year	Harvest	FW †	PD	ED	RI	FRM	APF	DPF	FT	FE
		‘Echota’
		Corvallis, Oregon
2019	Harvest 1	1.34 a	11.58 a	15.53 a	0.75 a	142.97 a	270.30 a	5.97 a	206.75 a	17.48 a
	Harvest 2	1.52 a	12.00 a	15.22 a	0.79 a	131.06 a	227.83 a	5.26 a	163.86 a	14.98 a
2020	Harvest 1	1.19 a	10.79 a	13.99 a	0.77 a	166.60 a	281.07 a	5.40 a	235.67 a	20.64 a
	Harvest 2	1.14 a	10.50 a	13.33 a	0.79 a	182.62 a	291.97 a	5.08 a	220.62 a	22.14 a
		Jackson Springs, NC
2019	Harvest 1	1.23 a	9.44 a	14.65 a	0.65 a	168.53 a	330.29 a	6.14 a	264.91 a	22.78 a
	Harvest 2	1.15 a	9.62 a	14.14 a	0.68 a	157.65 a	323.79 a	6.39 a	231.95 a	23.57 a
2020	Harvest 1	1.20 a	9.66 a	13.68 a	0.71 a	205.20 a	303.02 a	5.36 a	262.55 a	22.25 a
	Harvest 2	1.27 a	10.79 a	14.54 a	0.74 a	184.74 a	318.74 a	5.43 a	324.79 a	22.19 a
		‘Echota’
		Corvallis, Oregon
2019	Harvest 1	1.28 a	11.96 a	14.32 a	0.84 a	181.65 a	389.54 a	5.78 a	315.97 a	32.86 a
	Harvest 2	1.35 a	12.06 a	14.42 a	0.83 a	190.79 a	359.96 a	5.14 a	285.25 a	31.64 a
2020	Harvest 1	1.23 a	10.93 a	13.10 a ‡	0.84 a	226.74 a	374.00 a	5.25 a	303.14 a	28.89 a
	Harvest 2	0.79 a	8.64 a	10.44 b	0.83 a	204.19 a	350.23 a	5.70 a	285.49 a	41.82 a
		Jackson Springs, NC
2019	Harvest 1	1.20 a	10.77 a	13.94 a	0.78 b	154.16 a	272.86 a	5.20 a	220.53 a	19.83 a
	Harvest 2	0.75 b	8.43 a	9.70 b	0.87 a	186.32 a	367.22 a	4.99 a	273.83 a	39.64 a
2020	Harvest 1	1.47 a	12.68 a	16.32 a	0.78 a	195.61 a	316.28 a	5.16 a	275.67 a	19.49 a
	Harvest 2	1.37 a	11.80 a	15.29 a	0.77 a	179.64 a	276.51 a	5.30 a	207.55 a	18.24 a
		‘Echota’
		Corvallis, Oregon
2019	Harvest 1	0.89 b	8.63 a	10.46 a	0.82 a	226.29 a	491.59 a	5.61 a	453.87 a	57.43 a
	Harvest 2	1.05 a	9.67 a	11.51 a	0.84 a	239.00 a	600.44 a	5.34 a	464.11 a	63.38 a
2020	Harvest 1	1.11 a	10.79 a	13.06 a	0.83 a	273.17 a	528.53 a	5.63 a	482.06 a	49.23 a
	Harvest 2	0.59 b	7.54 b	8.64 b	0.87 a	255.82 a	654.81 a	4.71 a	433.02 a	96.83 a
		Jackson Springs, NC
2019	Harvest 1	0.98 a	9.80 a	12.10 a	0.82 b	218.29 a	437.75 a	4.96 a	318.57 a	44.86 a
	Harvest 2	0.79 a	9.32 a	10.47 a	0.90 a	193.19 a	376.66 a	5.14 a	271.20 a	40.96 a
2020	Harvest 1	1.05 a	9.99 a	12.53 a	0.80 a	184.30 a	344.19 a	5.90 a	301.33 a	34.85 a
	Harvest 2	1.09 a	10.01 a	13.25 a	0.76 a	207.19 a	343.81 a	6.04 a	288.90 a	34.36 a
		‘Echota’
		Corvallis, Oregon
2019	Harvest 1	1.52 a	11.36 a	16.58 a	0.69 a	144.93 a	227.42 a	4.96 a	173.15 a	20.06 a
	Harvest 2	1.36 a	11.00 a	14.68 a	0.75 a	170.80 a	286.73 a	4.83 a	219.61 a	26.20 a
2020	Harvest 1	1.31	11.23	15.54	0.72	178.34	278.19	4.43	166.12	24.97
	Harvest 2 §
		Jackson Springs, NC
2019	Harvest 1	1.35 a	11.19 a	15.47 a	0.73 a	148.54 a	302.83 a	6.56 a	274.01 a	27.24 a
	Harvest 2	1.34 a	11.39 a	14.59 a	0.78 a	163.04 a	304.46 a	6.16 a	281.30 a	26.92 a
2020	Harvest 1	1.57 a	12.02 a	17.06 a	0.71 a	200.45 a	304.16 a	5.76 a	232.78 a	25.64 a
	Harvest 2	1.37 a	11.38 a	15.45 a	0.74 a	153.27 a	286.77 a	6.66 a	235.62 a	25.14 a
		‘Sunrise’
		Corvallis, Oregon
2019	Harvest 1	1.27 a	11.67 a	15.06 a	0.78 a	152.75 a	328.55 a	6.19 a	285.58 a	28.21 a
	Harvest 2	1.21 a	10.28 a	13.50 a	0.76 a	146.72 a	315.38 a	6.61 a	262.09 a	30.89 a
2020	Harvest 1	1.07 a	9.93 a	12.78 a	0.78 a	198.41 a	315.17 a	5.07 a	220.69 a	31.87 a
	Harvest 2	0.84 a	9.01 a	11.59 a	0.79 a	221.32 a	396.62 a	4.74 a	324.88 a	47.02 a
		Jackson Springs, NC
2019	Harvest 1	1.04 a	10.39 a	14.58 a	0.71 a	171.01 a	310.55 a	5.10 a	320.05 a	30.02 a
	Harvest 2	1.03 a	9.50 a	12.36 a	0.77 a	190.07 a	350.86 a	5.52 a	303.74 a	37.36 a
2020	Harvest 1	1.26 a	10.59 a	14.58 a	0.73 a	175.99 a	275.45 a	5.48 a	233.75 a	26.26 a
	Harvest 2	1.25 a	11.06 a	14.96 a	0.74 a	180.86 a	280.46 a	5.81 a	207.68 a	25.52 a

^† FW = fruit weight (g); PD = polar diameter (mm); ED = equatorial diameter (mm); RI = roundness index; FRM = firmness (g·mm⁻¹); APF = absolute positive force (g); DPF = distance at absolute positive force (mm); FT = force at target (g); FE = fruit elasticity. ^‡ Means followed by the same letter(s) within a year under subheading of environment, and genotype are not significantly different using the least squares means (LSMEANS) Tukey HSD multiple comparisons procedure (p < 0.05) between harvests. ^§ ‘Summit’ did not have a second harvest in Corvallis, 2020.

and

Table 7. Harvest means for each of the five cultivars within a year, environment, and genotype for fourteen phenotypic traits using model 2 in PROC GLIMMIX (SAS v9.4).

Year	Harvest	L* †	HUE	CRM	TSS	TA
		‘Echota’
		Corvallis, Oregon
2019	Harvest 1	27.48 b ‡	264.06 a	3.87 b	11.75 a	0.87 a
	Harvest 2	29.24 a	263.94 a	4.54 a	12.80 a	0.57 a
2020	Harvest 1	31.32 a	264.43 b	5.58 a	12.20 a	0.88 a
	Harvest 2	28.64 b	267.13 a	4.86 b	13.35 a	0.99 a
		Jackson Springs, NC
2019	Harvest 1	32.15 a	264.54 a	5.78 a	12.80 b	1.24 a
	Harvest 2	33.06 a	263.30 a	5.71 a	15.15 a	1.34 a
2020	Harvest 1	32.93 a	268.89 a	5.93 a	10.95 a	1.08 a
	Harvest 2	29.12 b	265.90 b	5.29 b	9.75 a	1.39 a
		‘Echota’
		Corvallis, Oregon
2019	Harvest 1	24.03 a	262.80 b	2.92 a	13.45 a	0.43 a
	Harvest 2	27.47 a	266.43 a	3.73 a	13.95 a	0.32 a
2020	Harvest 1	26.61 a	269.84 a	3.97 a	13.70 a	0.64 a
	Harvest 2	23.54 a	268.47 a	3.33 a	13.39 a	0.68 a
		Jackson Springs, NC
2019	Harvest 1	32.28 a	264.18 a	4.73 a	15.3 a	0.28 a
	Harvest 2	27.11 a	265.15 a	4.09 a	16.1 a	0.42 a
2020	Harvest 1	22.77 a	263.43 a	2.83 a	11.1 a	1.04 a
	Harvest 2	24.83 a	265.48 a	4.87 a	13.8 a	0.56 b
		‘Echota’
		Corvallis, Oregon
2019	Harvest 1	30.37 a	268.04 a	4.27 a	16.10 a	0.58 a
	Harvest 2	31.05 a	268.10 a	3.71 b	16.55 a	0.40 b
2020	Harvest 1	26.64 a	266.8 b	4.62 a	15.35 b	0.49 b
	Harvest 2	23.69 b	271.95 a	3.42 b	17.89 a	0.76 a
		Jackson Springs, NC
2019	Harvest 1	31.06 a	262.18 a	4.77 a	16.30 a	0.64 a
	Harvest 2	28.42 b	264.05 a	4.36 a	15.05 a	0.69 a
2020	Harvest 1	27.69 a	271.64 a	4.84 b	11.55 b	1.17 a
	Harvest 2	29.48 a	263.31 b	5.84 a	15.45 a	1.13 a
		‘Echota’
		Corvallis, Oregon
2019	Harvest 1	29.12 a	264.98 a	5.03 a	14.6 a	0.22 a
	Harvest 2	28.39 a	265.32 a	4.78 a	14.0 a	0.39 a
2020	Harvest 1	27.87	264.74	5.19	16.5	0.46
	Harvest 2 §
		Jackson Springs, NC
2019	Harvest 1	31.61 a	267.43 a	4.46 a	12.05 a	0.84 a
	Harvest 2	34.41 a	263.06 a	4.86 a	13.15 a	0.57 a
2020	Harvest 1	29.01 a	269.73 a	4.75 a	10.50 a	0.65 a
	Harvest 2	27.41 a	266.69 a	4.58 a	10.80 a	0.35 a
		‘Sunrise’
		Corvallis, Oregon
2019	Harvest 1	28.27 b	263.07 a	4.19 a	14.25 a	0.42 a
	Harvest 2	30.69 a	264.68 a	3.75 a	14.70 a	0.50 a
2020	Harvest 1	27.26 a	265.72 a	4.25 a	12.35 b	0.51 a
	Harvest 2	24.80 b	267.06 a	3.14 b	14.84 a	0.58 a
		Jackson Springs, NC
2019	Harvest 1	30.34 a	264.98 a	4.97 a	12.70 a	1.34 a
	Harvest 2	30.69 a	264.57 a	4.55 a	14.00 a	1.18 a
2020	Harvest 1	29.30 a	265.47 a	5.00 a	10.65 a	0.93 a
	Harvest 2	27.40 b	264.08 a	3.79 b	13.20 a	0.80 a

^† L* = luminescence; HUE = hue angle; CRM = chroma; TSS = total soluble solids (%); TA = titratable acidity (%). ^‡ Means followed by the same letter(s) within a year under subheading of environment, and genotype are not significantly different using the least squares means (LSMEANS) Tukey HSD multiple comparisons procedure (p < 0.05) between harvests. ^§ ‘Summit’ did not have a second harvest in Corvallis, 2020.

). Comparing between genotypic averages, ‘Reveille’ had the highest RI, FRM, APF, FT, FE, and TSS values (Table 2); however, within genotype analysis between environments, our research indicated that ‘Reveille’ had the most significant variability in FRM, APF, and FE (Table 6). Specifically, the Corvallis environment had increased measures in all fruit’s firmness-related traits compared with Jackson Springs (Figure 1D,E). Overall, ‘Reveille’ fruit’s firmness-related traits were greater within a respective environment than those observed in other genotypes. In Corvallis, ‘Reveille’ had significantly greater values than other genotypes in eight of the measured traits (RI, FRM, APF, DPF, FT, FE HUE, and TSS). In contrast, in Jackson Springs, ‘Reveille’ had significantly greater values than other genotypes in all of the above traits except for DPF and HUE. At the Jackson Springs site, ‘O’Neal’ had significantly lower values in nine of the fourteen traits (FRM, APF, DPF. FT, FE, L*, HUE, CRM, and TA). Among these traits, only L*, CRM, and TA were statistically significant compared to other genotypes in either environment (Figure 1F,H, Supplementary Figure S3). Compared to other genotypes in Corvallis, ‘Summit’ had significantly lower FRM, APF, DPF, FE, and TSS values and higher FW, PD, and ED. ‘Summit’ traits including, FW, PD, ED, FRM, APF, and FE were consistent with statistical comparisons between genotypes at the Jackson Springs location (Table 6 and Table 7).

Overall, FW, PD, ED, and TA increased from year 1 to year 2 in ‘O’Neal’ and ‘Reveille’ while FRM and HUE increased across years for all genotypes. ‘Summit’ lacked a second harvest in the second year of study in Corvallis; however, averages of FW, PD, ED, and TA increased between years of study (Table 6). All genotypes showed decreased measures of TSS across years (Figure 2D,G, Supplementary Figure S4). ‘Echota’ and ‘O’Neal’ had the most significant variation among genotypes across years; ‘Echota’ decreased in FW, RI, and TSS and ‘O’Neal’ decreased in RI, FRM, APF, DPF, FE, and HUE from year 1 to year 2 (Figure 2C–E, Supplementary Figure S4). Statistical analysis between genotypes showed that ‘Reveille’ and ‘Sunrise’ were more consistent across the years and showed consistency for phenotypic trait values within a genotype across environments (Table 6 and Table 7).

3.3. Correlation between Agronomic Traits

Individual fruit measures were averaged within each clonal replicate (per genotype, location, harvest, and year) for each studied trait that utilized individual fruit (FW, PD, ED, RI, FRM, APF, DPF, FT, and FE). Pairwise Pearson’s correlation between each of the fourteen studied agronomic traits was performed using phenotypic measures of L*, HUE, CRM, TSS, and TA and averages per clonal replicates for FW, PD, ED, RI, FRM, APF, DPF, FT, and FE (Figure 3). We identified five strong positive (>0.70) pairwise correlations. Strongly correlated traits included FW-PD (r = 0.87), FW-ED (r = 0.94), PD-ED (r = 0.83), APF-FT (r = 0.76), and APF-FE (r = 0.92). Measures of deflection (FRM), initial resistance to puncture (FT), and force required for fruit puncture (APF) are all fruit measures relating to the firmness and innate resistance to bruising. Additionally, FE was calculated by DPF, the depression of the probe on the fruit between FT and APF, and ED; therefore, DPF is related to both fruit volume and fruit firmness measures.

Further, we identified twenty-six moderately correlated phenotypic traits (|0.70|> |r| >|0.40|). Of the identified moderately correlated traits, ten were positively correlated, including RI-APF (r = 0.43), FRM-APF (r = 0.63), FRM-FT (r = 0.48), RI-FE (r = 0.40), FRM-FE (r = 0.53), FT-FE (r = 0.69), L*-CRM (r = 0.69), APF-TSS (r = 0.45), FE-TSS (r = 0.51), and CRM-TA (r = 0.46). The sixteen traits with a moderate negative moderate correlation included FW-RI (r = −0.47), ED-RI (r = −0.62), FW-FRM (r = −0.43), ED-FRM (r = −0.44), FW-APF (r = −0.50), PD-APF (r = −0.45), ED-APF (r = −0.58), ED-FT (r = −0.41), FW-FE (r = −0.63), PD-FE (r = −0.60), ED-FE (r = −0.67), FW-TSS (r = −0.48), PD-TSS (r = −0.42), ED-TSS (r = −0.50), RI-TA (r = −0.47), and TSS-TA (r = −0.43) (Figure 3).

4. Discussion

As blueberry production expands, it encompasses increased geographic and climatic diversity. Regions with compatible climates for both northern and southern highbush genotypes accommodate extended harvest windows. However, genotype response to diverse environments is largely uncharacterized. The main objective of the present study was to determine the extent of variation on fruit quality traits of five blueberry genotypes resulting from different sources of variation, including environment, year, harvest, and their interactions. Understanding the sources attributed to phenotypic plasticity informs that trait’s genetic control and stability [10]. This research evaluates sources of variation significant to phenotypic stability and genotypic performance in select environments. Additionally, our analysis over multiple environments, years, and harvests was previously uncharacterized in blueberry for this range of phenotypic traits.

The selected environments saw differences in precipitation and air and soil temperatures. Corvallis had decreased net precipitation, air, and soil temperature, and increased evapotranspiration compared to the Jackson Springs environment in both years of study (

Table 8. Climatic conditions including air and soil temperature, precipitation, and evapotranspiration during blueberry fruit maturation and harvest season in Corvallis, OR and Jackson Springs, NC in 2019 and 2020.

Environment, Year	Air Temperature (°C)
	March		April		May		June		July		August		Total
	Max	Min	Max	Min	Max	Min	Max	Min	Max	Min	Max	Min	Max	Min
Corvallis, 2019	34.16	23.56	34.16	23.56	39.37	24.54	43.17	26.43	44.91	29.01	46.74	29.48	41.67	26.60
Corvallis, 2020	35.25	21.63	35.25	21.63	38.32	25.39	40.85	28.06	46.02	28.92	47.46	28.69	41.58	26.54
Jackson Springs, 2019	33.41	22.58	40.77	29.30	46.78	35.26	47.68	37.10	50.49	38.85	48.93	37.97	46.93	35.70
Jackson Springs, 2020	36.04	26.40	40.06	27.74	41.44	30.52	46.51	36.36	50.53	39.07	48.47	39.03	45.40	34.54
	Evaporation (EVP) (cm) and Precipitation (PRECP) (cm)
	March		April		May		June		July		August		Total
	EVP	PRECP	EVP	PRECP	EVP	PRECP	EVP	PRECP	EVP	PRECP	EVP	PRECP	EVP	PRECP
Corvallis, 2019	0.00	6.25	6.90	15.80	15.07	4.67	22.36	1.85	23.67	1.68	23.64	0.29	18.32	4.86
Corvallis, 2020	0.00	6.71	14.58	4.50	14.63	7.68	19.08	4.67	27.05	0.00	23.11	0.32	19.68	3.43
Jackson Springs, 2019	8.36	10.79	11.52	20.24	15.80	6.17	15.21	10.33	16.23	9.55	13.71	9.36	14.48	11.13
Jackson Springs, 2020	6.88	6.08	11.37	8.92	11.30	18.23	14.34	9.16	15.80	9.77	12.34	8.65	13.02	10.94
	Soil Temperature (°C)
	March		April		May		June		July		August		Total
	Max	Min	Max	Min	Max	Min	Max	Min	Max	Min	Max	Min	Max	Min
Corvallis, 2019	28.19	23.05	33.79	28.89	39.72	32.37	44.59	36.41	45.93	37.99	47.35	39.91	42.28	35.11
Corvallis, 2020	29.85	25.09	34.41	28.93	38.19	33.12	40.83	35.50	45.27	38.08	46.67	39.13	41.07	34.95
Jackson Springs, 2019	29.23	34.50	42.24	31.72	50.74	39.07	51.59	40.86	54.43	42.74	52.42	41.56	50.28	39.19
Jackson Springs, 2020	36.44	28.35	42.70	31.31	45.27	34.34	51.58	39.91	55.44	42.87	52.56	42.26	49.51	38.14

). However, our results indicated that E, G × E, G × E × Y, and G × E × H were not significant sources of variation for the majority of the phenotypic traits in either model. Yearly climatic differences in rainfall and temperature events suggested a greater influence on phenotypic expression. Precipitation and evapotranspiration are transitory effects that fluctuate year to year. Environmental differences are more evident in soil and air temperature and their effect on fruit maturation times. Air and soil temperatures warmed more gradually, and fruit matured later at the Corvallis location compared with Jackson Springs. While similar to Corvallis regarding maximum air and soil temperatures in March, the Jackson Springs environment experiences temperature increases in April that are not observed until May at the Corvallis location. The warmer air and soil temperatures at the Jackson Springs site agree with previous studies, which found that increased daytime soil and air temperatures advanced reproductive and vegetative growth of blueberry plants grown in high tunnels. Advancing the timeline of reproductive growth ultimately accelerated the fruit ripening by five weeks compared to blueberries under field conditions [11]. Further, higher temperatures have been reported to induce pigment accumulation and accelerate ripening in blueberries [25].

Yearly phenotypic variation is likely resultant from the differential precipitation and temperatures between years of analysis as observed within each environment—increased precipitation results in larger fruit in blueberry [6]. Across environments, fruit size metrics were not significantly different; however, they varied in accordance with annual precipitation events within a location. Within the Corvallis environment, increased precipitation in 2019 from rainfall events resulted in larger fruit, as exemplified by the 14% increase in FW in 2019 compared to 2020, Where precipitation in the Jackson Springs environment decreased by 13% between 2019 and 2020, FW was affected by a 26% decrease. The soil type in Jackson Springs is sandy, and it was amended by pine bark to acidify the soil and increase organic matter. The plants were irrigated by a linear movement irrigation system every 7–10 days. However, the soil did not have much capacity to retain moisture even after incorporating pine bark. Further, PD, ED, RI, and DPF decreased in the Jackson Springs environment from 2019 to 2020 and increased in Corvallis. Fruit volume metrics are likely to increase due to precipitation events within a location [6]. Conversely, water stress conditions, including high-temperature effects and drought conditions, negatively affect fruit volume metrics [16,26,27]. Fruit volume metrics increased at the Jackson Springs and decreased at the Corvallis locations between the years of the study, which correspond to increased temperatures during harvest. While high-chill northern blueberry cultivars are negatively impacted at temperatures ≥ 30 °C, noted by a reduced photosynthetic rate [28], low-chill-adapted blueberry cultivars have a higher estimated photosynthetic heat tolerance with no significant electrolytic leakage or superoxide radical accumulation observed at temperatures < 40 °C [29]. While high-chill northern blueberry cultivars are negatively impacted at temperatures ≥ 30 °C, noted by a reduced photosynthetic rate [28], low-chill-adapted blueberry cultivars have a higher estimated photosynthetic heat tolerance than northern high-chill-adapted cultivars with no significant electrolytic leakage or superoxide radical accumulation observed as a result of heat stress until at temperatures ≥ 40 °C [29]. In our findings, the maximum temperatures during harvest (May–June at the Jackson Springs location and July–August at the Corvallis location) exceeded 40 °C at both locations in both years of study. Overall, fruit volume metrics were higher in years with a lower mean temperature during harvest. The decreased fruit volume metrics in the NHB cultivars ‘Echota’ and ‘Sunrise’ were related to increased annual temperatures between years for both environments; these observations may support decreased photosynthetic rates due to of heat stress. Coinciding with years of lower maximum air temperatures, ‘Sunrise’ and ‘Echota’ had increased fruit volume metrics in Corvallis in the first year and Jackson Springs in the second year of study. However, it is noteworthy that the increased fruit volume metrics at the Jackson Springs location were likely significantly influenced by the three-fold increase in precipitation in the second year of study. Further research examining the effects between NHB and optimal fruit production temperature would establish a baseline for comparisons for fruit volume metrics.

FW and ED had a moderate negative correlation with fruit firmness metrics, including FRM, APF, and FE (r > 0.40). We observed decreases of 17%, 6%, and 15% in traits, respectively, from 2019 to 2020 in Corvallis. Conversely, fruit volume metrics of FW, PD, ED, and DPF increased 14%, 9%, 10%, and 10%, respectively, from 2019 to 2020 in Corvallis. Previous reports indicate that decreased precipitation corresponds to firmer fruit [6]. While Corvallis plants are consistently irrigated twice a week, climate data reported minimal rainfall in July and August in both years of study and increased net evapotranspiration. Corresponding to increased evapotranspiration during harvest, FW, PD, and ED all decreased while FRM, APF, and FE all increased at Corvallis in 2020. At the Jackson Springs location, only APF and FE had increased values coinciding with increased evapotranspiration in 2019 compared with 2020. Contrasting with the expected increase in precipitation and reduced evapotranspiration in May and June 2020, FRM measurements were significantly higher in Jackson Springs. At the same time, FW also increased from 2019 to 2020, likely due to the three-fold increase in precipitation. Thus, both FW and FRM simultaneously increased in Jackson Springs between 2019 and 2020. FRM increased in ‘Echota’, ‘O’Neal’, and ‘Summit’, all of which increased FW. In contrast, ‘Reveille’ showed an increase in FW but a decrease in FRM. Only ‘Sunrise’ had a lower FW and higher FRM between years of the study in Jackson Springs. Contrary to environmental influences between the first and second years of study, ‘Echota’, ‘O’Neal’ and ‘Summit’ exhibited higher fruit volume-related metrics and lower FRM, APF, FT, and FE values. ‘Reveille’ reflected expectations of increased fruit volume metrics and decreased fruit firmness measurements at Jackson Springs from 2019 to 2020. The overall moderate negative correlation between FW-FRM and ED-FRM indicates that denser fruit with higher weight values and fruit with a larger equatorial diameter had lower fruit firmness. These correlations and an inverse relationship with precipitation and evapotranspiration suggest minimal temporal influence over these phenotypic traits.

While increased water uptake due to increased irrigation or precipitation in blueberry fruit positively correlates with fruit volume metrics, it negatively correlates with TSS [5]. The Jackson Springs location faced lower precipitation in May in 2019 and higher rainfall in June than 2020. More importantly, evapotranspiration was nearly 40% and 6% greater in May and June 2019, respectively, compared with 2020 at Jackson Springs. Blueberry plant water loss through evaporation reduced fruit volume metrics, including FW, PD, and ED, correlated with increased TSS values (r > 0.40). Corvallis experienced increased evaporation in 2020, which had decreased FW, PD, and ED and increased TSS.

Additionally, the increased air temperatures during harvest at each location may positively affect fruit TSS values. For example, at Jackson Springs, TSS values averaged across cultivars in 2019 were 21% greater than 2020, when May and June maximum air temperatures were 12.8% and 2.5% greater in 2019 than 2020. Comparatively, there was only a slight increase in both maximum air and soil temperature (<3% in July and August) and TSS (<2%) values in fruit collected from Corvallis in 2020. Corvallis blueberry harvest in July and August was hotter with cooler nighttime temperatures in 2020 compared with 2019; however, TSS showed a minimal increase. Contrary to our results, previous studies in blueberry report a negative correlation between decreasing light level or photosynthetic active radiation and canopy temperature with reduced fruit water content and fruit weight, while TSS increased [30].

In contrast, another study found a higher air and soil temperatures due to growing blueberry plants in high tunnels had no reported effect on TSS compared with control blueberry plants in the field [11]. These studies did not separate the light level and temperature treatment effects. Thus, conclusions may differ as there have not been studies that isolate the direct effect of air temperature and GDD without impacting light level and air movement on pre-harvest phenotypic characteristics in blueberry fruit. Similar to blueberry, red grapes are non-climacteric fruit that accumulates pigment during ripening. In grapes, high temperatures inhibit pigmentation and anthocyanin accumulation in select cultivars [31]. In contrast to Jackson Springs blueberry observations, reduced temperatures in grape, with a mean temperature of 20.3 °C, reflected increased TSS values compared to a warmer (37 °C/32 °C day/night) temperature treatment [32]. In a similar study, TSS was highest in interspecific wine grape hybrid held at a constant 25 °C for 15 to 37 days and lower at 20 °C and 30 °C, and TA was highest in fruit harvested from plants kept at 20 °C [33]. A study on tomatoes also found that increasing temperature from 26° C to 32 °C corresponded with increased TSS values [8].

Our results found a moderate negative correlation between TSS and TA. The relationship between decreased precipitation and fruit volume is positively correlated with TSS. However, we did not find any correlation between fruit volume metrics and TA. Interestingly while it has been shown that plants under decreased precipitation and increased water stress had higher TA and TSS [5,6], our results showed that there was a negative correlation between TSS and TA.

The HUE angle of blueberry defines the gradient of the blue color of the fruit on a color spectrum; values closer to 270° are considered darker blue, and values approaching 220° are lighter blue. CRM, on the other hand, indicates the level of color intensity or saturation. From 2019 to 2020, CRM values decreased in Jackson Springs and increased in Corvallis, coinciding with warmer harvest temperatures at each location. Chromaticity values have been related to pigment accumulation in fruit [34]. Previous results found that post-harvest fruit and blueberry puree had higher chroma values at lower temperatures, agreeing with our findings at the Jackson Springs location; however, these results were not significant between 4 °C and room temperature, 75 °C, and 100 °C in a study of three NHB blueberry cultivars. Overall, L* and HUE decreased and increased, respectively, between years in our study and do not correspond to temperature or rainfall events in either environment. This observation was consistent across G × Y except for ‘Echota’, for which L* values slightly increased between years of study in Corvallis. Fruit color is a factor of fruit ripeness and maturation, which can be altered by harvest time and anthocyanin accumulation. As blueberry fruit ripens, anthocyanin content increases, turning from green to blue with intermediate red or blush hues. Further, the fruits of most blueberry cultivars produce a wax coating or bloom, which can preclude color determination. Fruit bloom affects L* in the light-dark color perception of the fruit and is impacted by fruit handling and weather, diminishing the bloom coating. As such, we would expect L* values to have increased variability.

As a whole, blueberry fruit ripening is indicated by fruit pigment accumulation for one month. Plant fruit ripening and the blueberry reproductive cycle are based upon genotype and environmental factors controlling budbreak. Fruit derived from the same plant but collected at a later harvest due to a later fruit maturation window will have accumulated more GDD and increased opportunities for rainfall events and temperature increases. As such, we would expect to see differences between harvests. Our results identified significant differences in both ‘O’Neal’ and ‘Reveille’ regarding fruit volume metrics between harvests, as well as differences in L*, HUE, CRM, and TSS in ‘Echota’, ‘Reveille’, and ‘Sunrise’ between harvests.

It was intriguing that we did not find any consistent effect between harvests within a genotype across environments or years for any phenotypic trait. However, we did observe that TSS increased between harvests for most genotypes on an environment and yearly basis. While previous findings reported blueberry fruit decreasing in CRM values and L* during ripening, our results did not find any consistent trend in either CRM or L* [25]. However, a major difference between the previous study [25] and ours is the evaluation of blueberry fruit ripening stages covering pigmentation. Our study looks at similarly ripe fruit and examines fruit variation between successive harvests. Variation in fruit measures within a plant between successive harvests signify differences in GDD and temperature events between harvest timepoints. As fruit matured and accumulated anthocyanins, Spinardi et al. (2019) found that fruit shifted in hue angle from 290° (violet) to 210° (blue) during ripening. In contrast, our findings did not see a shift in HUE between harvests, further supporting the uniform ripeness of fruit between harvests.

Intriguingly, we did not find any differences in fruit firmness metrics associated with either harvest timepoint. A previous study examined the metrics between blueberry fruit ripeness and canopy position found that fruit six days post-optimal-ripeness were significantly less firm than optimally-ripe fruit [15]. While our study only evaluated ripe fruit, it is important to note that increased GDD did not impact fruit firmness metrics. While average harvest values of FRM and FT decreased between harvests, the difference was not significant. Further, fruit firmness metrics of APF and FE increased. Additional interactions with genotype, environment or year were not a significant source of variation with harvest effects for any measured fruit firmness metrics.

Genotype is an expected significant source of variation. Accordingly, all studied phenotypic traits displayed significant differences between genotypes apart from DPF, and HUE. Overall, ‘Reveille’ performed highest in fruit firmness metrics (FRM, APF, FT, and FE) across genotypes within either environment or year of study. ‘Summit’ consistently had higher fruit volume metric averages (FW, PD, and ED) across the evaluated genotypes in both the environments and years of the study. Similarly, ‘Echota’ had the second-highest values for fruit volume metrics, and, as with ‘Summit’, lower fruit firmness metrics were observed, including FRM, APF, PF, FT, and FE. ‘Echota’, followed by ‘Summit’, had the highest L* and CRM values rating across genotypes. Higher L* values signify lighter fruit skin coloration as a result of pigmentation or genotypic bloom production. The current market preference in blueberry visual appeal coincides with higher L* values and light blue blueberry fruit in addition to larger fruit size characteristics [35] displayed by both ‘Echota’ and ‘Summit’. ‘O’Neal’ was the least acidic genotype in this study and had a rounder fruit similar to ‘Reveille’ yet had the lowest L* and CRM values. ‘O’Neal’ was average among the majority of studied traits, representing both medium volume and firmness. Similarly, ‘Sunrise’ represented the median value for nearly every phenotypic trait; this observation was consistent in both environments for all phenotypic traits except for DPF, wherein ‘Sunrise’ had the highest DPF value across genotypes in Corvallis.

5. Conclusions

The results of this study indicated the presence of significant interactions among genotypes, years, environments, and harvests for all measured fruit quality parameters. Genotype is a significant source of variation for the majority of phenotypic traits. ‘Reveille’ and ‘O’Neal’ phenotypic stability were consistent across environments and years; moreover, ‘Summit’ phenotypic data across years, locations, and harvests suggested similar stability. Clonal plant replicates within genotype and environment, and individual fruit measures were the most significant source of variability. This research established thirty-one moderate-to-strong correlations between multiple phenotypic traits, including fruit volume, firmness, color, soluble solids, and acidity. Air and soil temperature, precipitation, and evaporation were variable between years and suggested influence on phenotypic measures. Overall, these data provide a valuable foundation for breeding blueberries in different target locations and elucidating climatic effects on fruit characteristics.