Assessing Yield and Yield Stability of Hevea Clones in the Southern and Central Regions of Malaysia

Increased volatility in global rubber prices has led to declining Malaysian rubber production and smallholder income. Identifying rubber clones that can produce a consistently high yield in various environments is one of the potential measures to alleviate the impact of price fluctuations and improve smallholder livelihoods. In this study, we assessed rubber yields and yield stability of 37 rubber tree clones at two major production regions in Malaysia: Kota Tinggi (Southern region) and Sungai Buloh (Central region). In addition, we assessed relationships between climate data with rubber yields. Rubber yield and stability differed widely by clone, but showed relatively consistent trends across regions. Clones RRIM 2007, PB 260, and RRIM 2012 were high yielding in both regions and had high stability. Conversely, clone RRII 308 had the lowest mean yield across both regions and low stability. Mean annual yields showed a negative relationship with rising minimum temperatures, declining by ~3 g per tapping per tree (g t−1 t−1) per 1 ◦C rise in minimum temperature. Our findings highlight clones capable of achieving stable high yields. This information may be useful for breeders and agronomists in identifying germplasm and traits for further development. Further, this information can be used to assist clone recommendations to smallholders in these regions to mitigate the dual impacts of climate change and market volatility.


Introduction
Malaysia is the world's fifth-largest rubber producer, and together with Thailand, Indonesia, Vietnam, China, and India, it contributes more than 90% of global rubber production [1,2]. Rubber (Hevea brasiliensis) is the second largest commodity produced in Malaysia after oil palm and contributes 4.7% to the national gross domestic product [3]. In addition to supporting the national income, rubber plantations also play a substantial role in providing income for smallholders. This is because smallholders dominate more than 90% of Malaysia's rubber production [2]. The primary concern of smallholders is global rubber price fluctuations, which have trended downward in recent years [2,4]. Declining revenue drives smallholders to abandon farm management, resulting in a decline in national production [5,6]. One approach to address this issue is to ensure smallholder rubber producers are farming high yielding rubber clones with concomitant high yield stability. This would increase the consistency and productivity of rubber plantations and help secure smallholder income.
The average rubber yield produced by Malaysian smallholders is~1400 kg ha −1 yr −1 [2], which is far below the theoretically attainable yield of 7000-12,000 kg ha −1 yr −1 [7]. This yield gap can be driven by several factors, including crop genetics [8,9], environmental conditions [10,11], and agronomic practices [12,13]. One approach to improve smallholder productivity is to develop new varieties that are well-adapted to their environment and can produce stable and high yields. Rubber research The breeding program for the rubber clones was conducted at two MRB research stations, located at Kota Tinggi (herein referred to as the Southern region) and Sungai Buloh (herein referred to as the Central region) of Peninsular Malaysia (Figure 1). These research stations allocated clonal trial fields for MRB researchers to compare annual clonal growth and yield performance. In each region, a randomised complete block design (RCBD) was applied, which has three replicates of 30 trees per plot with 5 × 4 m spacing between rows and trees. Trees were tapped using the S/2 D/3 tapping system, which means tapping at a half spiral diameter of the trunk at three-day intervals, resulting in ~120 tapping days/year. The yield, in the form of latex, was collected in a plastic cup attached to each tree. Yield was expressed as grams of dry rubber per tapping per tree (g t −1 t −1 ). This is commonly used to assess clonal yield and yield stability performance, which is the focus of this study. However, g t −1 t −1 can also be used to indicate farm performance by converting to kg ha −1 yr −1 using Equation 1: where planting density is the number of trees per hectare and tapping days is the number of tappings per year. Throughout this paper, we use g t −1 t −1 because the planting density and tapping days in MRB research trials are constant across plots. Additionally, g t −1 t −1 is not affected by tapping days, which can be influenced by rainfall intensity, price fluctuations, and labour availability. We provide the equation in case readers wish to convert our data to kg ha −1 yr −1 values. Both Southern and Central regions have similar annual temperatures, while the Southern region is marked by greater mean The yield, in the form of latex, was collected in a plastic cup attached to each tree. Yield was expressed as grams of dry rubber per tapping per tree (g t −1 t −1 ). This is commonly used to assess clonal yield and yield stability performance, which is the focus of this study. However, g t −1 t −1 can also be used to indicate farm performance by converting to kg ha −1 yr −1 using Equation (1): where planting density is the number of trees per hectare and tapping days is the number of tappings per year. Throughout this paper, we use g t −1 t −1 because the planting density and tapping days in MRB research trials are constant across plots. Additionally, g t −1 t −1 is not affected by tapping days, which can be influenced by rainfall intensity, price fluctuations, and labour availability. We provide the equation in case readers wish to convert our data to kg ha −1 yr −1 values. Both Southern and Central regions have similar annual temperatures, while the Southern region is marked by greater mean annual rainfall ( Table 1). The Southern region has soil class III, while the Central region has soil class IV. Within Malaysia, soils have been split into five classes based on their suitability for growing rubber [47]; this is based on seven soil parameters including depth, slope, texture, rock percentage, drainage, pH, and nutrient level. Soil class I is considered most suitable for growing rubber, while class V is the least suitable.

