The Effect of Cold Periods on the Biological Cycle of Marchalina hellenica

Simple Summary The most important honeydew-producing insect in Greece is Marchalina hellenica (Coccoidea: Marchalinidae), which is a parasite on pine trees. The current work is part of an ongoing research project aiming to provide knowledge on honeydew-producing insects and the impact of critical factors (climate, beekeeping manipulations, honeydew-producing insect phenology) on honeydew honey production. Empirical evidence indicates that among the weather factors, the most important one, at least for spring honeydew secretions, appears to be temperature and, more specifically, the existence of cold winter days. Presently, we investigate the effect of cold periods in February on the life cycle of Marchalina hellenica. Our primary goal is to help beekeepers plan the timely exploitation of honeydew secretions of pine trees. Such a potential will be beneficial for beekeepers, the rural economy, and forest protection. It should be noted that current results highlight the impact of climate change in the field of entomology, and they indicate that the life cycle of Marchalina hellenica is expected to be drastically shorter. Abstract Climate change is considered a major factor affecting honeybees’ behavior and productivity with major consequences in both honey and agricultural production. Many research studies have expressed serious concerns about the mass losses of bee colonies and the role of bees as pollinators, while others have underlined important issues for the impact of the increase in temperature on honeybee abundance and honey yields. In the present work, we draw our attention to Marchalina hellenica, which is the most important honeydew-producing insect in Greece. A statistically significant forecasting model for the effect of cold periods in February on the life cycle of the insect is constructed, with the aid of the Cumulative Logit Model and the theory of runs. The forecasting model may help beekeepers plan the timely exploitation of honeydew secretions of pine trees, which will be beneficial for beekeepers, the rural economy, and forest protection. The new suggested model also indicates that, in view of the climate change scenarios seen in the literature, the life cycle of M. hellenica is expected to be drastically shorter.


Introduction
Climate change is considered a major factor affecting honeybees' behavior and productivity with major consequences in both honey and agricultural production [1]. Even though the main role of beekeeping is honey production [2,3], the maintenance of honeybees contributes significantly to the pollination activity. It should be noted that honeybees, as pollinators, provide valuable services to agricultural [4,5] and natural ecosystems with significant impacts on biodiversity and food security [6]. Therefore, many research studies have expressed serious concerns about the mass losses of bee colonies and the role of bees

Materials and Methods
M. hellenica is a normal non-cyst-forming coccid, with the life cycle of the females having three immature stages or instars and that of the males having four immature stages [19]. A detailed description of the life cycle of M. hellenica can be found in [19].
Regarding the study area, Rhodes is the largest Greek island in the Dodecanese and is located in the south-eastern part of the Aegean Sea. The main honey types produced are pine and thyme honey. Pine forests cover most of the island, reaching, in some places, up to the sea. Kalithies is at the north-eastern part of the island, while Embonas is at the center-west (Figure 1).

