Vector-borne diseases are one of the most significant public health burdens globally, with 80% of the total world population at risk [1] . Arboviruses, including dengue, Zika and chikungunya, are of particular concern due to the recent increase in global cases promoted by the rapid spread of both, their primary mosquito vector Aedes aegypti as well as their secondary vector Aedes albopictus [2] . Dengue infection can be asymptomatic but clinical presentations range from mild dengue fever (DF), a febrile illness similar to influenza, to the severe forms of dengue; dengue shock syndrome (DSS) and dengue haemorrhagic fever (DHF) [3] . Most Zika (ZIKV) infections are asymptomatic, with only approximately 20% of infections causing symptoms [4, 5] . The clinical presentations of symptomatic ZIKV can include Zika fever, congenital Zika syndrome and Guillain-Barré syndrome. Congenital Zika syndrome refers to a group of birth defects, notably Microcephaly, which have been associated with ZIKV during pregnancy [6, 7] . Infection with the chikungunya virus (CHIKV) is characterised by sudden onset fever, rash and arthralgia [8] .
Dengue causes an estimated 390 million infections per year and has a distribution that covers every continent of the world with the exception of Antarctica [9] . The number of global dengue cases reported to the WHO has increased 15 fold over the last 20 years, with deaths also seeing a significant increase (4-fold) [ challenging with a diurnal biting mosquito. A compound effect is the development of insecticide resistance in populations of Ae. aegypti reported in areas of Colombia [27] .
This study aims to investigate the epidemiology of these three arboviruses (dengue, Zika and chikungunya) co-circulating in a single vector species (Ae. aegypti) in three distinct ecosystems in Colombia between 2007 -2017. This is achieved using a multifactorial approach considering the potential correlations of several meteorological and socio-economic factors with disease incidence. We show that specific climatological factors are strong drivers for these arboviral diseases to which contextual socio-economical characteristics can act as modifiers. Importantly, we find a discriminatory pattern between these three diseases highlighting unexpected dynamics of transmission between Zika and dengue particularly in an area of high dengue circulation.
Three municipalities in Colombia were selected for analysis: Bello, Cúcuta and Moniquirá (Fig 1) . Cases of dengue, severe dengue (dengue shock syndrome and dengue haemorrhagic fever), chikungunya and Zika reported in each epidemiological week were obtained for the period of 2007 to 2017 from SIVIGILIA (National Public Health Surveillance System, Colombian National Institute of Health) [28] . This period of study (2007) (2008) (2009) (2010) (2011) (2012) (2013) (2014) (2015) (2016) (2017) [31] . The MPI at municipality level was obtained from DANE using data collected in the 2018 National Population and Housing Census and using the indicators and respective weightings listed in S2 Table [32] .
Patterns of disease incidence by location Differences in the total burden of all three Ae. aegypti borne viruses as well as the individual burden of dengue and severe dengue were investigated using the total number of cases from 2007-2017. Poisson Generalised Linear Models (GLMs) were initially carried out, revealing overdispersion (data variance greater than expected for the given model) statistics of 545, 542 and 123 for total disease, dengue and severe dengue respectively. To correct for the large overdispersion the GLMs were recalculated with negative binomial errors [33] using the glm.nb function of the R package MASS [34] . Differences between all three municipalities were tested with Tukey pairwise comparisons using the glht function of the multcomp R package [35] . Incidence of chikungunya in 2015 was initially modelled using a Poisson GLM which revealed an overdispersion statistic of 2.9. The standard errors were therefore corrected using quasi-GLMs where the variance was theta x mu. Where mu was the mean of the dependant variable and theta the dispersion parameter of the quasi-GLMs [33] . Quasi-GLMs were conducted using the glm function from the R package stats [36] . Zika incidence was modelled for the year 2017 only, initially a Poisson GLM was used and an overdispersion statistic of 99 was detected. As the overdispersion statistic was above 20 it was corrected for using negative binomial errors [33] . Total Ae. aegypti borne disease and dengue incidence were also modelled for 2015 and 2016 using quasi-GLMs to account for low level overdispersion except for total disease in 2016 which had a dispersion statistic of 90 and was therefore modelled with a negative binomial GLM. Population was used as an offset in all models in order to standardise disease incidence by population size.
For the pattern of disease over time we used total yearly incidence data. Poisson GLMs were initially used to model each disease in each municipality but overdispersion was detected in some models, hence error distributions were adjusted accordingly. For Bello incidence of both total disease (dengue, severe dengue, chikungunya and Zika) and dengue alone were modelled using quasipoisson GLMs, correcting for overdispersions of 3.69 and 3.63 respectively. For Cúcuta a negative binomial GLM was required for total disease incidence in order to correct for overdispersion of 32.52 and quasi-GLM was used for dengue incidence due to slight overdispersion of 4.87. Severe dengue incidence in Cúcuta was not found to be significantly overdispersed when modelled with a Poisson GLM (overdispersion statistic = 1.79), as the overdispersion statistic was <2. For Moniquirá both total disease incidence and dengue incidence alone were modelled using quasi-GLMs, correcting for respective overdispersion statistics of 2.39 and 2.34 [33] . Population was used as an offset in all models in order to standardise disease incidence by population size. All quasi-GLMs were conducted using the glm function from the R package stats [36] and negative binomial GLMs used the glm.nb function of the R package MASS [34] . Differences between the years were tested with Tukey pairwise comparisons using the glht function of the multcomp R package [35] .
The correlations between climatic variables and the total disease incidence (dengue, severe dengue, chikungunya and Zika) across all three locations were investigated using a generalised additive model (GAM). The weekly disease incidence and weather data for each municipality was converted into 4-week data, matching the dates of epidemiological months.
Combining the data into 4-week periods rather than individual weeks reduced zero inflation improving the reliability of the GAM outputs. All climate variables were lagged by plausible time lags for their effect on disease incidence, of 4 and 8 weeks. Square root transformations were used for total disease incidence and each weather variable due to non-normal distribution. Generalised additive models were chosen due to their ability to model non-linear relationships between a response variable (disease incidence) and multiple explanatory variables (climate variables) [37] . A quasi-maximum likelihood Poisson GAM was used in order to prevent possible overdispersion [33] . Population size was used as an offset to standardise disease incidence by population. Initially all climate variables with both 4 and 8-
week time lags were assumed to have a non-linear relationship and were therefore modelled as smoothed terms. Subsequent analysis of the effective degrees of freedom (edf) was used to identify variables with edf = 1.0, suggesting linearity. These variables were then included in the model as linear rather than smoothed terms. Generalised cross validation (GCV) was used to determine the most appropriate model. Generalised additive modelling and subsequent model validation was conducted using the R package mgcv [38] . Visualisation of GAM estimations were conducted using the mgcViz R package [39] .
Principle components analysis (PCA) was used for dimensional reduction to allow the inclusion of socio-economic data with previously compiled geographic and climate data. The
The three municipalities Bello, Cúcuta and Moniquirá were chosen due to their geographical separation, and distinct climate characteristics, demographics and burden of Ae. aegypti borne diseases (Table 1) . The spread of the data for the period analysed was assessed in 12-week periods, using smooth moving averages (SMA) (Fig 3) . 
Total disease incidence over the 11-year period was significantly lower in Bello than in Cúcuta (p = < 0.001) and Moniquirá (p = 0.005) ( Fig 4A) for all three diseases. The number of dengue cases were similarly high between Moniquirá and Cúcuta (p = 0.99). However, severe dengue incidence was significantly different across all three municipalities: Cúcuta had the highest burden of severe dengue when compared to both Bello (p = <0.001) and
Moniquirá (p = 0.005), and Moniquirá had a significantly higher burden of severe dengue when compared to Bello (p = 0.047) ( Fig 4A) .
The incidence data were also analysed separately for the years when the outbreaks of chikungunya and Zika were reported, 2015 and 2016, respectively (Fig 4B and 4C ). This allowed for a more directly and meaningful comparison of the burden represented by these three diseases. Cases of chikungunya in 2015 were not significantly different between any of the municipalities (Fig 4B) . However, Cúcuta had significantly higher incidence of Zika than is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint
The copyright holder for this this version posted November 20, 2020. ; https://doi.org/10.1101/2020.11.18.20233940 doi: medRxiv preprint both Bello (p= <0.001) and Moniquirá (p = <0.001) ( Fig 4C) . Interestingly, in the same year of the Zika outbreak (2016) each municipality had a significantly different number of dengue cases. Cúcuta had the lowest incidence of dengue, and Moniquirá the highest (Fig 4C) . While Cúcuta had the highest incidence of dengue in previous years, in 2016 the same location experienced the lowest incidence of dengue accompanied by the highest incidence of Zika ( Fig 4C) .
We compared the number of cases per 100,000 people in each year for dengue, severe Following the significantly high dengue incidence in Bello and Moniquirá in 2016, the incidence stabilised, and the incidence reported in 2017 were not statistically different to that of years prior to 2016. This was not the case however in Cúcuta where incidence of dengue continued to fall, with the lowest incidence of the study period observed in 2017 (Fig 5B) . In 2017 incidence of Zika was also much lower with only 3 cases per 100,000 people reported in Cúcuta in 2017 (Fig 5A and 5C ). is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint
The copyright holder for this this version posted November 20, 2020. ; https://doi.org/10.1101/2020.11.18.20233940 doi: medRxiv preprint
The geographical settings for the three locations studied here Bello, Cúcuta and Moniquirá (Table 1) We explored potential relationships between the time series climate data and total disease incidence in Bello, Cúcuta and Moniquirá at shorter 4 and 8-week time lags using a generalised additive model (GAM) as explained in Methods. The estimates from the quasipoisson GAM explained 57.6% of the variance in total disease incidence over time (Table 2 ). Effective degrees of freedom (edf) close to 1 represent relation close to linearity while high edf values for the smooth terms suggest that the relationship between climatic variables and disease incidence is non-linear. The GAM identified significant relationships between disease incidence and precipitation at 4 and 8-week lags, average humidity (4-week lag), minimum temperature (4 and 8-week lags), average temperature (8-week lag), maximum wind speed (4 and 8-week lag) and average wind speed (4-week lag) ( Table 2) . is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint
The copyright holder for this this version posted November 20, 2020. ; Table   2 were further investigated by plotting the smoothed variance of the latter against the ranges covered for each climate variable (Fig 7) . The analysis of the three most significant climatic contributors to total disease incidence -temperature, wind speed, and precipitationdelineated clear trends on how climate affected disease transmission. Increasing Tmin from 8 o C to 16 o C, at either 4 or 8-week, reduced the total disease incidence while the average temperatures (Tavg) up to 24 o C contributed to incremental levels of disease (Fig 7A-C) . This could be an indication of the need for a range of low minimal temperatures to adjust the average values at which disease transmission is more efficient. is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint
The copyright holder for this this version posted November 20, 2020. ; https://doi.org/10.1101/2020.11.18.20233940 doi: medRxiv preprint Increasing wind speed above 1 m/s was associated by an increase in disease incidence with a peak around 2 m/s after which disease incidence declined until around 3 m/s where a small rise in disease risk can also be seen. This relationship between maximum wind speed and disease incidence was observed after both 4 and 8-week time lags (Fig 7D-E) . The decisive influence on wind speed was substantiated by the negative effect on disease incidence at incremental average wind speeds (Fig 7F) . Precipitation showed a positive relationship with disease incidence above 90 mm at 4-week time lags and above 25 mm for the data with 8week time lag (Fig 7G-H) . The high level of non-linearity shown for the relationship between average humidity and disease incidence (edf = 7.5) ( Table 2) is detailed in Fig 7I. Disease incidence increased as humidity increased between 55-60%. This was followed by a slight decrease and plateau between 60-70%, a more rapid increase was observed between 70 and 75% above which disease incidence begins to decline (Fig 7I) .
We followed a holistic approach by further including socio-economic data for these three locations in the investigation for modifiers to the disease transmission of dengue, Zika and chikungunya. The overall multidimensional poverty index (see Methods) was lowest in Bello at 14.2. Cúcuta and Moniquirá had similar MPIs at 25.7 and 27.1 respectively (Table 3) . The overall MPIs in Cúcuta and Moniquirá were similar. However, there were differences in specific poverty measures relevant to the transmission of vector borne viral diseases. Cúcuta had higher rates of overcrowding (16.4%), barriers to both childhood and youth services (2.2%) and healthcare services (5.2%) and inappropriate exterior wall material (6.1%) ( Table   3 ). The findings also pointed to other socio-economic indexes that also affect general health and well-being in Moniquirá: inadequate excreta disposal (sanitation) (12.1%), no access to an improved water source (20.4%), illiteracy (15.3%) and low education achievement (59.4%) were all highest in Moniquirá (Table 3) .
In order to introduce the socio-economic data into the analyses undertaken with the epidemiological and climatic data we carried out a dimensionality reduction and correlation with a principal component analysis (PCA). Importantly there was a clear separation of the three municipalities along both dimensions PC1 and PC2 that together integrates 89.9% of the compiled parameters (Fig 8) . This approach made apparent a discriminatory set of factors both climatic and socio-economic for all three locations. Cúcuta had an extensive combination of climate factors (i.e. wind speed and temperature) that together with school absence, dependency, overcrowding, wall material and school failure are potential modifiers is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint
The copyright holder for this this version posted November 20, 2020. ; https://doi.org/10.1101/2020.11.18.20233940 doi: medRxiv preprint of dengue, Zika and chikungunya risk in this municipality. Interestingly, Moniquirá showed mainly socio-economic variables (i.e. water source, sanitation, illiteracy, low education, flooring material, child labour, high MPI, informal work) to be potential modifiers for disease risk. On the other hand, Bello had mainly climate variables as potential modifiers of disease transmission -average humidity, precipitation and elevation -with only health insurance as a socio-economic factor. is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint
The copyright holder for this this version posted November 20, 2020. ; https://doi.org/10.1101/2020.11.18.20233940 doi: medRxiv preprint
This study set out to determine the longitudinal dynamics of three Aedes arboviral diseases co-circulating in three regions of Colombia over a 11-year period between 2007 and 2017.
We found significant differences in the burden of viruses among the three municipalities studied. Bello had the lowest level of disease incidence across all diseases. Cúcuta had the highest incidence of severe dengue (2007-2017) and Zika (2016) and the highest overall disease incidence. In addition to climatic factors the burden of these vector borne diseases in is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint
The copyright holder for this this version posted November 20, 2020. ; https://doi.org/10.1101/2020.11.18.20233940 doi: medRxiv preprint cases across the Americas was lower than any of the 10 previous years [54] with a 73% decline between 2016 and 2017 alone [53] .
Changes in epidemiological surveillance systems can cause the identification of inaccurate patterns of disease incidence. However, we did not observe any indication of significant changes to the surveillance system used to report Aedes borne viruses to SIVIGILIA, the database used here. However, the circulation of multiple viruses in the same localities at the same time does provide challenges for surveillance systems. Clinical presentation of Zika is very similar to that of dengue [22, 55] and this can cause cases to be misidentified when laboratory testing is not conducted. We note that although the cases analysed in this study are all confirmed cases, confirmation is not always done by laboratory testing but also by epidemiological links. Misidentification could therefore explain the increase in dengue that was observed in Bello and Moniquirá in 2016, where low incidence of Zika was reported. Coinfection of the primary vector Ae. aegypti with multiple arboviruses (i.e. DENV, CHIKV, ZIKV) has been reported following laboratory exposure [56] [57] [58] [59] [60] [61] . Aedes mosquitoes have also been shown capable of transmitting more than one arbovirus in a single biting event [57, 61, 62] . Although coinfection has yet to be found in wild Ae. aegypti [62] . Coinfection of multiple Aedes borne viruses has been reported in mammalian hosts including humans [63] [64] [65] [66] [67] [68] [69] . In Colombia patients have been diagnosed with DENV-CHIKV [48, 70, 71] , DENV-ZIKV [48, 71] , CHIKV-ZIKV [48, 71] and DENV-CHIKV-ZIKV [22, 49] coinfections.
However, the frequency of DENV-ZIKV co-infections seems low at 8.8% [71] . DENV and ZIKV co-transmission in mice through the bite of Ae. aegypti mosquitoes showed preferential transmission of ZIKV [61] . is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint
The copyright holder for this this version posted November 20, 2020. ; https://doi.org/10.1101/2020.11.18.20233940 doi: medRxiv preprint Host cross-immunity of ZIKV and DENV could have been a contributing factor in the dengue declines observed in this study. The observed decline in dengue cases following Zika outbreaks reported within this study and in others across the Americas suggest that there may be cross-immunity between ZIKV and DENV in humans [51] [52] [53] . Flavivirus immunity involves a T cell response and studies have reported cross-reactivity of CD4+ and CD8+ T cells to both DENV and ZIKV [99] [100] [101] [102] [103] . Cross-reactivity of antibodies and T cells and cross-immunity from Zika, although not necessarily conferring cross-protection, has been presented as the most probable reason explaining the decline of dengue across the Americas in 2017 [53].
Specific climatic factors associated here significantly affected disease incidence. This was particularly evident for Cúcuta, the municipality with the highest disease burden. We found significant co-relationship between average temperature and wind speed with disease transmission, with a peak at around 2m/s, consistent with findings of previous studies [104] [105] [106] [107] . Ae. aegypti has a small flight range of 200 m and high wind speeds reduce mosquito flight distances while low winds mean reduced dispersion of mosquitoes. We found the optimum wind speed to be around 2 m/s which is in line with current knowledge of mosquito flight [104, 108, 109] . Climate variables can be used to build predictive models to anticipate when outbreaks of dengue, chikungunya and Zika are likely to occur [110] [111] [112] [113] . This is useful in the prioritisation of vector-control resources. Recent modelling studies have reported an increase in the incidence and geographical spread Ae. aegypti borne viruses when using climate change simulation models [114] [115] [116] [117] [118] . This highlights the importance of consideration of environmental factors when assessing risk of vector-borne disease [119] .
We report differences in measurements of socio-economics between Bello, Cúcuta and Moniquirá. Bello, the municipality with the lowest burden of Aedes borne viruses also had the lowest poverty index, whilst Cúcuta and Moniquirá were much higher in both disease is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint
The copyright holder for this this version posted November 20, 2020. ; https://doi.org/10.1101/2020.11.18.20233940 doi: medRxiv preprint incidence and multidimensional poverty. Higher incidence of Aedes borne disease has been associated with lower socio-economic status and higher poverty levels [120] [121] [122] [123] [124] [125] [126] [127] . Cúcuta had the highest rate of critical overcrowding. Overcrowding has been reported to be an important contributing factor to dengue incidence [119] [120] [121] 128] . Inadequate sanitation, and no access to improved sources of water are both well-known contributing factors in increasing burden of Aedes borne disease due to the ecology of Aedes mosquitoes [26, 119, [123] [124] [125] 129] . These socio-economic risk factors were highest in Moniquirá where there were also high levels of low educational achievement and illiteracy. Illiteracy and low educational level have previously been associated with increased vulnerability to dengue in Colombia and Brazil [130, 131] .
Heading into the COVID-19 2020 global pandemic, 2019 saw the highest number of dengue cases in The Americas for twenty years [132] . The factors that can exasperate the spread of COVID-19 (e.g. poor sanitation, lack of access to clean water, overcrowding) are also contributing factors to arboviral diseases. Disease surveillance and control programs are being significantly impacted by the pandemic, and the consequences of this will not be clear for some time but the pandemic also places an extra burden on already fragile health systems.
Lockdown and social distancing measures will be detrimental for Aedes control which often involves community focused activities (e.g. door to door) led by community health workers and vector control experts. The World Mosquito Programme for example halted all control activities involving community interaction (e.g. releasing Wolbachia positive Ae. aegypti) for three months. Social distancing can be difficult or impossible in poorer areas with high population densities. It is therefore more imperative than ever that we more fully understand the dynamics of these important diseases which are likely to escalate over the coming years.
Having different ecosystems Bello, Cúcuta and Moniquirá presented a valuable opportunity to explore longitudinal arboviral disease incidence over 11 years that encompassed a Zika is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint
The copyright holder for this this version posted November 20, 2020. ; epidemic. Chikungunya was the only disease for which incidence did not significantly differ between the three municipalities. Cúcuta had the greatest disease incidence having the most favourable climatic factors and greater poverty index but as it borders with Venezuela mass movement of people is also suggested to be a contributing factor. Climatic factors associated with disease incidence were precipitation, average humidity, temperature and wind speed.
Co-transmission of dengue and Zika during the epidemic led to a significant reduction of dengue cases in Cúcuta where dengue had previously been high. This significant finding warrants further investigation. Where the poverty index was low, as in Bello, so was the disease incidence. Socio-economic factors such as barriers to health and childhood services, inadequate sanitation, poor housing and poor water supply were implicated as drivers of disease transmission. Aedes aegypti and Ae. albopictus are increasing their geographical range and climate change is predicted to alter the distribution of these vectors and hence disease risk. Arboviral epidemiology is further complicated by humanitarian crises (e.g.
Venezuela) and the COVID-19 pandemic which reinforces the urgency for understanding the dynamics of these global health problems. Context dependent and actionable understanding of the drivers for disease transmission that consider local dynamics, both climatic and socioeconomic, should contribute to the design of more effective vector mosquito control programmes. Socio-economic variables: no access to improved water source (Water_Source), inadequate disposal for excreta (Sanitation), illiteracy, low educational achievement (Low_Ed), inappropriate flooring material (Flooring_Material), child labour, multidimensional poverty index (MPI), informal work, school absence, dependency rate (Dependency), barriers to health services (Barriers_HS), critical overcrowding, inappropriate wall material (Wall_Material), barriers to early childhood services (Bariers_CS), no health insurance (Health_Insurance). Climate variables: elevation, maximum wind speed (WSmax), average wind speed (WSavg), maximum temperature (Tmax), minimum wind speed (WSmin), average temperature (Tavg), precipitation (Pre) and average humidity (Havg). The length of the arrows represents the contribution of each variable.
Supporting information S1 is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint
The copyright holder for this this version posted November 20, 2020. is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint
The copyright holder for this this version posted November 20, 2020. is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint
The copyright holder for this this version posted November 20, 2020. is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint
The copyright holder for this this version posted November 20, 2020. ; https://doi.org/10.1101/2020.11.18.20233940 doi: medRxiv preprint 