AFLA-PISTACHIO: Development of a Mechanistic Model to Predict the Aflatoxin Contamination of Pistachio Nuts

In recent years, very many incidences of contamination with aflatoxin B1 (AFB1) in pistachio nuts have been reported as a major global problem for the crop. In Europe, legislation is in force and 12 μg/kg of AFB1 is the maximum limit set for pistachios to be subjected to physical treatment before human consumption. The goal of the current study was to develop a mechanistic, weather-driven model to predict Aspergillus flavus growth and the AFB1 contamination of pistachios on a daily basis from nut setting until harvest. The planned steps were to: (i) build a phenology model to predict the pistachio growth stages, (ii) develop a prototype model named AFLA-pistachio (model transfer from AFLA-maize), (iii) collect the meteorological and AFB1 contamination data from pistachio orchards, (iv) run the model and elaborate a probability function to estimate the likelihood of overcoming the legal limit, and (v) manage a preliminary validation. The internal validation of AFLA-pistachio indicated that 75% of the predictions were correct. In the external validation with an independent three-year dataset, 95.6% of the samples were correctly predicted. According to the results, AFLA-pistachio seems to be a reliable tool to follow the dynamic of AFB1 contamination risk throughout the pistachio growing season.


Introduction
Pistachio (Pistacia vera L.) is a dioecious wind-pollinated tree cultivated for its highly nutrient nuts in many regions of the world under different climate conditions, but especially in subtropical and temperate zones [1]. Warm summer and mild winter conditions are preferable to achieve optimum yields and a high-quality product. It is well known that in order to attain high yields, chilling accumulation is essential during winter months, specifically from early October until late March [2,3], and hot conditions from April until harvest [4]. Several models exist to quantify the chilling accumulation, such as the Dynamic Model (in Chill Portions (CP)), which was originally developed for warm winters. It assumes that winter chill results from a process that takes place in two great impact on farmers' incomes. The significance of the economic impact of aflatoxin occurrence in agricultural products has been highlighted by a social network analysis of commercial exports and global trade patterns [23].
Good management practices at the postharvest level are clearly defined and they can guarantee no aflatoxin increase if properly applied. Therefore, during recent years, research has been focused on the prevention of aflatoxin contamination in pistachio fields using biological or chemical plant protection products [24][25][26][27]; nevertheless, no predictive models are available to support farmers' decisions in pistachio crop management. Predictive models have widely been used in Integrated Pest Management (IPM) systems of plant diseases in the frame of sustainable crop production, and several examples are available for mycotoxin management, such as the ochratoxin A (OTA) contamination of grapes [28], the fumonisin contamination of maize [29], the aflatoxin contamination of peanuts [30], and the deoxynivalenol contamination of winter wheat [31][32][33]. In the context of the present work, the most interesting model from the literature is AFLA-maize, a mechanistic model to predict A. flavus growth and AFB 1 contamination risk in maize [34,35]. It is based on two sub-models, one accounting for the host crop phenology and the other for the A. flavus infection cycle. The rationale behind the present study was to transfer the available knowledge, intended as the predictive maize model of A. flavus and aflatoxins, to a different pathosystem, A. flavus-pistachio. In particular, the objectives of this study were to (i) develop a predictive model for pistachio phenology in order to be combined with the sub-model of A. flavus infection cycle; (ii) define the susceptible pistachio growth stage for A. flavus infection; (iii) develop and validate the resulting predictive model with the aflatoxin data collected from Greek pistachio orchards.

Collection of Meteorological Data
The collected meteorological data were: the daily mean temperature (T, • C), maximum daily temperature (T max , • C), minimum daily temperature (T min , • C), daily relative humidity (RH, %), and daily precipitation (R, mm) from the 1st of January until the 31st of December for the years 2014-2019. The collected meteorological data were used to compute the pistachio growth stages (2014-2016) and to develop and validate the predictive model, AFLA-pistachio (2017-2019).

Computation of Pistachio Growth Stages
The island of Aegina is a warm and dry area where the daily minimum temperature during winter is above 7 • C, so there is no chilling accumulation (www.meteo.gr). As a result, the CH calculation for pistachio trees using the Crossa-Raynaud model was proved to be not applicable for the region of Aegina island due to high winter temperatures.
The GDD model was applied for the determination of the growth stages of the crop. The GDD based on the field observations for 2014, 2015, and 2016 were calculated for all the crucial growth stages of the crop (Table 1). For all the three years, the computation of GDD indicated that flowering occurs after at least 579 GDD and varies up to 836 GDD, with a mean daily temperature (T) in this period of 17.9 • C (mean of three years); 2015 had the lowest temperatures, while 2014 and 2016 were warmer years. Consequently, in 2015 all the growth stages occurred later than in 2014 and 2016. Harvesting takes place around 3000 GDD for the current location investigated, specifically 3013 GDD, 2845 GDD, and 3167 GDD for 2014, 2015, and 2016, respectively. The existence of small period windows leads to the result that the differentiation between the years studied occurs within a small frame.