Materials and Methods
M. hellenica is a normal non-cyst-forming coccid, with the life cycle of the females having three immature stages or instars and that of the males having four immature stages [19]. A detailed description of the life cycle of M. hellenica can be found in [19].
Regarding the study area, Rhodes is the largest Greek island in the Dodecanese and is located in the south-eastern part of the Aegean Sea. The main honey types produced are pine and thyme honey. Pine forests cover most of the island, reaching, in some places, up to the sea. Kalithies is at the north-eastern part of the island, while Embonas is at the center-west (Figure 1). We used the data collected by Gounari [20] regarding the biological cycle of the insect in the two areas of interest (for a detailed presentation of sampling procedure and data collection, see [20]). We should reiterate that temperature constitutes a parameter of major importance when observing the life cycle of M. hellenica. Therefore, the winter temperature in the areas of interest, Kalithies and Embonas, during the period from 2014 to 2019 was monitored. The climatic variable of temperature for the two regions was recorded in two different weather stations of IERSD/NOA (Institute for Environmental Research, National Observatory of Athens). More specifically, interest was focused on the minimum daily winter (December to February) temperature. The basic descriptive statistics of temperature are presented in Table 1. Next, the process of completion of adult insects was observed using the available data from 2014 to 2019 for the island of Rhodes. The fortnight of integration as well as other characteristics are presented in Table 2. More specifically, the fortnight of adult M. hellenica was categorized in = 5 categories. The first category (1) corresponds to the occasion when M. hellenica becomes an adult during the second fortnight of March, category (2) when the completion takes place during the first fortnight of April, category (3) when We used the data collected by Gounari [20] regarding the biological cycle of the insect in the two areas of interest (for a detailed presentation of sampling procedure and data collection, see [20]). We should reiterate that temperature constitutes a parameter of major importance when observing the life cycle of M. hellenica. Therefore, the winter temperature in the areas of interest, Kalithies and Embonas, during the period from 2014 to 2019 was monitored. The climatic variable of temperature for the two regions was recorded in two different weather stations of IERSD/NOA (Institute for Environmental Research, National Observatory of Athens). More specifically, interest was focused on the minimum daily winter (December to February) temperature. The basic descriptive statistics of temperature are presented in Table 1. Next, the process of completion of adult insects was observed using the available data from 2014 to 2019 for the island of Rhodes. The fortnight of integration as well as other characteristics are presented in Table 2. More specifically, the fortnight of adult M. hellenica was categorized in J = 5 categories. The first category (1) corresponds to the occasion when M. hellenica becomes an adult during the second fortnight of March, category (2) when the completion takes place during the first fortnight of April, category (3) when the completion takes place during the second fortnight of April, category (4) when the completion takes place during the first fortnight of May, and category (5) when the completion takes place during the second fortnight of May. Due to some differences between the two regions regarding the time when M. hellenica completes its biological cycle, an independent samples t-test to check whether this difference was statistically significant was performed. The results showed that there was no statistical difference (p > 0.050) in the fortnight of completion between the two regions, indicating that the data should not be differentiated per region for the purposes of this analysis.