Data Collection
Data for rubber yield and weather were compiled from July 2003 to June 2010 (Southern region; seven years) and July 2003 to June 2007 (Central region; four years). The MRB recorded yield on a monthly basis as an average of the 30 trees per clone per plot (Department of Crop Protection and Improvement). Daily climate data for the periods mentioned were obtained from the NASA Langley Research Center (LaRC) POWER Project (funded through the NASA Earth Science/Applied Science Program). The NASA POWER Project data is designed to provide agroclimatological data for input into crop models. NASA POWER Project climate data has been found to compare favourably with local weather station data [48,49] and several studies have reported close agreement between crop yields estimated from NASA POWER climate data and observed yields in the field [49,50]. Climate data for the relevant time periods was collected for both locations and consisted of minimum, maximum and mean temperatures, and rainfall and relative humidity (RH). Vapour pressure deficit (VPD) was subsequently calculated using Equation (2) [51]: where VP sat is the saturated vapour pressure deficit. Daily VP sat values were calculated based on daily air temperature values using Equation (3) [51], where T is the daily air temperature: Annual climate data were calculated as within-year averages (July to June) based on daily data. While local weather station data would be preferable to the NASA POWER data, the nearest weather stations were determined to be not representative of the study sites as they are located >80 km from the site and have different altitudes and annual rainfall.

