Seasonal Phenology and Climate Associated Feeding Activity of Introduced Marchalina hellenica in Southeast Australia

Simple Summary Giant pine scale, Marchalina hellenica Gennadius (Hemiptera, Marchalinidae), is a sap sucking insect native to the Eastern Mediterranean Basin. In 2014, giant pine scale was detected for the first time in Victoria, Australia, feeding on a new host, Pinus radiata. We studied the life cycle and feeding activity of giant pine scale in Victoria over 32 months, with the aim of drawing comparisons between exotic and native populations. Australian life stages of this pest emerged during similar months to Greek seasonal equivalents, although the timing of Australian life stages differed between years. Insect density and feeding activity on infested trees differed among locations and between generations. Strong evidence was found in support of density and feeding intensity being explained by climatic conditions. The value of Pinus radiata as a food source for this insect may fluctuate with climate conditions. Our findings aim to inform future management efforts for this scale insect, including surveillance strategies and optimal seasons for release of biocontrol agents. Our findings also suggest that the impact of this pest in Australia may be exacerbated by climate change. Abstract Invasive insects pose an increasing risk to global agriculture, environmental stability, and public health. Giant pine scale (GPS), Marchalina hellenica Gennadius (Hemiptera: Marchalinidae), is a phloem feeding scale insect endemic to the Eastern Mediterranean Basin, where it primarily feeds on Pinus halepensis and other Pinaceae. In 2014, GPS was detected in the southeast of Melbourne, Victoria, Australia, infesting the novel host Pinus radiata. An eradication program was unsuccessful, and with this insect now established within the state, containment and management efforts are underway to stop its spread; however, there remains a need to understand the insect’s phenology and behaviour in Australia to better inform control efforts. We documented the annual life cycle and seasonal fluctuations in activity of GPS in Australia over a 32 month period at two contrasting field sites. Onset and duration of life stages were comparable to seasons in Mediterranean conspecifics, although the results imply the timing of GPS life stage progression is broadening or accelerating. GPS density was higher in Australia compared to Mediterranean reports, possibly due to the absence of key natural predators, such as the silver fly, Neoleucopis kartliana Tanasijtshuk (Diptera, Chamaemyiidae). Insect density and honeydew production in the Australian GPS population studied varied among locations and between generations. Although insect activity was well explained by climate, conditions recorded inside infested bark fissures often provided the weakest explanation of GPS activity. Our findings suggest that GPS activity is strongly influenced by climate, and this may in part be related to changes in host quality. An improved understanding of how our changing climate is influencing the phenology of phloem feeding insects such as GPS will help with predictions as to where these insects are likely to flourish and assist with management programs for pest species.