Statistical Analysis
During recent decades, a wide range of problems in several research areas such as Agriculture, Biology, Finance, Meteorology, Quality Control, and Systems' Reliability (see e.g., [22,23]) has been modeled by classifying the experimental trials in two exclusive categories, considering sequences of binary trials (with values 0 or 1) and studying the sequences of outcomes. This study usually involves searching for the greatest concentration of outcomes of a specific type. Such a concentration is traditionally measured by the enumeration of runs of (k ≥ 1) 1s. In the present work, an ordinal regression model employing runs was proposed. Depending on the special features of the specific application under study, several ways for counting runs have been proposed in the literature. Currently, we should draw attention to the Type I counting scheme [22,24]. According to the Type I counting scheme, a run of length k is registered when k consecutive 1 s are observed, and the (k + 1)th consecutive 1 is considered as the first 1 of the next run.
As an example of the Type I counting scheme for runs of length k, let us assume that by examining 18 consecutive days, the sequence 011111101110111110 of days that have been classified as cold or not-cold ones (1 stands for a cold day and 0 stands for a not-cold day according to prespecified criteria) has been observed. Then, we may count, for example, six different runs of length k = 2, four different runs of length k = 3, or two different runs of length k = 4 (see . the completion takes place during the second fortnight of April, category (4) when the completion takes place during the first fortnight of May, and category (5) when the completion takes place during the second fortnight of May. Due to some differences between the two regions regarding the time when M. hellenica completes its biological cycle, an independent samples t-test to check whether this difference was statistically significant was performed. The results showed that there was no statistical difference (p > 0.050) in the fortnight of completion between the two regions, indicating that the data should not be differentiated per region for the purposes of this analysis.

Statistical Analysis
During recent decades, a wide range of problems in several research areas such as Agriculture, Biology, Finance, Meteorology, Quality Control, and Systems' Reliability (see e.g., [22,23]) has been modeled by classifying the experimental trials in two exclusive categories, considering sequences of binary trials (with values 0 or 1) and studying the sequences of outcomes. This study usually involves searching for the greatest concentration of outcomes of a specific type. Such a concentration is traditionally measured by the enumeration of runs of (k ≥ 1) 1s. In the present work, an ordinal regression model employing runs was proposed. Depending on the special features of the specific application under study, several ways for counting runs have been proposed in the literature. Currently, we should draw attention to the Type I counting scheme [22,24]. According to the Type I counting scheme, a run of length k is registered when k consecutive 1 are observed, and the (k + 1)th consecutive 1 is considered as the first 1 of the next run.
As an example of the Type I counting scheme for runs of length , let us assume that by examining 18 consecutive days, the sequence 011111101110111110 of days that have been classified as cold or not-cold ones (1 stands for a cold day and 0 stands for a notcold day according to prespecified criteria) has been observed. Then, we may count, for example, six different runs of length = 2, four different runs of length = 3, or two different runs of length = 4 (see .    the completion takes place during the second fortnight of April, category (4) when the completion takes place during the first fortnight of May, and category (5) when the completion takes place during the second fortnight of May. Total 12 Due to some differences between the two regions regarding the time when M. hellenica completes its biological cycle, an independent samples t-test to check whether this difference was statistically significant was performed. The results showed that there was no statistical difference (p > 0.050) in the fortnight of completion between the two regions, indicating that the data should not be differentiated per region for the purposes of this analysis.

Statistical Analysis
During recent decades, a wide range of problems in several research areas such as Agriculture, Biology, Finance, Meteorology, Quality Control, and Systems' Reliability (see e.g., [22,23]) has been modeled by classifying the experimental trials in two exclusive categories, considering sequences of binary trials (with values 0 or 1) and studying the sequences of outcomes. This study usually involves searching for the greatest concentration of outcomes of a specific type. Such a concentration is traditionally measured by the enumeration of runs of (k ≥ 1) 1s. In the present work, an ordinal regression model employing runs was proposed. Depending on the special features of the specific application under study, several ways for counting runs have been proposed in the literature. Currently, we should draw attention to the Type I counting scheme [22,24]. According to the Type I counting scheme, a run of length k is registered when k consecutive 1 are observed, and the (k + 1)th consecutive 1 is considered as the first 1 of the next run.
As an example of the Type I counting scheme for runs of length , let us assume that by examining 18 consecutive days, the sequence 011111101110111110 of days that have been classified as cold or not-cold ones (1 stands for a cold day and 0 stands for a notcold day according to prespecified criteria) has been observed. Then, we may count, for example, six different runs of length = 2, four different runs of length = 3, or two different runs of length = 4 (see Figures 2-4).     the completion takes place during the second fortnight of April, category (4) when the completion takes place during the first fortnight of May, and category (5) when the completion takes place during the second fortnight of May. Total 12 Due to some differences between the two regions regarding the time when M. hellenica completes its biological cycle, an independent samples t-test to check whether this difference was statistically significant was performed. The results showed that there was no statistical difference (p > 0.050) in the fortnight of completion between the two regions, indicating that the data should not be differentiated per region for the purposes of this analysis.

Statistical Analysis
During recent decades, a wide range of problems in several research areas such as Agriculture, Biology, Finance, Meteorology, Quality Control, and Systems' Reliability (see e.g., [22,23]) has been modeled by classifying the experimental trials in two exclusive categories, considering sequences of binary trials (with values 0 or 1) and studying the sequences of outcomes. This study usually involves searching for the greatest concentration of outcomes of a specific type. Such a concentration is traditionally measured by the enumeration of runs of (k ≥ 1) 1s. In the present work, an ordinal regression model employing runs was proposed. Depending on the special features of the specific application under study, several ways for counting runs have been proposed in the literature. Currently, we should draw attention to the Type I counting scheme [22,24]. According to the Type I counting scheme, a run of length k is registered when k consecutive 1 are observed, and the (k + 1)th consecutive 1 is considered as the first 1 of the next run.
As an example of the Type I counting scheme for runs of length , let us assume that by examining 18 consecutive days, the sequence 011111101110111110 of days that have been classified as cold or not-cold ones (1 stands for a cold day and 0 stands for a notcold day according to prespecified criteria) has been observed. Then, we may count, for example, six different runs of length = 2, four different runs of length = 3, or two different runs of length = 4 (see Figures 2-4).    In the present work, the theory of runs was employed along with an ordinal regression model, and the number of cold periods in February is a significant parameter to the life cycle of M. hellenica.
Ordinal regression models are used when the response variable is ordinal with more than two categories. The present work used the Proportional Odds Model (POM), which assumes that the relationship between predictor and response variables is independent of the response variable's category [25].
Let Y be a response variable with J categorical outcomes, denoted by 1, 2, . . . , J, and let x be a vector of covariates (independent/predictor variables). Let π 1 (x), π 2 (x), . . . , π J (x) also be the probabilities of the J-ordered categories of the response variable when covariates have the value x. Then, the Proportional Odds Model is defined as follows: where a j is the intercept in category j and β is a vector of regression coefficients corresponding to x (we denote by β the transpose of the matrix β).
In logit form, the Proportional Odds Model transforms as follows: The a j are increasing in j, as Pr( Y ≤ j|x) increases in j for fixed x, and the logit is an increasing function of this probability [26].
It must be noted that the POM utilizes the cumulative probabilities: and the regression coefficients vector β does not depend on j, meaning that the relationship between Y and x is independent of j.
In addition, the cumulative odds ratio is independent of j and depends only on the difference between the covariate values, In the present work, the vector of independent variables x consists of a single variable, X k , which is obtained through the theory of runs of length k, for different values of the parameter k. As a result, the regression coefficients vector β consists of a single coefficient, β 1 , as well. Based on these, the model becomes, and in logit form, Determining the Variables Gounari et al. [16] employed an ordinal regression model to predict the fortnight of completion of the biological cycle of the honeydew-producing insect, M. hellenica, based on the number of cold days in February. The threshold for classifying a day of February as a cold one was set to the 20% percentage (P 20 ) of the dataset of the minimum daily winter temperature, which was calculated to be equal to 7.3 • C (see Table 1), i.e., any day with a minimum temperature less than 7.3 • C was classified as a cold day. The current study also considered another statistically significant criterion, the number of cold periods in February. In order for the definition of a cold period to be possible, a different classification must be made. Any day in February is classified as a relatively cold day (rcd) if the daily minimum temperature is lower than the 30% percentage (P 30 ) of the dataset of the minimum daily winter temperatures, which, in this case, is 8.5 • C (see Table 1). Therefore, any day in February with a minimum temperature less than 8.5 • C is classified as a rcd. It should be noted that there is strong empirical evidence that the number of days in February when the minimum daily temperature is lower than a specific threshold has an analogous effect on the prolongation of the completion of the biological cycle of M. hellenica, and both temperature thresholds of 7.3 • C and 8.5 • C coincide with the practitioners' experience. A run of k consecutive rcds is referred to as a cold period. It is proven that the number of cold periods in February is a significant parameter to the life cycle of M. hellenica, a predictability that is of vital importance to beekeepers. It should be noted that the criterion of the number of cold days, employed by Gounari et al. [16], is now replaced by the number of runs of consecutive rcds (cold periods). The daily temperature threshold is now relaxed as a run has a cumulative effect. The new results of this work regarding the effect of the number of k (k = 2, 3) consecutive runs of rcds in February on the life cycle of M. hellenica is now presented. Table 3 describes in detail the dependent variable (Y) and the respective independent variable (X k ) that was used in the one-dimensional ordinal regression model for the values k = 2, 3. Moreover, Table 4 gives the realizations of X 2 and X 3 in the experimental data every year from 2014 to 2019.

Results
In this section, the results regarding the three different statistically significant criteria to the life cycle of M. hellenica (number of cold days in February, number of runs of two consecutive rcds in February, number of runs of three consecutive rcds in February) are presented. Before the presentation of the new results regarding the criteria of runs of rcds (cold periods), the results of Gounari et al. [16] regarding the effect of the number of cold days in February (Criterion 1) are presented.

Criterion 1: Number of Cold Days in February
As it is already mentioned in the Introduction, Gounari et al. [16] proved the effect of temperature, and more specifically, the number of cold days in February, to the life cycle of M. hellenica. They employed the POM to prove that for each additional cold day in February, the Odds Ratio (OR) of M. hellenica to complete its cycle in a specific fortnight decreased by 52.15%. We now proceed to the new results regarding the criteria of cold periods in February.

Criterion 2: Number of Runs of Two Consecutive Rcds in February
The POM was employed with the dependent variable Y and the independent variable X 2 . The assumption that all the logit surfaces, lines in this case, are parallel (i.e., proportional odds assumption) was tested and the null hypothesis of parallelism of the logit lines with respect to x was accepted (p > 0.050). The POM was tested, and it was statistically significant (p < 0.050) ( Table 5). The odds ratio estimate, OR = 0.363, indicated that the odds of M. hellenica to complete its cycle in a specific fortnight were multiplied by 0.363 for each additional run of two consecutive rcds in February. In other words, for each additional run of two consecutive rcds in February, the odds of M. hellenica to complete its cycle in a specific fortnight decreased by 63.7%. Therefore, the number of two consecutive rcds was a more critical factor to the completion time of the insect's biological cycle than the number of cold days. A diagram of the observed cumulative frequencies for the fortnight of completion of the biological cycle of M. hellenica ( Figure 5) is also given. Each curve of the plot corresponds to the number of runs of rcds, from zero to six, which was recorded for the month of February. The cumulative probability curves on the left side of the plots correspond to a lower number of runs of rcds. More specifically, these curves express the low probability that M. hellenica will become an adult in one of the subsequent reference fortnights. Figure 5 reveals that as the number of runs of two consecutive rcds increased, the probability of completion of the biological cycle of M. hellenica in a following fortnight also increased. Considering Equation (2) and the estimated coefficients in Table 5, the cumulative probabilities of adult M. hellenica biological cycle completion and, consequently, the probability of every specific fortnight were computed (Table 6).     Note: The column of CDF indicates the probability of the biological cycle completion of M. hellenica lower or equal to a specific fortnight j (γ j ). The column PMF indicates the probability of the biological cycle completion of M. hellenica equal to a specific fortnight j and can be calculated as: γ j − γ j−1 .

Criterion 3: Number of Runs of Three Consecutive Rcds in February
The POM was employed with the dependent variable Y and the independent variable X 3 . The assumption of parallelism was successfully tested (p > 0.050). The POM was proved to be statistically significant (p < 0.050) ( Table 7). According to Table 7, the addition of a run of three (consecutive) rcds in February decreased the OR of M. hellenica to complete its biological cycle within a specific fortnight by 79%. In comparison with the previous case, the increase in the run's length per rcd led to a further decrease in the corresponding OR by 15.3%.  Figure 6 and Table 8 are now given. In Figure 6, the cumulative probability curves corresponding to six runs of three rcds indicated that the probability of completion of the biological cycle of M. hellenica in one of the initial reference fortnights was almost zero.    Note: The column of CDF indicates the probability of the biological cycle completion of M. hellenica lower or equal to a specific fortnight j (γ j ). The column PMF indicates the probability of the biological cycle completion of M. hellenica equal to a specific fortnight j and can be calculated as: Let µ i denote the average minimum temperature in day i (i = 1, . . . , 28) of February for the area of Embonas of Rhodes Island and all years from 2014 to 2019, and W k denote the number of runs of k consecutive values of µ i less than 8.5 • C. We may calculate that W 2 = 7 and W 3 = 4. In the Discussion, these values are compared to the respective values under climate change scenarios and the impact of climate change is highlighted.

Discussion
Current work is part of an ongoing research project aiming to provide additional knowledge on honeydew-producing insects and the impact of critical factors (climate, beekeeping manipulations, honeydew-producing insect phenology) on honeydew honey production. Its primary goal is to help beekeepers plan the timely exploitation of honeydew secretions of pine trees. Such a potential results in several personal and social benefits. More specifically, the timely exploitation of honeydew secretions of pine trees:

•
Reduces production costs of honey with the simultaneous increase in competitiveness, for a large part of beekeepers, who exploit honeydew insect secretions to collect honey, through streamlining the transportation of beehive colonies. • Reduces transportation costs by reducing beekeepers' transfers for beehives' transport and on-site inspections. It also reduces the stress of vehicles and the likelihood of an accident on the road. • Reduces labor costs due to transportation reduction and on-site inspections.

•
Reduces bee colony losses. Part of the population of bees are lost and/or dying during transport. Frequent inspections and inappropriate harvesting can cause looting, which may result in the loss of entire flocks or in the killing of Queens.

•
Increases production per hive. The choice of region and the timely transfer and removal of bees from one honeydew increases the overall performance of honeydew production.
Presently, the recent work of Gounari et al. [16] is expanded by the addition of extra criteria, which include the number of cold periods in February (number of runs of rcds of length k = 2 and k = 3). The new results of this work were obtained by the combination of ordinal regression and the theory of runs. It should also be noted that, as far as the authors know, this is the first time in academic research that the theory of runs has been combined with ordinal regression models. The authors strongly believe that such a combination can be profitably applied to other scientific fields as well, such as biomedical engineering (see [27]).
The impact of temperature on the biological cycle of M. hellenica can be readily seen in Tables 7 and 8. For instance, in the absence of runs of two consecutive rcds in February or three consecutive rcds in February, M. hellenica will have completed its biological cycle by April's second fortnight with a probability 97.17% and 96.86%, respectively (employing Criterion 1, the respected probability was calculated to be 98.35%, see [16]). When February has four runs of two consecutive rcds, the probability that M. hellenica will have completed its biological cycle by April's second fortnight decreases by about 60% compared to the case when there are no consecutive rcds. On the other hand, if February has four runs of three consecutive rcds, the respective probability now decreases by about 91%. It is also worth noting that when February has two runs of two consecutive rcds, the most probable fortnight of M. hellenica's biological cycle's completion is the first fortnight of April with a probability of around 52%. When February has two runs of three consecutive rcds, the most probable fortnight of M. hellenica's biological cycle's completion is the second fortnight of April with a probability of about 31%.
Therefore, it may be seen that the combination of runs and ordinal regression highlighted and further quantified the effect of temperature to the life cycle of M. hellenica. Moreover, the new theoretical tools enable an estimation of the effect of climate change to the field.
The atmospheric temperature is projected to continue to rise throughout the 21st century in most areas of the planet (see [28] and the references therein). In particular, on the basis of average results over a number of climate model simulations, the average atmospheric temperature is expected to increase by 1.8-4 • C by 2100, depending on developments in the concentration of greenhouse gas emissions [28]. The increase in temperature is expected to be greater at higher latitudes and in continental regions as opposed to the oceans. Therefore, let us consider the following three climate change scenarios: Scenario (i) the average minimum atmospheric temperature is expected to increase by 1.8 • C (optimistic scenario); (ii) the average minimum atmospheric temperature is expected to increase by 2.9 • C (moderate scenario); (iii) the average minimum atmospheric temperature is expected to increase by 4 • C (pessimistic scenario). Let µ Therefore, the values of W (s) k for k = 2, 3 may now be calculated, which are all equal to zero. It may be concluded that the value of W 2 has decreased by seven units and the value of W 3 has decreased by 4, even by assuming the optimistic climate change scenario. In view of the previous results in this section, this drastic decrease in the number of the expected runs indicates that, even in the optimistic scenario, the life cycle of M. hellenica will be significantly shortened. We recall that, according to our model, if, for example, the number of runs of two consecutive rcds reduces by one, the odds of M. hellenica to complete its cycle in a specific fortnight are multiplied by 1/0.363, or they increase by 175%.
Therefore, the new suggested model enables the incorporation of the results in the climate change scenarios seen in the literature. It is currently shown that the life cycle of M. hellenica is expected to be drastically shorter, which will affect the forests, the rural communities, and society in general.