Yield Stability Analysis
Several metrics have been developed to assess crop yield stability, including regression coefficients [52], coefficient of variation (CV) [53], and Taylor's Power Law Residuals (POLAR) [54]. The appropriate metric to use has been shown to be dependent on the structure of the data [54]. For this reason, we calculated all three metrics to determine the most appropriate one for our dataset.
The regression slope approach was conducted through the development of linear mixed effects models (LMEs) using the 'nlme' package of R software [55]. LME are an extension of simple linear models that allow both fixed and random effects and are particularly useful when there is nonindependence in the data, such as arises from a hierarchical structure [56]. For example, when there are multiple levels in a structure, such as yield performance of multiple clones in the same region, the variability in the outcome can be interpreted as being either within-clones or between-clones. We compared annual population rubber yields (defined as Environmental Index, EI) for each of the two locations over 2003-2010 (Southern region) and 2003-2007 (Central region) in relation to variation in individual clonal annual yields over the same period. Throughout our LME analyses, the most parsimonious model was selected based on the Akaike Information Criteron (AIC) and validated by plotting model residuals.
Each model contained a clone-level random intercept and slope structure, and the EI was fitted as the fixed effect. The equation was represented by the following model structure in R: clonal mean rubber yield~EI, random =~EI | clone. Individual clonal minimum yield potential (MYP) values were calculated from LME predictions of rubber yield at the minimum EI for each location. MYP is the minimum yield a clone produced under the observed worst environmental conditions (for each site this was determined to be the year with the lowest yield). High MYP is important because it indicates clones that provide the greatest relative yields and therefore income under poor environmental conditions. The regression slope for each clone was obtained as the best linear unbiased predictors (BLUPs) from the LME. The value of the coefficient defines yield stability, where <1 is stable, 1 is the average population stability (i.e., matches 1:1 with the EI), and >1 is less stable and thus more sensitive to environmental fluctuations [57]. The CV for each clone within a location was calculated as the quotient of rubber yield standard deviation over mean rubber yield.
POLAR was calculated as the residuals of a linear relationship of log variance and log yield. POLAR is a recently developed yield stability metric based on Taylor's Power Law and has been presented as an alternative to the regression slope and CV approaches, which are considered to have several drawbacks [54]. For the regression slope approach, if variances scale with environmental means, the confidence interval for yields of individual varieties becomes broader, and the regression slope of residuals from the environmental slope increases as the environmental mean increases. On the other hand, the CV can either increase or decrease with an increasing yield based on the slope values (b > 2 for the former, and b < 2 for the latter), or constant if b = 2. These arguments impose that slope and CV may be strongly correlated with yield, and therefore may not be reliable and independent yield stability metrics. However, the drawbacks of POLAR include lower visibility of relationships to yield when comparing treatments or genotypes with similar means and requires many measurements of yield to observe the relationships [54].
We compared the suitability of these three stability metrics for assessing rubber clone yield stability by conducting a Pearson correlation analysis of each metric with yield. The stability metric with the lowest correlation with yield was chosen as the most suitable one. We then plotted our chosen yield stability metric against clonal mean yield and superimposed a quadrant centered on the median values of clonal mean yield and our yield stability metric ( Figure 2). Clones that are high yielding with high yield stability are located in the bottom right quadrant ( Figure 2). Hence, use of these clones would help ensure smallholder productivity is both greater and more stable.

Rubber Yields and Climatic Data
In this analysis, we excluded yield data from the first year of harvesting as it is typically associated with lower yields [46] that can confound interpretation of the impact of climate on yields. This process reduced our data sets to six years (Southern region) and three years (Central region); because of the limited data set for the Central region, only data from the Southern region was analysed. The relationship between rubber yields and climate variables was analysed using LME. The average annual value for each weather variable, covering maximum and minimum temperatures, precipitation, humidity and VPD, was used in relation to yield.

Mean Yields and Yield Stability
Clonal annual yield averages (EI and individual clones) were ranked from the 'worst' to 'best' performing years to quantify clonal yield stability [52] (Figure 3). Figures 3A (Southern region) and 3B (Central region) show the performance of each clone against the EI compared to the population mean. The overall average yields were 51 g t −1 t −1 and 54 g t −1 t −1 for the Southern and Central regions, respectively. The results show that mean yields at both locations were similar, but that rubber clones varied widely in terms of productivity and stability. Interestingly, the top three clones in the Central region outperformed all clones in the Southern region, despite the Central region having a lower soil quality rating.

Rubber Yields and Climatic Data
In this analysis, we excluded yield data from the first year of harvesting as it is typically associated with lower yields [46] that can confound interpretation of the impact of climate on yields. This process reduced our data sets to six years (Southern region) and three years (Central region); because of the limited data set for the Central region, only data from the Southern region was analysed. The relationship between rubber yields and climate variables was analysed using LME. The average annual value for each weather variable, covering maximum and minimum temperatures, precipitation, humidity and VPD, was used in relation to yield.