Introduction
Invasive insects pose a serious threat to forest health, ecosystem stability, and even public health in areas where they establish [1][2][3][4]. Field observations of an exotic insect's life cycle can be important to understanding their population dynamics in the new environment and when control efforts could be maximised [5,6]. This would be particularly important for understanding species establishing in both the northern and southern hemispheres, due to the inverse of calendar season. Phenological studies might include documenting insect activity patterns and emergence of developmental stages (e.g., [7,8]) and, in the case of herbivorous insects, matching these to fluctuations in host quality [9,10]. For phloem feeding insects, mapping the timing of feeding life stages and feeding activity can be used to predict when impacts are greatest [11][12][13]. In Australia, insect forest pests are being detected at an accelerating rate [14,15], many of which are phloem-feeding insects.
In 2014, the giant pine scale, Marchalina hellenica Gennadius, (Hemiptera, Marchalinidae) was first detected in Victoria, Australia, feeding on the novel host Pinus radiata, or California Monterey pine [16]. Marchalina hellenica is endemic to the eastern Mediterranean basin [17][18][19]; its establishment on P. radiata (which originates from North America) poses a significant threat to Australian softwood plantations. Pinus radiata represents the vast majority of Australia's more than 1 million hectares of softwood plantation [20]. Despite the widespread use of P. radiata as a plantation species and an ornamental, in public gardens and farm windrows, across Australia, the documented distribution of M. hellenica is currently restricted to southeast Melbourne [16].
Little is known about how the life cycle of M. hellenica in its new Australian habitat compares to its endemic range. The insect's three immature nymphal instars acquire nutrients for development from the sap of pine trees [21]. The nymphal instars feed on phloem by inserting syringe-like stylets into twigs, branches, trunks, and exposed roots [22][23][24]. As with many phloem-feeding insects, excess carbohydrates and undesirable compounds such as insecticides [25] that are ingested must be excreted as honeydew [26,27]. Marchalina hellenica does not feed on host trees all year round, and only the nymphal instar stages possess a functioning stylet and produce honeydew [28]. Marchalina hellenica likely originates from Mount Carmel in Northern Israel, the origin of its primary host, P. halepensis [29]. Honeydew produced by M. hellenica is collected by honeybees across the Eastern Mediterranean and accounts for 60% (~15,000 tonnes) of annual honey produced in Greece [30,31]. While this is important for apiculture in Greece, it also highlights the phloem feeding capability of this insect. Its impact as an introduced pest to Australian forestry could be severe, particularly if feeding stages are prolonged and/or populations flourish in the absence of important native predators, such as the predatory fly, Neoleucopis kartliana Tanasijtshuk (Diptera: Chamaemyiidae) [32]. In addition to exuding honeydew, M. hellenica secretes cotton-like wax filaments [21] termed flocculent, which may provide a hydrophobic layer and microclimate to protect scales from desiccation and flooding of feeding sites [33][34][35].
Insect density and availability of honeydew are key indicators of the resources expropriated by M. hellenica. Explaining variation in these indicators is important for understanding the insect's activity patterns. As an ectothermic phloem feeder, M. hellenica activity is likely related to climate via effects on host nutritional quality [36,37]. As an insect sheltering under bark and flocculent, climate conditions directly experienced by M. hellenica may therefore provide inferior explanations for variation in insect activity than the climate experienced by host plants [9,38]. Understanding how M. hellenica responds to the climate may help predict the potential spread of this pest in Australia or the impact as the climate in the Mediterranean [39,40] and Pacific region [41,42] trend towards more unpredictable conditions through climate change.
In this study, we used field-sampling to document the seasonal phenology of M. hellenica in southeast Melbourne, Victoria, across two and a half insect generations. As these represent the only known M. hellenica populations in the southern hemisphere or outside Europe, we explored whether major life cycle milestones (e.g., nymph emergence Insects 2023, 14, 305 3 of 22 or moults) in Melbourne aligned with equivalent seasonal periods to East-Mediterranean conspecifics (i.e., 6 calendar-month difference). We also investigated whether the timing of instar stages and honeydew production (as the proportion of honeydew producing insects, HPI) followed the same seasonal pattern as Mediterranean conspecifics. As a univoltine (one generation per year) and semelparous (mortality following oviposition) insect, density of M. hellenica was predicted to reach the generational minimum during the adult stage and maximum at nymph emergence. Finally, we explored the prediction that variation in M. hellenica activity would be best explained by the environmental conditions experienced by the host tree, rather than the conditions directly experienced by M. hellenica.