Field Sampling and Aflatoxin Occurrence Data
Twenty-nine pistachio samples were analyzed for AFB 1  Contamination data are presented analytically in Table 2; the contamination data from each orchard sample were compared to the prediction of the model for that specific weather station and year. The percentage of field samples exceeding the legal limit of AFB 1 contamination set at 12 µg/kg was 40% for 2014, 31% for 2015, and 2% for 2017, and no samples exceeded the legal limit in 2016, 2018, and 2019.

Predictive Model
The relational diagram ( Figure 1) that describes the infection cycle of A. flavus on pistachio starts with the inoculum source, which is not quantified in the model. The produced spores (SoI) are dispersed (DispR) and they arrive on pistachio bunches (DSoP). The spores' flux is described by a dispersal rate (DispR). When the environmental conditions are suitable, the spores germinate, and the fungus can grow (GoP) and infect (IP) early-developed nuts according to the germination and growth rate (respectively, GeR and GrR). If pistachio is in a susceptible growth stage (GS), from flowering to splitting the fungal growth continues in pistachio nuts, as well as AFB 1 production (AFB 1 -I), regulated by the AF production rate (AFB 1 R).

Probability of Aflatoxin Contamination
According to the logistic regression, parametrized with dataset 1 using the AFB1-I calculated by the predictive model as the independent variable (Table 3), the probability p = 0.5 corresponded to AFB1-I = 1153. Comparing the predicted not contaminated samples with the observed data for the internal validation, intended as each field sample collected in the 3-year period 2014-2016 (43 samples in total, Table 4), 55.8% of the prediction was correct (predicted 0 and observed 0), while 16.3% was a false negative prediction (predicted 0 and observed 1). Regarding the prediction of contaminated samples, 16.3% were correctly predicted (predicted 1 and observed 1) and 11.6% were overestimated (predicted 1 and observed 0). The global accuracy of the model was 72.1% (55.8% + 16.3%). The accuracy refers to the percentage of individual samples per year as the model denoted the potential risk for contamination. The robustness of the model is indicated through the percentage of correctly predicted samples; the higher the number of correctly predicted samples, the higher the robustness of the model. Table 3. Parameters (b and c) and statistics of the logistic regression applied to predict the probability of having pistachio samples contaminated above 12 μg of aflatoxin B1 per kg of nuts as a function of the output of the predictive model (AFB1-I = aflatoxin index). The probability p = 0.5 corresponded to AFB1-I = 1153; AFB1-I < 1153 corresponded to p < 0.5 and the related pistachio samples were predicted with a contamination lower than 12 μg/kg (0 = not contaminated). On the contrary, p ≥ 1153 samples were considered contaminated (1= contamination above 12 μg/kg).

Parameters
Parameters  The AFB 1 cumulative index AFB 1 -I, intended as the AFLA-pistachio output, was computed using the meteorological data as input from flowering to crop harvest ( Table 2)

Probability of Aflatoxin Contamination
According to the logistic regression, parametrized with dataset 1 using the AFB 1 -I calculated by the predictive model as the independent variable (Table 3), the probability p = 0.5 corresponded to AFB 1 -I = 1153. Comparing the predicted not contaminated samples with the observed data for the internal validation, intended as each field sample collected in the 3-year period 2014-2016 (43 samples in total, Table 4), 55.8% of the prediction was correct (predicted 0 and observed 0), while 16.3% was a false negative prediction (predicted 0 and observed 1). Regarding the prediction of contaminated samples, 16.3% were correctly predicted (predicted 1 and observed 1) and 11.6% were overestimated (predicted 1 and observed 0). The global accuracy of the model was 72.1% (55.8% + 16.3%). The accuracy refers to the percentage of individual samples per year as the model denoted the potential risk for contamination. The robustness of the model is indicated through the percentage of correctly predicted samples; the higher the number of correctly predicted samples, the higher the robustness of the model. Table 3. Parameters (b and c) and statistics of the logistic regression applied to predict the probability of having pistachio samples contaminated above 12 µg of aflatoxin B1 per kg of nuts as a function of the output of the predictive model (AFB 1 -I = aflatoxin index). The probability p = 0.5 corresponded to AFB 1 -I = 1153; AFB 1 -I < 1153 corresponded to p < 0.5 and the related pistachio samples were predicted with a contamination lower than 12 µg/kg (0 = not contaminated). On the contrary, p ≥ 1153 samples were considered contaminated (1= contamination above 12 µg/kg).