Mean Yields and Yield Stability
Clonal annual yield averages (EI and individual clones) were ranked from the 'worst' to 'best' performing years to quantify clonal yield stability [52] (Figure 3). Figure 3A (Southern region) and Figure 3B (Central region) show the performance of each clone against the EI compared to the population mean. The overall average yields were 51 g t −1 t −1 and 54 g t −1 t −1 for the Southern and Central regions, respectively. The results show that mean yields at both locations were similar, but that rubber clones varied widely in terms of productivity and stability. Interestingly, the top three clones in the Central region outperformed all clones in the Southern region, despite the Central region having a lower soil quality rating. MYP and mean yields were positively correlated (Southern region: r = 0.99, p < 0.001; Central region: r = 0.99, p < 0.001; Figure S1A, S1B; supplementary material), demonstrating that clones that performed better on average were also able to maintain higher yields under the worst environmental conditions observed during the study. Consequently, it may be preferable to plant clones with high MYP to secure a minimum level of productivity under suboptimal growing conditions. In both regions, LME regression slopes showed a strong positive relationship with mean yield ( Figure 4A-C (Southern region), 4D-F (Central region)), demonstrating that LME regression slopes are not independent of mean yield. This violates the compatibility of this metric to measure yield stability for our dataset [54].
In contrast, CV and POLAR showed no correlation with mean yields (for CV, Southern region: r = −0.22, p = 0.19; Central region: r = −0.09, p = 0.6; and for POLAR, Southern region: r = 0.14, p = 0.4; Central region: r = 0.14, p = 0.4; Figure 4). This demonstrates that CV and POLAR are independent of mean yield, indicating they are robust yield stability metrics for our data [54]. The more stable clones have both a lower CV and POLAR, which means there is less variation of yield across the harvesting period of that clone. In our subsequent analyses to identify stable and high yielding clones, we have chosen to present our results using POLAR due to its better presentation of independence from mean yield compared with CV and LME slopes (Figure 4). MYP and mean yields were positively correlated (Southern region: r = 0.99, p < 0.001; Central region: r = 0.99, p < 0.001; Figure S1A,B; Supplementary material), demonstrating that clones that performed better on average were also able to maintain higher yields under the worst environmental conditions observed during the study. Consequently, it may be preferable to plant clones with high MYP to secure a minimum level of productivity under suboptimal growing conditions. In both regions, LME regression slopes showed a strong positive relationship with mean yield ( Figure 4A-C (Southern region), Figure 4D-F (Central region)), demonstrating that LME regression slopes are not independent of mean yield. This violates the compatibility of this metric to measure yield stability for our dataset [54].
In contrast, CV and POLAR showed no correlation with mean yields (for CV, Southern region: r = −0.22, p = 0.19; Central region: r = −0.09, p = 0.6; and for POLAR, Southern region: r = 0.14, p = 0.4; Central region: r = 0.14, p = 0.4; Figure 4). This demonstrates that CV and POLAR are independent of mean yield, indicating they are robust yield stability metrics for our data [54]. The more stable clones have both a lower CV and POLAR, which means there is less variation of yield across the harvesting period of that clone. In our subsequent analyses to identify stable and high yielding clones, we have chosen to present our results using POLAR due to its better presentation of independence from mean yield compared with CV and LME slopes (Figure 4).
The categorisation of each clone into yield stability quadrants is shown in Figure 5A (Southern region) and Figure 5B (Central region). In this paper, we only present the results of the quadrants for mean yields due to the high correlation identified between the mean yield and MYP. The clonal mean yield and stability ranking (derived from the bottom-right quadrants) for the Southern and Central region is presented in Table S2; Supplementary material.  The categorisation of each clone into yield stability quadrants is shown in Figure 5A (Southern region) and Figure 5B (Central region). In this paper, we only present the results of the quadrants for mean yields due to the high correlation identified between the mean yield and MYP. The clonal mean yield and stability ranking (derived from the bottom-right quadrants) for the Southern and Central region is presented in Table S2; supplementary material.  Clones that have higher yields (>50 g t −1 t −1 ) and stability than the population median in both the Southern and Central regions include RRIM 2007 and PB 260 ( Figure 5). This finding is comparable with earlier research, which investigated yield performance based on three harvesting years in a large-scale clone trial [46]; however, the research did not consider the impact of different regions on the clonal yield performance. For RRIM 2007, our findings are consistent with the [46] trial, where it produced the highest yield (an average of 74 g t −1 t −1 ) compared to other clones. RRIM 2007 was reported to possess excellent secondary traits, including vigour at opening tapping panel, disease resistance, and yield potential [46]. In contrast, PB 260 was low yielding (27 g t −1 t −1 ) in the [46] trial. PB 260 was found to have poor resistance toward tapping panel dryness and virgin bark at opening tapping panel, and moderate resistance toward disease and wind damage, which may have contributed to the low yield [46]. Our results indicate that not only do both these clones produce high yields in the Central and Southern regions, their yields were also stable across the harvesting period. Additional high performing clones (stable and high yielding) in the Southern region were RRIM 928, RRIM 2001 and RRIM 2009, which yielded higher than the population median yields (>50 g t −1 t −1 ). These results are similar to what was reported from the large-scale clone trial [46]. These clones were reported to have excellent vigour and virgin bark at opening tapping panel, and resistance toward tapping panel dryness and wind damage [46]. Additional high performing clones identified in this study in the Central region were PB 372 and PB 374 (mean yields for both are~73 g t −1 t −1 and~81 g t −1 t −1 , respectively); the yield performance of these two clones has not been evaluated in previous studies.
In both regions, clones RRII 308, RRII 109, and RRIM 2014 showed low mean yields (17 g t −1 t −1 , 22 g t −1 t −1 , and 40 g t −1 t −1 at Southern; and 22 g t −1 t −1 , 24 g t −1 t −1 , and 39 g t −1 t −1 at Central) and low stability, and hence are not recommended for smallholders in these regions. This result is consistent with the findings presented by [58], where significantly higher CV (lower stability) and lower mean yields were observed in clone RRII 308 compared to clone RRIM 600 (Malaysian clone). However, [58] also revealed that despite RRII 308 having low mean yields, it produced competitive yields during winter months, indicating the importance of accounting for seasonal variation in clone yield performance. In the large-scale clone trial conducted by [46], RRIM 2014 had low mean yield (39 g t −1 t −1 ) [46], which is consistent with our finding. RRII 109, on the other hand, has not been studied previously in these regions, but based on this study appears unsuitable for planting. These findings highlight the role of cultivar origin in impacting clonal yield performance (e.g., clones from India generally performed poorly in the two regions in Malaysia in this study) and thus the importance of local clone breeding to ensure good yield performance.