Study Sites
Infested M. hellenica branches were collected every second week (fortnightly) from January 2019 to October 2021, or weekly during moulting periods, from two field sites in south-east Melbourne, Australia: Cardinia Reservoir Park (CRP) in Cardinia (managed by Parks Victoria, a government agency of the state of Victoria, Australia) ( . Insects were considered to be in a moulting period when more than one life stage was concurrently detected. The two field sites were selected due to (i) the reliable accessibility, (ii) trees not being removed during the three-year survey period, and (iii) differences in habitat structure, specifically tree density and size. Pinus radiata trees at CRP were planted in 1973 to help stabilise soil around the reservoir's banks [43]. The stands of M. hellenica infested P. radiata were planted in plantation style rows and are of a fairly uniform age. The understory is smothered in dead P. radiata needles, with little to no grass or other ground-covering vegetation and an open mid-story consisting of native Pittosporum undulatum. Forest surrounding these P. radiata stands were dominated by Eucalyptus and Acacia species. By contrast, DR is a small publicly accessible park located in Harkaway. There are only a handful of ornamental P. radiata at DR, all of which were found to be infested. The DR pine trees were much larger than those sampled at CRP. The ground of this reserve is covered by introduced grasses. Other large ornamental trees are present, such as Eucalyptus leucoxylon and Acacia mearnsii; the canopies of P. radiata and other trees at DR do not form a closed canopy.

Environmental Conditions
Local temperature and relative humidity (RH %) at CRP were recorded using data loggers (accurate to 1.0 • C and resolved to 0.125 • C; DS1923, iButton ® , Whitewater, WI, USA). At each site, a data logger was either affixed to the surface of an infested tree trunk (exposed) or wedged into a P. radiata bark fissure and under naturally present flocculent, where sufficient flocculent was present to shelter the data logger (fissure) as close to 1.5 m above the south-facing base of the tree as feasible. Data loggers recorded temperature and relative humidity every 30 min, with the exception of a few recordings missed during monthly data retrieval and logger resets. Vapour pressure deficit (VPD), recorded in kPa, was calculated from temperature and relative humidity via an established adaptation of the Magnus equation for vapour saturation point [44][45][46]. Exposed records were recorded over 707 days (1 November 2019-7 October 2021) and bark fissure records over 469 days (26 June 2020-7 October 2021).
Temperature ( • C) and solar exposure (MJ m −2 ) records were sourced from the Bureau of Meteorology and recorded at Ferny Creek weather station [47]. The Bureau of Meteorology sourced records are referred to here as 'ambient' climate variables. Daily rainfall (mm) was sourced from Melbourne Water [48] and recorded in Emerald, which was closer to field sites than Ferny Creek weather station and at a similar elevation. Rainfall, solar exposure, and ambient temperature records span the entirety of the honeydew and density data; 934 days (19 March 2019-7 October 2021. This also allowed temperature from iButtons to be compared with ambient temperatures captured by weather stations. Daily climate maximums were converted to fortnightly means for analyses. Rainfall was recorded as cumulative rainfall per fortnight. The three climate record types covered different temporal ranges. To better compare climate record type, ambient and exposed climate recordings were truncated to match the time period recorded by fissure records (26 June 2020-7 October 2021; 469 days).

Seasonal Phenology
Seasonal phenology was documented at a minimum interval of fortnightly from January 2019 to October 2021. Pinus radiata infested with M. hellenica, since December 2018, at each field site were recorded. Pinus radiata trees were randomly selected for sample collection through a series of dice rolls. P. radiata were excluded from this random selection if their stem diameter was less than 50 cm or if branches were not accessible from groundlevel. Sampled P. radiata at CRP were approximately 75-metres apart and 140-metres apart at DR, meaning the foliage of these trees at each site did not touch. Excised branches were collected from two trees infested with M. hellenica at each site across the study period. Samples were limited to two samples of 30 cm of infested branches per tree per week, with a diameter from 0.5 cm to 2 cm. Infested branches were those that had fresh flocculent or wandering adults (for September to October). Lengths of 9 cm to 11 cm subsections of infested pine twigs and M. hellenica were inspected in the laboratory under a dissection microscope (Wild M3Z, Wild Heerbrugg, Heerbrugg, Canton of St. Gallen, Switzerland) ( Figure 1) and illuminated using a halogen gooseneck lamp (KL 2500 LED, Leica, Wetzlar, Hesse, Germany). Surface area of sampled twigs, calculated from twig length and diameter, was used to determine the density of M. hellenica, expressed as M. hellenica per square centimetre (GPS [cm −2 ]) [49].
The developmental stage of M. hellenica present on host material was determined based on a combination of antennal segment count, presence of body setae or dermal spines, and sclerotized body part size (key to M. hellenica life stages: [28]). Life stages were identified as honeydew producing insects (HPI) via visual confirmation of the presence of honeydew after lightly squeezing their abdomen [50]. The number of HPI was converted to a percentage of the sampled population of M. hellenica (HPI %). To account for combined variation in GPS [cm −2 ] and HPI [%], GPS [cm −2 ] was multiplied by HPI [%] to obtain the density of honeydew producing insects, expressed as HPI [cm −2 ]. This was used as an estimation of the density of M. hellenica producing honeydew, accounting for the decline in GPS [cm −2 ] as each generation progresses [51] and seasonal fluctuations in HPI [%] [50]. Honeydew was collected from HPI from each phenological sample. When less than 10 HPI were available, all HPI were recorded. During the honeydew squeeze test of these HPI subsets, any honeydew was collected using microcapillary tubes of known volume, enabling honeydew volume (HV µL) to be calculated. GPS

Statistical Analysis
Statistical analyses were carried out in SPSS (version 28.0.1.0, IBM Corp., Armonk, NY, USA; access provided by La Trobe University: Department of Ecology, Environment and Evolution; α = 0.05) [52]. Life cycle milestones for recorded generations were reported descriptively, including the first detection of each life stage within generations. The proportions of instar, adult, and egg stages represented in sampled populations were plotted over time, with HPI % overlayed to facilitate visual interpretation of M. hellenica phenology in relation to honeydew availability [50,53] (Dunnett T3). Dunnett T3 is a non-parametric post-hoc test, and thus, does not assume equal distribution or variance among groups [54]. Dunnett T3 was unable to resolve differences in insect density between any month-based groups. To indicate which months were most dissimilar, the Tukey's Honestly Significant Differences (Tukey's HSD) post-hoc test was used to explore differences in density between months [55], although Tukey's HSD outcome must be interpreted with caution (Tukey's HSD assumes data are parametric).  The developmental stage of M. hellenica present on host material was determined based on a combination of antennal segment count, presence of body setae or dermal spines, and sclerotized body part size (key to M. hellenica life stages: [28]). Life stages were identified as honeydew producing insects (HPI) via visual confirmation of the presence of honeydew after lightly squeezing their abdomen [50]. The number of HPI was converted to a percentage of the sampled population of M. hellenica (HPI %). To account for combined variation in GPS [cm −2 ] and HPI [%], GPS [cm −2 ] was multiplied by HPI [%] to obtain the density of honeydew producing insects, expressed as HPI [cm −2 ]. This was used as an estimation of the density of M. hellenica producing honeydew, accounting for the Differences in fortnightly means of M. hellenica activity between site and generation groups were tested via ANOVA. Where ANOVA was used, Tukey's HSD was used for post-hoc comparison of differences among groups. Where data were determined to be non-parametric via Levene's test of equal variances among groups, differences between site and generation groups were tested via KW. Group comparison via pairwise independent samples t-test where KW was used did not assume equal variances between groups. This process was done for GPS [ were ln(x) transformed prior to regression analysis, as this provided stronger climate models that were a better fit for the tested data. Climate variables included temperature ( • C), RH (%), VPD (kPa), as recorded by exposed and bark fissure data loggers, and ambient temperature, cumulative rainfall (mm), and solar exposure (MJ m −2 ). Statistics associated with these models may be found in the Appendix A Tables A9-A12. Indicators of data fit, correlation strength, and statistical clarity for regression models were used to draw comparisons between models: for instance, if the same climate variable would explain variation in insect activity when recorded from different data sources (i.e., bark fissure, exposed trunk, weather station). Models for predicting HPI [cm −2 ] or HV [µL] were plotted for each of the climate variables sources from the data recorder that produced a statistically significant regression model with the strongest correlation strength. The relationship between HPI [cm −2 ] or HV [µL] and solar exposure (MJ m −2 ) was also plotted, although solar exposure (MJ m −2 ) was only sourced from the Bureau of Meteorology. Approximate Nonlinear Durbin-Watson test was used to check for autocorrelation among HPI [cm −2 ] or HV [µL] and climate variables [56,57].
Climate variables explaining a particular M. hellenica activity were included in multiple regression analysis with the goal of obtaining residual values for all climate variables included. Since climate variables were measured with differing units of measurement, standardised residual values were extracted. These standardised residual values were then included in a linear regression analysis to confirm that a given M. hellenica was responsive to the combined effect of the temperature ( • C), RH (%), VPD (kPa), and solar exposure (MJ m −2 ) records that most strongly correlated with M. hellenica activity. These standardised residual plots, with 95% confidence intervals of the mean, are visualized in Figures A1-A4.

Seasonal Phenology
The onset and duration of seasonal feeding behaviour was visualized from January 2019 ( (Table A1), although not for the adult stage. In 2019, the first instars began moulting in February, but the first instars began moulting in late January and early February in 2020 and 2021. The persistence of the second instar appeared to shorten each year; it was present for 17 weeks in 2019, 16 weeks in 2020, and 12 weeks in 2021. The third instars were always initially detected in April. However, the time from the first time that the third instar was detected until 100% of insects sampled were third instars, was less in 2021 than earlier generations. The third instar stage was the longest persisting stage, lasting up to 7 months from April until October. Final ecdysis usually occurred around mid-September when the third instars moulted into adults. By mid-October 2019 and 2020, the third instars were absent.

Feeding Activity
Data collection of honeydew producers in samples did not commence until 21 March 2019; hence, a sudden spike in honeydew availability was visualized in Figure 2A. Insects emerging in November 2019 were not detected as producing honeydew until early January 2020 ( Figure 2B), compared to nymphs emerging November 2020, which commenced honeydew production in early December 2020 ( Figure 2B). HPI [%] was significantly correlated with the proportion of second instar nymphs in a population (Pearson correlation, r = 0.747, t = 7.368, p < 0.001). As a percentage of the sampled population, the second instar stage was found to be most the likely to produce honeydew compared to all other life stages observed (Table A2). The exception to this was a two-week period from late December 2020 to early January 2021 when 86.80% of the first instar nymph population were producing honeydew ( Figure 2C). HPI [%] did not exceed this percentage again in 2021, with 69.35% being the next highest in 2021, recorded from second instar nymphs in early Population progression also appeared to accelerate each year for adults. For populations sampled in mid-September, 8% of M. hellenica were adults in 2019, 17.24% in 2020, and 37.79% in 2021. Two adults were detected in mid-late August 2020, although none were found the next fortnight or in August of 2019 and 2021. Despite most adults being present during October, they may persist into November. Adults oviposit throughout late September and October, with eggs developing for approximately 4-6 weeks before emerging into first instar nymphs around late-November. Adult M. hellenica are semelparous organisms, meaning that they die after a single reproduction event. As a result, by the time the first instars emerge, the adults have almost entirely perished.

Feeding Activity
Data collection of honeydew producers in samples did not commence until 21 March 2019; hence, a sudden spike in honeydew availability was visualized in Figure 2A. Insects emerging in November 2019 were not detected as producing honeydew until early January 2020 ( Figure 2B), compared to nymphs emerging November 2020, which commenced honeydew production in early December 2020 ( Figure 2B). HPI [%] was significantly correlated with the proportion of second instar nymphs in a population (Pearson correlation, r = 0.747, t = 7.368, p < 0.001). As a percentage of the sampled population, the second instar stage was found to be most the likely to produce honeydew compared to all other life stages observed (Table A2). The exception to this was a two-week period from late December 2020 to early January 2021 when 86.80% of the first instar nymph population were producing honeydew ( Figure 2C). HPI [%] did not exceed this percentage again in 2021, with 69.35% being the next highest in 2021, recorded from second instar nymphs in early April. As insects reach their second ecdysis, honeydew production decreased sharply in all generations. Third instar nymphs were much less likely to produce honeydew than second instar nymphs, with third instar HPI [%] never exceeding 50% during winter months. In all generations, there was a small but noticeable rise in HPI [%] from late August to early September of 2019 and 2021.
Analysis of M. hellenica density across 32 months revealed differences across temporally separated populations (Kruskal-Wallis test of independent samples (KW): p = 0.023, F 31,85 = 48.590). Variability in mean density was visualized in Figure 3. Density was highest during November 2020, when M. hellenica first instar nymphs were emerging from eggs. Tukey's HSD provided evidence of density-related differences between specific temporal groups (Table A8), but since this test assumes a normal distribution and equal variances among groups, the outcome must be interpreted with caution. Most months could not be distinguished, with only November 2020 differing from a few months of 2020 and October 2019.
Insects 2023, 14, x FOR PEER REVIEW 9 April. As insects reach their second ecdysis, honeydew production decreased sharp all generations. Third instar nymphs were much less likely to produce honeydew second instar nymphs, with third instar HPI [%] never exceeding 50% during wi months. In all generations, there was a small but noticeable rise in HPI [%] from late gust to early September of 2019 and 2021. Analysis of M. hellenica density across 32 months revealed differences across tem rally separated populations (Kruskal-Wallis test of independent samples (KW): p = 0 F31,85 = 48.590). Variability in mean density was visualized in Figure 3. Density was hig during November 2020, when M. hellenica first instar nymphs were emerging from e Tukey's HSD provided evidence of density-related differences between specific temp groups (Table A8), but since this test assumes a normal distribution and equal varia among groups, the outcome must be interpreted with caution. Most months could no distinguished, with only November 2020 differing from a few months of 2020 and Oct 2019.  Table A3). First and third instar stages dem strated the most substantial difference in honeydew produced per insect (Table A4). instars produced a mean volume of 0.019 μL, and the third instars produced a mea 0.240 μL, which is approximately twelve-times more honeydew. Ln(x) transformed volumes of honeydew produced per insect (HV [µL]) were found to increase with each instar stage (Kruskal-Wallis: p < 0.001, F 2,454 = 189.592, Figure 4). Pairwise comparisons of each instar stage demonstrated significant differences between all group pairings, with HV [µL] production exhibiting an exponential increase as each successive instar stage is reached (Figure 4; Table A3). First and third instar stages demonstrated the most substantial difference in honeydew produced per insect (Table A4). First instars produced a mean volume of 0.019 µL, and the third instars produced a mean of 0.240 µL, which is approximately twelve-times more honeydew.

Insect Activity Response to Climate
Climate variables used to predict honeydew producer insects per cm −2 (HPI [cm −2 ]) generally provided stronger and more statistically significant models compared to either GPS [cm −2 ] Tables A9-A12). AND test did not reveal evidence of autocorrelation between climate variables and ln (HPI [cm −2 ]), or climate variables and ln (HV [µL]), for statistically significant regression models (Tables A11 and A12).

Insect Activity Response to Climate
Climate variables used to predict honeydew producer insects per cm −2 (HPI [cm −2 ]) generally provided stronger and more statistically significant models compared to either GPS [cm −2 ] (Tables A11 and A12).
Statistically significant trends were found clarifying that ln (HPI [cm −2 ]) was at least weakly correlated with most climate variables (Table A11). As temperature increased, HPI [cm −2 ] increased exponentially. Ambient temperature provided the strongest model correlation of this relationship (linear regression: p < 0.001, F1,29 = 35.231, r 2 = 0.549; Figure 6A) and indicated that increases in ambient temperature may explain approximately 53.3% of observed variation in ln (HPI [cm −2 ]). Statistical evidence found that ln (HPI [cm −2 ]) was negatively related to exposed RH, which explained 27.9% of the variation observed in HPI [cm −2 ] (linear regression: p = 0.0013, F1,29 = 12.858, r 2 = 0.303; Figure 6B). Natural log trans- Increases in exposed temperature records were moderately associated with exponential decreases in mean ln (HV [µL]) (linear regression: p < 0.001, F 1,27 = 40.061, r 2 = 0.597; Figure 7A). Here, 58.2% of observed variation in HV [µL] could be explained by variation in exposed temperature. Regression modelling revealed that exponential increases in HV [µL] were moderately correlated with lower RH. This was only the case when the regression model utilised exposed RH records (linear regression: p < 0.001, F 1,27 = 15.360, r 2 = 0.363; Figure 7B), which indicated that approximately 33.9% of variation observed in log transformed HV [µL] could be explained by exposed RH. Evidence was found to show HV [µL] was moderately correlated against exposed VPD (linear regression: p < 0.001, F 1,27 = 31.823, r 2 = 0.541; Figure 7C), with an estimated 52.4% of variation observed in ln (HV [µL]) explained by exposed VPD. Very strong statistical significance was found in support of the observation that ln (HV [µL]) was negatively associated with mean solar exposure (linear regression: p < 0.001, F 1,27 = 18.976, r 2 = 0.380; Figure 7D). This model had a moderately strong correlation, indicating that approximately 36% of observed variation in ln (HV [µL]) may be explained by variation in fortnightly mean solar exposure. This was a weak correlation, with variation in solar exposure explaining approximately 10.3% of observed variation in HPI [cm −2 ]. Figure 6. Relationship between natural log transformed fortnightly density of honeydew excreting nymphs (ln (HPI [cm −2 ])) to (A) exposed temperature, (B) exposed relative humidity, (C) exposed vapor pressure deficit, and (D) solar exposure. Regression model is plotted as continuous black line, and curved lines indicate confidence intervals (95%) of the mean. Strong evidence was found that multiple climate variables correlate with, and explain, observed variation in ln (HPI [cm −2 ]) (Table  A11). Increases in exposed temperature records were moderately associated with exponential decreases in mean ln (HV [μL]) (linear regression: p < 0.001, F1,27 = 40.061, r 2 = 0.597; Figure 7A). Here, 58.2% of observed variation in HV [μL] could be explained by variation in exposed temperature. Regression modelling revealed that exponential increases in HV [μL] were moderately correlated with lower RH. This was only the case when the regression model utilised exposed RH records (linear regression: p < 0.001, F1,27 = 15.360, r 2 = 0.363; Figure 7B), which indicated that approximately 33.9% of variation observed in log transformed HV [μL] could be explained by exposed RH. Evidence was found to show HV [μL] was moderately correlated against exposed VPD (linear regression: p < 0.001, F1,27 = 31.823, r 2 = 0.541; Figure 7C), with an estimated 52.4% of variation observed in ln (HV [μL]) explained by exposed VPD. Very strong statistical significance was found in support of the observation that ln (HV [μL]) was negatively associated with mean solar exposure (linear regression: p < 0.001, F1,27 = 18.976, r 2 = 0.380; Figure 7D). This model had a moderately strong correlation, indicating that approximately 36% of observed variation in ln (HV [μL]) may be explained by variation in fortnightly mean solar exposure. Figure 6. Relationship between natural log transformed fortnightly density of honeydew excreting nymphs (ln (HPI [cm −2 ])) to (A) exposed temperature, (B) exposed relative humidity, (C) exposed vapor pressure deficit, and (D) solar exposure. Regression model is plotted as continuous black line, and curved lines indicate confidence intervals (95%) of the mean. Strong evidence was found that multiple climate variables correlate with, and explain, observed variation in ln (HPI [cm −2 ]) (Table A11). )) to (A) exposed temperature, (B) exposed relative humidity, (C) exposed vapor pressure deficit, and (D) solar exposure. Regression model is plotted as continuous black line, and curved lines indicate confidence intervals (95%) of the mean. Strong evidence was found that multiple climate variables correlate with, and explain, observed variation in ln (HV [μL]) (Table A12).

