Comparing sterile male releases and other methods for integrated control of the tiger mosquito in temperate and tropical climates

The expansion of mosquito species worldwide is creating a powerful network for the spread of arboviruses. In addition to the destruction of breeding sites (prevention) and mass trapping, methods based on the sterile insect technique (SIT), the autodissemination of pyriproxyfen (ADT), and a fusion of elements from both of these known as boosted SIT (BSIT), are being developed to meet the urgent need for effective vector control. However, the comparative potential of these methods has yet to be explored in different environments. This is needed to propose and integrate informed guidelines into sustainable mosquito management plans. We extended a weather-dependent model of Aedes albopictus population dynamics to assess the effectiveness of these different vector control methods, alone or in combination, in a tropical (Reunion island, southwest Indian Ocean) and a temperate (Montpellier area, southern France) climate. Our results confirm the potential efficiency of SIT in temperate climates when performed early in the year (mid-March for northern hemisphere). In such a climate, the timing of the vector control action was the key factor in its success. In tropical climates, the potential of the combination of methods becomes more relevant. BSIT and the combination of ADT with SIT were twice as effective compared to the use of SIT alone.


B Undesirable increase in female abundance
When started late in the mosquito season, SIT with late releases of less than 1100 males per hectare, has been found to cause a temporary increase in female abundance (see Results), which is not a desirable effect of control. This undesirable increase of female abundance has been reported in several studies [2][3][4] . We further investigated the causes by running two sets of three simulations for SIT, BSIT and a control, one with larval competition and the other without (i.e., removing the competition term 1 + P k P (1−a r prev ) , eq. (2)). To allow comparison while avoiding the burst in population induced by the removal of larval competition, the total numbers of adult males and females were caped at 200 individuals.
An increase in female abundance was observed for SIT, which disappeared with the removal of larval competition from the model (Figure 2), suggesting that the sterile males releases induce a reduced larval competition, leading temporarily to more emergence. This increase is not observed with BSIT, either with or without larval competition, probably because the pp-contamination of breeding sites heavily lower down the larval population and overcome the effect of any competition process.

C Model validation in the temperate climate
The simulated population dynamic without any sterile male releases (eq. (2) : λ M s = λ M sc = 0 and a = 0) has been validated on entomological data recorded by the "Entente Interdépartementale de Démoustication (EID) Méditerranée", the public agency in charge of the mosquito control in the southeast France.
Ovitraps were placed in 5 residential areas (Table 1) located in Montpellier (Vert-bois, Pompigane Sud and Aiguelongue) or in the adjacent city of Castelnau-Le-Lez (Volhe and Chevalerie) ( Figure 3). Studied areas were described as discontinuous medium to dense urban fabrics by the reference CORINE, an European reference describing the land-cover structure based on Figure 1. Best release starting date in the temperate climate for a maximum SIT and BSIT effect. (A) Effectiveness at different starting dates for SIT (blue) and BSIT (red). Periods reaching at least 95% and 99% of the maximum effect are respectively filled in blue and red. (B) SIT and BSIT effect on female dynamics over time. The release starts at its optimal day (March, the 24 th ) in year 2015. The scenario consists of 18 releases of 1,000 males/ha during a period of 4 months. Grey background represents the mosquitoes favorable period, i.e.,. when no diapause occurs. satellite imagery 5 . Ovitraps were mostly placed in sites shaded by vegetation. They consisted of 3L plastic black buckets filled with 2L of tap water, in which a floating polystyrene square (5 x 5cm) acted as an ovipositing support. Biolarvicide Bacillus thuringiensis israelensis was added into the bucket to prevent the production of mosquitoes from the trap.
The number of trapped eggs were monitored weekly between April and November, in 2015, 2016 and 2017 (Table 1). The trapping network in each parcel was composed of 8 to 25 ovitraps. We considered that a reasonable-enough trapping effort was performed to get an accurate description of the local Ae. albopictus population dynamics. Hatched and unhatched eggs were removed from the trap each week for identification and count at the laboratory. These observed data were then compared to the predictions of the model using a Spearman correlation test. The predicted numbers of new eggs produced per week were independently computed as (γ gc (β n F n + β p F p )). Model simulations were highly consistent with the dynamics recorded in the field, with a significant correlation of 61% (p value 10 −27 ). However, the model tends to underestimate the number of new eggs produced at the beginning of the favorable period (April). An early egg peak is observed in 2017 in all studied areas ( Figure 4). This burst is not captured by our simulations and could 2/8 Figure 2. Population dynamics during late implementation of SIT or BSIT, with or without pupal competition. The release period is filled in grey. We bounded the adult females abundance to avoid exponential growth of population when the density-dependent (i.e., the competition) process are removed. not be explained by temperature or rainfall fluctuations. It is probably responsible for the lower correlation (52%, p value 0.01) observed in the Aiguelongue area where the few records (24) provided were all performed in 2017. Eggs measured in Chevalerie, Volhe and Vert-bois areas, were followed for more than one year and show higher correlation coefficients: 0.75 (p value 10 −11 ), 0.69 (p value 10 −9 ) and 0.65 (p value 10 −11 ), respectively. The best fit to the entomological data was observed in the South Pompignane neighborhood with a significant correlation of 86% (p value 10 −9 ); this area is the only one where no records were performed in year 2017.