Rubber Yields and Climate
Climatic factors have significant impact on rubber yield. For example, temperature and relative humidity influence tree stomata regulation, which affects latex flow [34,59,60], and rainfall intensity affects smallholders tapping activity [61]. In our study, the fixed effects are represented by the weather variables of minimum, maximum and mean temperatures, rainfall, relative humidity and VPD, while the fitted random effect was clone. The best-supported model indicates minimum temperature as the only significant variable and is negatively related to clonal mean yield ( Figure 6).
Although three of the 37 clones showed increases in yields as minimum annual temperatures increased, the other 34 clones showed a negative relationship (r = −0.97, p < 0.001). Across all clones, an increase of 1 • C minimum temperature reduced yield by~3 g t −1 t −1 (Figure 6). We found that the highest yield occurred at 23 • C and decreased as the minimum temperature increased, with the lowest yield at 26 • C. This finding is consistent with the literature, where rising minimum temperatures have been found to have a negative relationship with rubber yield [11,33]. While daytime temperatures of 27-33 • C are optimal for photosynthesis [62], temperatures above 22 • C are unfavourable for latex flow due to the increased evapotranspiration rates [63]. Further, decreased yield with higher minimum temperature suggests that higher night-time respiration reduces moisture availability within the tree, resulting in lower yield [34,64]. This finding highlights the potential impact of climate change on rubber yield. Current climate models indicate a continuous increase in minimum temperatures [65][66][67], which can lead to greater declines in latex yields due to an increase in tree respiration [68]. In addition to moisture loss, increasing respiration causes reductions in tree carbon-use efficiency, i.e., it increases losses of photosynthetically fixed carbon as CO 2 at the expense of conversion to proteins, carbohydrates and other compounds that comprise latex [61]. However, our analysis is based on NASA POWER data, which has a 0.5 • latitude, longitude grid cell spatial resolution. Consequently, while our analysis points to a significant negative influence of rising minimum temperatures on rubber tree productivity, local long-term climate data is needed to confirm this trend.
Agronomy 2020, 10, x 2 of 15 temperatures [65][66][67], which can lead to greater declines in latex yields due to an increase in tree respiration [68]. In addition to moisture loss, increasing respiration causes reductions in tree carbonuse efficiency, i.e., it increases losses of photosynthetically fixed carbon as CO2 at the expense of conversion to proteins, carbohydrates and other compounds that comprise latex [61]. However, our analysis is based on NASA POWER data, which has a 0.5° latitude, longitude grid cell spatial resolution. Consequently, while our analysis points to a significant negative influence of rising minimum temperatures on rubber tree productivity, local long-term climate data is needed to confirm this trend. Previous studies have found an inverse relationship between latex yields and VPD through changes to turgor pressure in tree laticifers [7,69,70]. High VPD also causes tree trunk shrinkage [71]. Thus, the absence of a relationship between yield and VPD in our study (r = −0.002, p = 0.7) was surprising. A previous study [71] reported that tree trunk shrinkage was only affected when VPD is >1 kPa. In our dataset, the average VPD value was 0.25 kPa, which may explain why no relationship was observed.