Parameters
Parameters

External Model Validation
In dataset 2 (2017-2019, Table 4), 98.9% of the samples were classified as not contaminated (AFB 1 content in pistachio nut lower than 12 µg/Kg), and the remaining 1.1% were classified as contaminated. According with the probabilistic equation developed, 96.7% of the samples were predicted as not contaminated and 3.3% as contaminated. The contingency table shows that 95.6% of the samples were classified as not contaminated (predicted 0 and observed 0) and none as contaminated (predicted 1 and observed 1). Almost 3% of pistachio nut samples were overestimated (predicted 1 and observed 0), and only 1.1% of samples were underestimated (predicted 0 and observed 1). The global accuracy of predictions with dataset 2 not used for the development of the probabilistic function was 95.6%.

Discussion
Predictive models are nowadays very useful tools supporting modern IPM for managing plant pathogens and the mycotoxins they produce [36]. Pistachio nuts are one of the commodities contaminated by aflatoxins [12], and the meteorological and contamination data suggest that the Mediterranean basin is an area highly prone to aflatoxin contamination. Global warming is also significantly altering the distribution of temperature variability and extremes and precipitation patterns, with a significantly higher frequency of exceptionally unfavorable years for several crops. Therefore, it is very important to organize a predictive system for AF risk in pistachios especially under climate change scenarios, as recently stressed by Battilani and Camardo Leggieri [34] for aflatoxins in maize. Therefore, research is focusing on preharvest, this being the growing period of crops crucial to determine "if" mycotoxin contamination will take place and "how high" it will be at harvest [37]. As a result, supporting producers via predictive models is the ultimate goal of the new generation decision support agricultural systems.
Despite the fact that pistachio cultivation is considered to be a minor crop in Greece as well as in several other regions of the world, the importance of the crop is underlined by the high added product value for producers. However, pistachio farmers often face an aflatoxin contamination problem. As a result, this problem not only leads to uncertainty in the consumer perception of consuming pistachios, but also creates additional production costs and income loss for producers, distributors, and other stakeholders.
In the frame of the current study, the phenology of pistachio nut cultivation was investigated. The crop phenology of pistachio cultivation appears to have been insufficiently studied. The lack of recorded data imposes the need for research regarding the tree's phenology.
The growth degree day approach, as described in the literature [38,39], was selected to describe the tree's phenology in the current study using on site observations combined with meteorological data. A model for crop phenology was developed; it worked reasonably well for several consecutive years for Aegina island, with a day interval of 1.4-4.0 days in the first 2 years and one a bit wider in the third year (4.0-6.8 days). Therefore, the pistachio phenology can be properly predicted as support for the predictive model.
The pistachio variety Aeginis, widely cultivated in Greece, appears to be better adapted to higher winter temperatures and lower chilling accumulation when compared to other pistachio varieties such as Kalle-Ghuchi and the Tunisian cultivar Mateur [39][40][41]. This is also in accordance with data regarding the chilling requirements of pistachio trees in California's Central Valley, where trees are not adapted to such high winter temperatures as in Greece [42].
In the current work, the developed model for the phenology of the crop was combined with an already existing model for the prediction of the risk of A. flavus infection in maize and aflatoxin contamination above the legal limit [33]. Specifically, the model accuracy in the model development and internal validation was 72.1%, with 16.3% underestimates and 11.6% overestimates. This is comparable with the predictive model performances in other pathosystems, such as Fusarium spp in wheat [31] and A. flavus in maize [33,43].
The model performance exceeds the accuracy of 95.6% correctly classified samples in the external validation, with 3.3% of pistachio nut samples overestimated and only 1.1% of samples underestimated. This result is very good, considering that AFLA-pistachio is a mechanistic model using as input only meteorological data. However, it is important to highlight that a higher number of samples contaminated above the legal limit should be used in order to confirm the robustness of the developed model.
The AFB 1 -I values given by the model varied from year to year, which leads to the conclusion that the index is sensitive to changes in weather conditions. The model demonstrated the highest AFB 1 -I values, specifically 1014 and 1043, in years 2014 and 2015, respectively, whereas in 2016 the AFB 1 -I was 879. The warmest year of the three of dataset 1 was 2016. On the other hand, the data for the period 2017-2019, which have been used as dataset 2 for external validation (comparison between the model predictions and real observations of AF contamination), led to lower AFB 1 -I values, varying from 349 to 900. The coldest year of all the 6 years studied was 2019, and it also had the lowest model AFB 1 -I value (AFB 1 -I = 349). Based on these numbers, it is clear that the perception of warmer year alone is not sufficient to support a high aflatoxin contamination risk; the predictive model confirms its role in supporting the risk prediction.
The threshold for the classification of the samples is set at 12 µg/kg of AFB 1 , which is the legal limit for pistachio nuts destined for human consumption. Based on a three-year external validation dataset, the AFLA-pistachio model performance (95.6%) is considered to be high compared to other existing developed models for mycotoxin contamination prediction, like AFLA-maize with 73% accuracy [33] or 84% accuracy for the prediction of deoxynivalenol in wheat [44], taking into consideration that the validation was conducted based on a 3-year contamination dataset. It has also been well demonstrated that several factors, such as the crop variety and its characteristics, can affect aflatoxin contamination [45]. In the current research, such factors were not taken into consideration and were not included in the development of the model because the variety was common to all the orchards assessed. Although the existing inoculum of A. flavus in each pistachio orchard may vary from field to field, this factor was neither taken into consideration nor included in the model development. This is because the presence of A. flavus inoculum in the field is a guarantee that contamination is possible and also because of the difficulty and uncertainty in its assessment in each orchard. Additionally, it is important to underline that agricultural practices can significantly impact the fungal growth and mycotoxin accumulation in the complex interactions in the host plant-fungi-environment, as stated for maize [46] and wheat [47]. High interaction between all the agricultural practices applied can make it very difficult to understand the effects of each operation and even harder to model it with mathematical functions in order to be included in a mechanistic model [37]. In fact, also in the most studied pathosystem related to mycotoxin, Fusarium head blight in small grains, "agronomic programmes" or "arrays of growing techniques" are mentioned, more than single actions, to explain the risk of contamination [48,49]. Empiric approaches were instead used to account for the cropping system in maize [50].
Depending on the agronomical characteristics of each variety cultivated, the early split of pistachio nuts can be controlled up to a certain level by the application of good agricultural practices such as appropriate irrigation. After the naturally occurring split shell, during the maturity stage of the nuts the crop enters a susceptible period where the nuts are more prone to be aflatoxin-contaminated. As the model output provides a daily prediction, it is clearly shown that the date of harvest is a crucial factor affecting the aflatoxin contamination of pistachio nuts. The model predicts a high-risk condition, and therefore it indicates the appropriate time of harvest, which is a very useful suggestion for farmers. This can reduce the aflatoxin contamination of pistachio nuts at preharvest and also indicate the lots with the highest risk of contamination due to the field prediction which require more careful postharvest management.
As such a type of support has been becoming more available to farmers during recent years, the application of agronomical practices and plant protection products can be guided by predictive models. This guidance can be achieved with Decision Support Systems (DSS) that exploit the output of predictive models and provide farmers with valuable information [51]. The model developed in this study could be the core of a DSS, aimed at improving pistachio nut cultivation and reducing the inflow of chemicals in the food chain production. This is crucial, nowadays, in the widest context of IPM, which is fundamental for ensuring agricultural productivity while maintaining economic and environmental sustainability [52,53].
As far as we know, this is the first attempt to transfer such a mechanistic model from one crop to another sharing problems caused by a common fungus. The results are very encouraging, with a model accuracy comparable to the original model developed for A. flavus in maize. The external validation was excellent, but both the limited number of official samples obtained for the validation in the 3-year period 2017-2019 (87 in total) and their distribution respect to the fixed threshold (all except 1 below the contamination threshold) make this result not conclusive. Nevertheless, the model transfer apparently works; this preliminary validation must be supported by further data to state the robustness of the model before its delivery as a support tool for farmers. Further, model transfer may lead, in the future, to the association of the A. flavus sub-model with other crops of interest for aflatoxin contamination. This should be considered a valuable asset in the prevention of the mycotoxin contamination of crops and food products, with the practical application of preharvest models also in areas with a high potential risk of mycotoxin contamination but poor resources for research and modelling.