D Full model, modified from Haramboure et al. 1
The original model was developed for Reunion Island (southwest Indian Ocean, tropical climate) by Haramboure et al. 1 . The present model is an adaptation that allows for the temperate climate, using Montpellier (southern France, Mediterranean 3/8 Figure 3. Study area. Colored points represent the two meteorological stations (blue for Prades-Le-Lez and red for Montpellier's airport) of the area. The same color code was used to link each parcel with its closest meteorological station. Entomological data were provided for 5 parcels : Vert-Bois, Aiguelongue, South Pompignade, Chevalerie and Volhe (filled parcels). The maps were produced using QGIS 3.10 (http://www.qgis.org). The shapefile of the regions of France is taken from the open street map (https: //www.data.gouv.fr/fr/datasets/contours-des-regions-francaises-sur-openstreetmap/). temperate climate, i.e., hot and dry summers vs. mild and humid winters) as a reference: Constant parameters are in Greek letters. They include mortality rates (µ X being the mortality rate of compartment X, 11 compartments in total, see text), and the transition rates from emerging to nulliparous females (γ F em ), and from nulliparous to parous females (γ gc ). γ gc is then multiplied by the number of eggs laid per ovipositing nulliparous (β n ) and parous (β p ) females to compute the egg-laying rate. An additional mortality rate is added for female compartments, due to their host-seeking behaviour (µ r ). σ is the sex ratio at emergence. Each male being able to fertilize several females, we assumed that the mating probability for females was 1, i.e., independent of the male density.
Latin letters are weather-dependent functions that act on the mortality rates (m x ), the transition rates ( f x ) and the environmentdefined carrying capacities (k x ) for the aquatics stages. Carrying capacities drive the density-dependent mortality at larval stage m L 1 + L k L and the pupae density-dependent success for adult emergence exp −µ F em 1 + P k P . Temperature impacts m x and f x , whereas k x is rainfall-dependent.
During diapause (d start ,d end ), observed only in the temperate climate, the transition between the eggs and larvae compartments is stopped by the z parameter that represents the egg dormancy process 6 .
When SIT, or BSIT, starts (T start ), a quantity of sterile males (λ X , with X = M s or M sc , resp.), is periodically released into the population (∆ t ) during a given period of time τ. The sterile males have a probability to mate with an emerging females (F em ) depending of their relative availability compared to the overall number of males, i.e., their proportion adjusted by the competition parameter ω: p = ωX ωX+M . No eggs produced by sterile females will hatch. Emerging females that encounter a sterile male thus become sterile (F s ), while the other become nulliparous (F n , probability 1 − p).
Moreover, for BSIT the released sterile males (M sc ) are also coated with pyriproxyfen (PP). Males attempt to mate with any females of the population (F tot = F em + F n + F p + F s ), even if already fertilized. Each encounter of a contaminated male with a female makes it loose a part of its PP coating until total decontamination after κ F matings : ω(M sc +M s )+M (F n + F p + F s )γ gc , but only to the non-contaminated fraction 1 − B c B tot , with B c the number of contaminated sites among all breeding sites (B tot ). Females loose part of the contaminant at each oviposition, until total decontamination after κ Bc gonotrophic cycles. Pupae emergence rate is affected by the probability of the pupae to 1) grow in a contaminated breeding site and 2) to survive PP with a probability φ 1 − B c B tot (1 − φ ) . Finally, breeding site become decontaminated at a rate ν. Bold sections of the equation (eq. (2)) model the vector control methods that are not based on sterile males releases : ADT (O T , S T ), ovitraps (O T ), BGS-traps (BGS) and prevention (r prev ) through breeding site destruction.