Conclusions
In this study, yield stability analysis of 37 Hevea clones in two different environmental conditions represented by the Southern and Central regions in Malaysia was conducted. We used three yield stability metrics to assess the clonal yield stability; we found that POLAR was the most suitable metric due to its lower correlation with mean yield and better presentation of independence from mean yield compared with CV and LME slopes. Clone RRIM 2007 was identified as the most stable and highest yielding with the greatest MYP in both the Southern and Central regions. The mean yields of this clone at the Southern and Central regions were 78 g t −1 t −1 and 75 g t −1 t −1 , respectively. These yields are substantially higher than the current mean yields of each region (27 g t −1 t −1 at Southern and 28 g Previous studies have found an inverse relationship between latex yields and VPD through changes to turgor pressure in tree laticifers [7,69,70]. High VPD also causes tree trunk shrinkage [71]. Thus, the absence of a relationship between yield and VPD in our study (r = −0.002, p = 0.7) was surprising. A previous study [71] reported that tree trunk shrinkage was only affected when VPD is >1 kPa. In our dataset, the average VPD value was 0.25 kPa, which may explain why no relationship was observed.

Conclusions
In this study, yield stability analysis of 37 Hevea clones in two different environmental conditions represented by the Southern and Central regions in Malaysia was conducted. We used three yield stability metrics to assess the clonal yield stability; we found that POLAR was the most suitable metric due to its lower correlation with mean yield and better presentation of independence from mean yield compared with CV and LME slopes. Clone RRIM 2007 was identified as the most stable and highest yielding with the greatest MYP in both the Southern and Central regions. The mean yields of this clone at the Southern and Central regions were 78 g t −1 t −1 and 75 g t −1 t −1 , respectively. These yields are substantially higher than the current mean yields of each region (27 g t −1 t −1 at Southern and 28 g t −1 t −1 at Central [2]). Thus, greater planting of RRIM 2007 in these regions may help improve smallholder productivity and income. The importance of this is underscored by the increased volatility in global rubber prices, which increases risk for smallholders. By planting the recommended clones, smallholders are more likely to ensure high yield during poor environmental conditions and maintain higher yield across the harvesting period, thereby minimising the impact of price fluctuations on their revenues. This quadrant-based analysis can also be used in other regions to identify high yielding and stable clones.
We also found that an increase of 1 • C minimum temperature was associated with a reduction in latex yields of~3 g t −1 t −1 . Globally, minimum temperatures are projected to continue to increase, based on general circulation models. This further highlights the importance of identifying clones that can perform well under poor environmental conditions and facilitating policy on planting those identified clones to minimise the impact of climate change on rubber production.
Supplementary Materials: The following are available online at http://www.mdpi.com/2073-4395/10/5/643/s1, Figure S1: Relationship between minimum yield potential (g t −1 t −1 ) and mean yields (g t −1 t −1 ), Table S1: Rubber tree clones examined in this study, Table S2: Clonal mean yield and stability ranking derived from the optimal quadrant for the Southern and Central region.