Collection of Meteorological Data
The daily meteorological data of air temperature (T, • C) and rain (R, mm) were collected from the website www.meteo.gr, supported by the National Observatory of Athens for Aegina island, a region with considerable PDO pistachio production in Greece. The daily relative humidity (RH, %) data were kindly supplied by the Greek National Meteorological Service. The island is located close to Athens in the prefecture of Piraeus, region of Attica, in the Saronic Gulf. The geographical coordinates of Aegina city, the capital of the island, are 37 • 44 51.42 ' and 23 • 25 44.97 '. The reference period was 2014-2019.

Description and Computation of Growth Stages
The chilling hours were calculated based on the Crossa-Raynaud Model [1], the most commonly used model for the computation of pistachio growth stages; the data from 1st October to 31st March in the seasons 2013-2014, 2014-2015, and 2015-2016 were used.
A very crucial point when using this model for the description/prediction of crop phenology is the definition of the temperature threshold; for pistachio trees, the threshold is defined at 7 • C and a temperature below 7 • C is required during winter for the break of dormancy for most varieties [40,54].
The chilling period is then followed by heat requirements from April until harvest [4]. The Growth Degree Days (GDD) for each growth stage of pistachio were calculated in this study according to [7]. Briefly, the useful GDD (T > 7 • C; [40,54]) were added daily throughout the growing season starting from April.