Discussion
Broad similarities were identified between Mediterranean and Australian M. hellenica populations in the timing of life stages in local seasons [31,53]. As expected, major developmental milestones typically exhibited a six-month delay between Mediterranean and  HV [µL])) to (A) exposed temperature, (B) exposed relative humidity, (C) exposed vapor pressure deficit, and (D) solar exposure. Regression model is plotted as continuous black line, and curved lines indicate confidence intervals (95%) of the mean. Strong evidence was found that multiple climate variables correlate with, and explain, observed variation in ln (HV [µL]) (Table A12).

Discussion
Broad similarities were identified between Mediterranean and Australian M. hellenica populations in the timing of life stages in local seasons [31,53]. As expected, major developmental milestones typically exhibited a six-month delay between Mediterranean and Australian populations. In Australia, the occurrence of life stages between years of study was broadly similar, although the developmental stages showed progressive changes each year; for example, instar stages were detected earlier with successive generations. By contrast, in Greece [50,58], the timing of the occurrence of different instars varied between generations; however, initial detection of each life stage did not occur earlier or later with successive generations. Phenology and development can be delayed or accelerated by changes in climate, a trend observed in other scale insects [10], other Hemiptera [12], and other insect orders [7,53]. Given the relatively short time it has been present, Australian M. hellenica may still be adjusting to the effects of climate drivers on their activity. The current study was limited by the fragmented distribution of exotic P. radiata in Australia. Relating the seasonal phenology of M. hellenica in other climates to that in Australia would provide valuable attestation of the relationships identified here.
The densities of M. hellenica on Australian P. radiata were typically higher than Mediterranean conspecifics feeding on Pinus halepensis during seasonally equivalent months [51]; however, the proportion of honeydew producing insects did not reach proportions as high as those previously documented in Greece [50]. We found evidence that variation in the density of honeydew producing insects and the quantity of honeydew produced could be explained by changes in temperature, humidity, vapour pressure in the atmosphere, and solar exposure. Depending on where climate was recorded (i.e., weather station, exposed tree trunk or bark fissure), the degree of association M. hellenica demonstrated with climate variables varied greatly in this study. The insect's flocculent and preference for crevices and non-sun facing aspects [59] realises a microclimate distinct from broader forest conditions. Microhabitat climate was found to provide weaker predictions of insect activity, as it directly records the climate experienced by sheltered ectothermic M. hellenica [37], rather than the experience of the host. With regards to this study, the climate, as experienced by the host, most often provided the best predictive models.
In Australia, variation in the density of M. hellenica nymphs and proportion of honeydew producing insects (HPI [%]) were associated with warmer and drier conditions. HPI [%] in Australia was often higher over winter than reported in Greece [50]; however, unlike honeydew availability in Greece, there were no instances where 100% of M. hellenica sampled in Australia were excreting honeydew. Broadly, Australian HPI [%] followed similar seasonal patterns to Mediterranean observations, including an increase in honeydew availability following overwintering. Since the phenological study from Greece occurred from 2001 to 2003, it is unclear whether the difference in findings between the current study and previous observations are attributable to underlying geographic differences, inherent generational variation, or host species. The insect density recorded in the current study, from 2019-2020, was only separated by six months from insect densities recorded from Greek conspecifics in 2018-2019 [51]. However, the density of Australian M. hellenica was generally higher than the density of the concurrent Greek conspecifics. Insect density, HPI [%], and HPI [cm −2 ] showed a positive relationship with temperature and atmospheric vapour pressure deficit (VPD) and negative association with rainfall and RH. The volume of honeydew produced per insect had the opposite associations with these same climate variables. These relationships with environmental conditions may explain why M. hellenica in Australia achieved greater densities than Mediterranean populations. Marchalina hellenica in Australia may occupy a more solar exposed, warmer, and drier climate than Mediterranean conspecifics, which are conducive to Pinus growth [60,61] and, by association, tree quality as a host for herbivores [62,63]. Unlike previous studies showing microclimate better reflects insect ontogeny, compared to macroclimate records [37], climate conditions within the microhabitat of settled M. hellenica nymphs in this study were usually not the best indicator of changes in insect activity. This may be due to strong climate effects on the host's activity [64,65], and therefore, the availability of resources for M. hellenica. The influence of host quality on phloem-feeder activity may therefore be amplified because settled M. hellenica are sheltered from direct effects of forest climate condition and depend on climate-driven phloem enrichment [66]. Indeed, phloem-feeder honeydew production and density can fluctuate in association with daily cycles of phloem enrichment [13] and chemical defensive responses [51].
The highest recorded density for Australian M. hellenica in this study was over 5 M. hellenica per cm 2 , compared to the highest density recorded in Greece of approximately 1.4 M. hellenica per cm 2 [51]. Both of these maximum density records coincided with first instar emergence, and at other times, Greek M. hellenica densities were close to zero for most months, and almost always less than 0.25 M. hellenica per cm 2 . By contrast, Australian density presented here was 0.5 per cm 2 or higher across most months. One explanation for elevated M. hellenica density in Australia compared to Greece may be that the insect has escaped its natural enemies. The notable enemy is the predatory fly, Neoleucopis kartliana, a monophagous predator and prospective biocontrol agent [53,67]. Neoleucopis kartliana is recorded as being the most effective and widely distributed predator of M. hellenica in Turkey [32]. With N. kartliana absent, M. hellenica populations in Australia may proliferate with relatively little impact from natural enemies [68][69][70]. If N. kartliana is deemed suitable for introduction into Australia as a classical biocontrol agent, inoculative mass-releases of this fly may be guided by the seasonal patterns of Australian M. hellenica documented here. For example, the effectiveness of biological control agents can be heavily dependent on the presence of particular life stages, not just any individual, from a given target species [71]. This may be critical to any future release of N. kartliana, as the timing of its phenology matches the seasonal phenology of M. hellenica [53].
Host species is also likely important for M. hellenica growth and survival. In the Mediterranean, Pinus halepensis is the primary host for M. hellenica, a tree that shares a long history of hosting M. hellenica, which provide honeydew as an alternative to floral nectar for apicultural purposes [17]. This relationship dates back to the late Roman and early Byzantine Empires along the east coast of the Mediterranean Sea [19,29,72]. A long co-evolutionary history may also be implied by P. halepensis possessing targeted defence responses that coincide with seasonal nymph emergence [51,73]. With M. hellenica recorded as acquiring the novel host in 2014 [16], P. radiata in Australia may be lacking specific defence adaptations that have evolved in P. halepensis. Even if P. radiata has a similar seasonal-specific defence, the shifting phenology of M. hellenica in Australia may allow the insect to escape a temporally restricted defence by emerging earlier. A similar process is documented in other insects (e.g., aphids, bark beetles, Heteroptera), notably involving interactive effects with regional climate [74][75][76] or broadscale climate change [5,42,77].
This study's findings are similar to previous investigations linking climate to changes in both plant and herbivore activity [76,[78][79][80][81]. Climate effects may therefore affect the quality of M. hellenica as a host for predatory insects, as has been observed in parasitic/predatory Diptera [82,83]. Marchalina hellenica occurs in areas of similar climate conditions to its endemic Greek distribution [67] and is largely dependent on P. radiata phloem feeding-the quality of which is influenced by climate conditions [74,84]. For example, VPD is known to be a driver of daily cycles of transpiration and photosynthesis in Pinus [40,85,86] and responses to novel climate conditions may alter phloem quality for herbivores. In Pinus, these responses are a result of tight regulation of stomatal pore closure and, therefore, strict responses in transpiration, photosynthetic activity, or enrichment of phloem elements [87][88][89]. With such rapid responses to environmental changes, M. hellenica, by expropriating carbohydrates, amino acids, and other macromolecules [24], may illicit these strict responses [90] and potentially exacerbate stressful climate effects on the host tree. As global climate conditions deviate from previous means, this insect may experience changes in phenological timing and feeding activity [5,7,91]. Investigation of how M. hellenica life cycle and feeding activity may respond to changes to environmental conditions in its endemic geographic range or on native hosts would provide valuable comparisons with the relationships identified in Melbourne populations.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.
Acknowledgments: Thank you to Breanna Butera, Hannah Young, Loretta Gibilisco, and Sam Heywood for support and during initial field site assessments, collections, and sample inspections. Thank you to James Martin who assisted with field the bulk of field work during Melbourne's COVID-19 lockdown. Thank you to Greg Lefoe and Umar Lubanga for providing feedback on late-stage drafts. Thanks to Angus Carnegie for inviting M.J.S. to propose this PhD project.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript; or in the decision to publish the results.
Appendix A Table A1. Summary statistics of first detections of notable life history milestones and feeding activity across Marchalina hellenica instars of a given generation. ± represents one standard deviation.