E Model parameters
Parameters and weather-dependent function have been adjusted according to expert knowledge and literature (Table 2 and 3).
The density of breeding sites B tot was estimated from the urbanization level in each delimited subarea, as the mean number of breeding sites per household (10 − 20) multiplied by the number of households per hectare (25 − 30) and the surface of the parcel. All were estimated from field observations, as well as the distinction between permanent and rainfall-dependent breeding sites.

F Sensitivity analysis in temperate climate
A sensitivity analysis (SA) was performed on the model outputs to identify key parameters of the SIT and BSIT in temperate climate. Four model outputs were studied : the reduction rate, the resilience, the reduction rate of BSIT compared SIT and the undesirable female augmentation. SIT parameters (λ M s , ω, µ M s ), BSIT parameters (λ M sc , ν, κ Bc , κ F , µ M sc ) and parameters describing the release period (T start , τ, ∆ t ) were varied one by one through 100 trajectories to complete a Morris Sensitivity analysis 11 with the sensitivity R package. The grid of parameters space was discretized by 10 levels and the walking step was defined as 10 2 , following Morris's recommendations.
Each parameter space (Table 1) was determined from literature or expert knowledge, except for T start , which was just estimated. Table 3. Functions adjusted for temperate conditions 6 . R is the rainfall and T, the temperature at time t.

Function Definition
Expression f E Transition function from egg to larvae (T (t) − T E )/T DD E f L Transition function from larvae to pupae −0007T 2 + 0.0392T − 0.3911 f P Transition function from pupae to emerging adult 0.0008T 2 − 0.0051T + 0.0319 m L Larva mortality µ L + 0.0007exp(0.1838(T − 10)) m P Pupa mortality µ P + 0.0003exp(0.2228(T − 10)) m F Female mortality µ F + 0.0003exp(0.1745(T − 10)) R norm Normalized weekly rainfall amount ∑ t i=t−14 R(i) /100 k L Environmental carrying capacity for larvae k L f ix + min (k L var R norm (t), k L var ) k P Environmental carrying capacity for pupae k P f ix + min (k P var R norm (t), k P var ) As previously found for the tropical climate 1 , the main parameters influencing the mosquito population dynamics are essentially the same between SIT and BSIT ( Figure 5). The elementary effect variance increases for most of the parameters according to the absolute value of the elementary effect. Most parameters have few or no influence on the output values. Figure 5. Keys parameters of the SIT and BSIT in temperate climate. Sensitivity analysis performed on the SIT (blue triangles) and the BSIT (red points). Four outputs were studied: the resilience (A), the reduction rate (B), the reduction rate of the BSIT compared to the SIT (C), and the increase in the female population (D). The variance (σ ) of each elementary effect is plotted against its absolute value (µ * ). Parameters without influence (σ = 0 and µ * = 0) were removed. The bottom-left region (low µ * and σ ) of the reduction rate (B) and the female population increase (D) plots have been magnified (black boxes).
For 3 outputs out of 4, the date for the start of the releases (T start ) is by far the most influencing factor ( Figure 5 (B, C, D)). Its elementary effect absolute value, describing its overall influence on the outputs, is 2 to 100 times higher than the effect of other parameters. Its participation to the output has either a non-linear relationship with the value, or a strong interaction with other inputs (high σ ). The length of the release (τ) and the sterile male mortality (µ s ) are parameters that can significantly 7/8 modify the population resilience after SIT or BSIT ( Figure 5 (A)).
Even if the starting day of release remains by far the major influential parameter, BSIT is less sensitive than SIT to this parameter for the reduction rate and resilience outputs ( Figure 5 (A, B)). The reduction rate of BSIT compared to SIT is mostly influenced by the starting day of releases, and, to a much lesser extent, by sterile males competitiveness (ω) and the number of sterile males released (λ M s ) ( Figure 5 (C)). Male competitiveness, starting date of releases and the quantity of sterile males released drive an undesirable augmentation of the females number in SIT, while none has been observed for BSIT ( Figure 5  (D)). These parameters may act in interaction (high σ ). Boost related parameters (e.g. the decontamination rate, ν, the larval resistance to PP, φ and the PP transfers, κ Bc and κ F ) have few or no influence in the outputs ( Figure 5).