Field Sampling and Aflatoxin Occurrence Data
The field sampling was performed according to European protocols for sampling (Commission Regulation (EC) No 401/2006) by the Agricultural Association of Pistachio Growers of Aegina from 2014 to 2019 between September and October. All the pistachio nuts were initially dehulled and properly dried. The determination of aflatoxin B 1 (AFB 1 ) contamination was carried out by authorized analytical laboratories that trade in food safety issues such as AgroLab RDS, SkyLab Med, and Tsakalidis Analysis and Testing. High Phase Liquid Chromatography (HPLC) with a fluorescence detector was used, and the limit of detection (LOD) was <0.1 µg/kg according to the certified external laboratories that conducted the analysis. The field data were organized in two different datasets to be used for (i) the development of the probability of AFB 1 contamination and (ii) for the model validation, as follows:

Predictive Model
The relational diagram was developed and adapted based on an already existing diagram [33] and was named AFLA-pistachio ( Figure 1). Briefly, the relational diagram was organized according to the principles of "systems analysis" [55]. The status of the fungus is represented by boxes, intended as state variables. The flow from one state to the following is driven by constants/parameters or intermediate variables, basically driven from weather data or crop data. The rate variables are represented by "valves", described by mathematical functions. Variables operating in the predictive model are finally linked in a coherent mathematical framework to calculate the final output of the model-i.e., the AFB 1 cumulative index (AFB 1 -I). AFLA-pistachio is a weather based predictive model; meteorological data (T, R and RH) were the input data used both to predict the crop phenology and A. flavus behavior.

Probability of AFB 1 Contamination
The AFB 1 -I computed for the meteorological station and year was associated with the AFB 1 contamination detected in the related pistachio samples, according to Battilani et al. (2013) [33] for the dataset 1 (2014-2016). Contamination data on AFB 1 in pistachio at harvest were shared in 2 exclusive groups based on the threshold 12 µg/kg, the legal limit set for pistachio nuts to be subjected to sorting, or other physical treatment, before human consumption or use as an ingredient in foodstuffs (European Commission Regulation (EU) no. 165/2010 amending Regulation (EC) no. 1881/2006). A binary logistic regression (p = 1/(1 + exp −(c+b*AFB1-I) ) was developed using as the dependent variable the AFB 1 contamination of the samples and as the independent variable the AFB 1 -I generated as the output of AFLA-pistachio [33]. This approach estimates the probability of an event to occur (AFB 1 contamination below/above the legal limit, binary data) based on an independent variable (the index given by the model at the end of the growing season). The probability (P) can range from 0 to 1; the event is considered as occurring when p > 0.5 and not occurring when p ≤ 0.5. The logistic regression module of IBM SPSS Statistics (version 25.0) was used to estimate the parameters (b and c, Table 3) of the logistic equation.

Validation
To confirm the ability of AFLA-pistachio to accurately predict the occurrence of AFB 1 in pistachios, the estimated calculated probabilities were compared to the observed aflatoxin contamination field data as a measure of the goodness of fit of the predictions. The performance of the prediction was evaluated using both sets of field data; dataset 1 as the input data to estimate the parameter of logistic regression (internal validation) and the independent dataset 2 (years 2017-2019, external validation) [56].