Rio de Janeiro city, Brazil, is a dengue and chikungunya endemic city with complex, heterogeneous and disordered urban structure, in which high-income areas and slums are in vicinity to each other. In this scenario, Ae. aegypti larval habitats have been shown to vary according to urban characteristics such as the availability of piped water, human density and housing type (Maciel-de-Freitas et al., 2007b; David et al., 2009). Field collections were conducted in four neighborhoods: Curicica (CUR), Prainha (PRA), Tubiacanga (TUB) and Vila Valqueire (VILA), representing the range of socio-demographic and economic conditions across Rio de Janeiro (Table 1). Therefore, it is expected that the containers sampled are representative of Ae. aegypti larval habitats in Rio. At each site, agents from the Rio de Janeiro Department of Health randomly visited and inspected the houses until they found approximately 10 negatives and 50 positives water-holding containers for Ae. aegypti immature forms (larvae and/or pupae). Dengue and chikungunya epidemic season in Rio usually goes from December-April. Data were collected April–August 2010.

Containers were classified according to type, location (indoors or outdoors) and the following parameters were registered in loco: dissolved O₂, water temperature (oxymeter SevenGo Pro, Mettler Toledo, Ohio, EUA), pH (pHmeter model 826 pH Mobile, Metrohm, São Paulo, Brazil), conductivity (codutivimeter SevenGo, Mettler Toledo, Ohio, EUA), air temperature, wind speed, light incidence (CA810 Luxmeter, Chauvin Arnoux, Paris, France) and the presence of insect predators (e.g., Toxorhynchites and Odonata immature stages). Water volume was calculated from the length, width and height of the water column in large containers (>10 L) or directly measured by transferring it to a graduated beaker. The surface of larval habitats was calculated using container size measurements.

A sample of the water and sediment from each container was collected using an algae sampler, which consisted of a kitchen sink plunger with a brush inside connected to an plastic syringe for sample collection (Loeb, 1981). The sediment was sampled by pressing the plunger firmly to the substrate, which was scrubbed with the brush. The contents of the plunger cup were then collected with a syringe. Two sizes of algae samplers were used according to containers surface: one using an 20ml syringe (7 cm² of sampling area) for containers with <50cm² and one with an 80ml syringe (44.2 cm² of sampling area) for containers >50cm² of surface area. The water/sediment samples from containers were transferred to plastic bags, properly identified and transported in a cool box to the laboratory, where they remained frozen (–20°C) for about 2 months.

Unfiltered frozen samples (10 ml) of water/sediment were processed for Total Phosphorus via orthophosphate determination in a spectrophotometer according to Valderrama (1981) and Total Nitrogen through nitrate detection with ion chromatography with UV detection (Metrohm, São Paulo, Brazil).

Water samples were filtered with GF/F 47 mm Whatman filters, preserved with 15 μL of phosphoric acid and kept in the fridge at 4°C for no longer than 10 days. DOC was determined by an elemental carbon analyzer (model HiperToc, Thermo Scientific, Waltham, EUA) for oxidation with sodium persulfate and determination with UV light. Samples were analyzed in analytical triplicate.

The sediment was analyzed for total organic carbon (sediment organic carbon, SOC) levels after drying in porcelain crucibles in an oven with forced air circulation for 30 min at 60°C. The remaining dry material was weighed, and the organic matter was dosed by calcination in a muffle furnace at 400°C for 12 h.

All larvae and pupae were collected with a strainer and plastic pipette, transferred to a plastic bags containing water from the larval site and stored in a cool box. For ethical issues, all the remaining water was discarded from containers. When larval site elimination was not possible (e.g., water for domestic use), the water holding container was covered or treated with larvicide by local health agents.

Samples were transported to the laboratory where immature mosquitoes were counted, and the pupae reared to adults in dechlorinated tap water at 28 ± 2°C and 70 ± 10% relative humidity. After taxonomic identification (Consoli and de Oliveira, 1994), adult Ae. aegypti mosquitoes were killed with ethyl acetate and dried in a heat oven at 50°C by 15 h. Subsequently, specimens were weighed in a precision scale (Denver Instrument, New York, EUA) and the wing length measured by taking the distance from the axillary incision to the apical margin excluding the fringe (Harbach and Knight, 1980).

Statistical analysis was conducted with R environment (R Development Core Team, 2011). Container types were classified into five categories according to the criteria used by the Brazilian Ministry of Health (Ministério da Saúde, 2009; Table 1). A list of all positive and negative containers types sampled per neighborhood is available in Supplementary Tables S1, S2. Exploratory analyses were first performed to identify environmental differences between neighborhoods, as well as variations in physicochemical characteristics according to larval habitat categories.

Principal Component Analysis (PCA) and Random Forest (RF) classification (Liaw and Wiener, 2002) were employed to test whether larval habitats exhibited specific physicochemical profiles according to container category or geographical origin. Differences between positive and negative larval habitats for Ae. aegypti were also graphically evaluated through PCA. RF classifying models were trained with a random sample of ∼70% of the data (N = 139 larval habitats) and then validated with the remaining 30% (N = 60 larval habitats) using the “caret” R package. The overall model accuracy and the Kappa index, as well as the identification of most influencing parameters for classification, were obtained after applying the RF models to the validation dataset. PCA was applied on the correlation matrix after scaling variables to have unit variance since they were on different scales and therefore exhibit different variances. In this case, variables were standardized to mean equal to zero and variance equal to one. The environmental parameters “air temperature,” “container surface,” and “location (indoors or outdoors)” were removed from data analysis due to strong correlation (r > 0.7) with “water temperature,” “water volume,” and “light incidence,” respectively. The presence of predators was excluded from the analyses since only a single predatory larva of Toxorhynchites sp. was observed in one bromeliad from VILA. The environmental variables considered in RF and PCA analyses are listed in Table 3.

In order to elucidate possible associations between larval habitats characteristics and mosquito abundances, the total count of immature per recipient was included as dependent variable in generalized linear mixed models (GLMMs) using the “lme4” R package (Bates et al., 2015). We considered the total number of Culicidae immature forms collected, since only the pupae were reared to the adult stage for taxonomic identification, and also positive and negative containers for Ae. aegypti. The independent variables considered in the GLMMs are listed in Table 3. In order to improve the quality of model fit, the variables “volume,” “light incidence,” “conductivity,” “DOC,” and “SOC” were included in the models at the log scale. The count of Culicidae immature per larval habitat was analyzed with a GLMM with negative binomial distribution with the neighborhood as random effect in order to incorporate the variation between collection sites in the models. A regression with negative binomial distribution was preferred over the traditional Poisson distribution because data exhibited over-dispersion (i.e., variance was larger than the mean), confirmed by the Pearson Chi² test and Dispersion Statistic >1, calculated using the “msme” R package (Hilbe and Robinson, 2013). The consistency of data with the negative binomial distribution was verified using the goodness of fit test “Minimum Chi-squared” (Pearson χ² = 240.7, df = 233, p = 0.35) from the “goodfit” command, implemented in the “vcd” R package (Meyer et al., 2020). After adjusting independent variables to univariate models, those with significant effects (α = 0.1) on immature Culicidae per larval habitat were used to build GLMMs. The most informative and parcimonious model was selected through second-order Akaike’s information criterion scores (AICc) and level of support (AICcWt), calculated using the “AICcmodavg” package (Mazerolle, 2019). Collinearity between independent variables was checked in the best model through VIFs (Variance Inflation Factors) (Zuur et al., 2010). The assumptions of the best model were examined by checking heteroscedasticity, residuals dispersion and the presence of outliers with the R package DHARMa (Hartig, 2020). The conditional Pseudo-R-squared for GLMMs (R²c), calculated using the MuMIn R package (Bartón, 2014), was employed as goodness-of-fit metric and a measurement of variance explained by the entire model, including both fixed and random effects.

Regarding the body size analyses, Ae. aegypti average wing length and dry weight were calculated for each larval habitat according to mosquito sex. Since there was a significant, but moderate correlation between these two variables (Pearson’s correlation female: t = 11.7, p < 0.001, r = 0.69, r² = 0.48; male: t = 10.4, p < 0.001, r = 0.60, r² = 0.36; Supplementary Figure S1), both were fitted as dependent variables at the Linear Mixed Models (LLMs) with the neighborhood as random effect (Table 3). As data was considered normally distributed, a gaussian distribution was applied (One-sample Kolmogorov-Smirnov test p > 0.05). The independent variable “wind speed” was removed from these analyses once aquatic mosquito larvae and pupae are submerged, therefore would not be directly exposed to variations in this environmental condition. Model selection and validation followed the previously mentioned procedures.

The dry weight and wing size were compared between Ae. aegypti originating from larval habitats with and without Ae. albopictus (or any other mosquito species). Differences in mosquito body size were tested between groups using the Kruskal-Wallis test as male wing size was not considered normally distributed (One-sample Kolmogorov-Smirnov test p < 0.05).

Larval surveys are routinely performed by health agents from the Rio de Janeiro City Health Department (Secretaria Municipal de Saúde do Rio de Janeiro, SMS-RJ). Access to containers inside dwellings were obtained after consent of householder approving the entrance of the research team and the health agent who routinely worked at that neighborhood.

In total, 240 water holding containers were selected at the four field sites, of which 204 were positive for Ae. aegypti immature presence. The number of sampled containers according to category and study site are shown in Table 4 and a total of 3,184 mosquito immatures were collected, of which 1,168 were pupae. Pupae survival to the adult stage allowed the taxonomic identification of 1,064 mosquitoes, 692 males and 372 females. Although most of adult mosquitoes were identified as Ae. aegypti (Table 5), Ae. albopictus, Aedes (Ochlerotatus) sp., Culex sp., and Wyeomyia sp. were also identified in the inspected containers (Table 5). Of the 981 adult Ae. aegypti, 638 were males and 343 were females which provided the measure of body weight and wing length as a proxy for body size.

The accuracy and Kappa index of RF models fitted to classify larval habitats using physicochemical parameters of water and sediment were higher according to geographical location than container category (Table 6). Overall model accuracy for neighborhood was 78.3%, with Classification sensitivity (true positive rate) varying between 63 (VILA) to 92% (TUB). Water temperature, volume and light incidence were the most important predictors for sample discrimination into neighborhood of collection. Median water temperature in TUB containers was 25.7°C meanwhile it ranged between 21.4 and 21.6°C in the other field sites. Container volume tended to be higher in TUB and contrasted to the median volume of larval habitats sampled elsewhere (TUB: median of 7.2 L; other sites: median of 0.3–0.6 L). On the other hand, median light incidence was 63.5 lx in TUB larval habitats while it ranged between 113 and 154.1 lx in the other sampled neighborhoods. A less efficient classification was observed for container categories, with the container volume as the most important predictor and overall accuracy of 50%. Raw larval habitat data can be found in Supplementary Table S3.

According to the PCA, the first two axes explained 27.3% of total variation. The first axis accounted for 13.8% of data variance, with DOC, light incidence and water temperature as the variables that most contributed to explain the dataset variation (Supplementary Figure S2A). Dissolved O2, SOC and water temperature were the main variables explaining the second component (accounting for 13.5% of variance) (Supplementary Figure S2B). TUB larval habitat formed a more heterogeneous cluster, compared to CUR, PRA and VILA (Figure 1A). No discrimination was noted between container types (Figure 1B) nor between positive and negative Ae. aegypti larval habitats (Figure 1C), corroborating RF classification results. Considering the higher heterogeneity between TUB larval habitats in comparison to the other field sites, GLMMs and LMMs for abundance and body size investigations were adjusted with neighborhood as a random effect.

Univariate GLMM analyses indicated a significant association of immature Culicidae abundance in larval habitats with volume, total nitrogen, total phosphorous and DOC concentrations (Supplementary Table S4). Subsequently, GLMMs were fitted using these four independent variables and compared to the null model. The most informative GLMM had volume, total nitrogen, total phosphorous and DOC as independent variables, but only volume and DOC exhibited significant effects on the number of immature Culicidae. Model R²c was 0.37. These results suggest that the number of immature Culicidae increased with larval habitat volume and DOC concentration (Table 7). For example, TUB, the neighborhood in which container volume tended to be higher, exhibited 31.3 immature per larval habitat in average, while it varied between 4.7 and 8.8 in the other field sites. Descriptive analysis of immature abundance data can be found in Supplementary Table S5.

Female Ae. aegypti originated from 126 larval habitats. Regarding average dry weight and wing length per container, the univariate LMMs presented only water temperature as significant independent variable, with R²c of 0.15 and 0.06, respectively. The models output indicates that both female weight and wing length increased with water temperature (Table 8). Additionally, neighborhood variance estimates were near zero, suggesting there was virtually no among-neighborhood variance left to be explained by the random intercepts (Table 8). Indeed, in TUB, the neighborhood with higher water temperatures, female Ae. aegypti had slightly larger wings, with 2.56 mm in average, while it ranged between 2.39 and 2.51 in the other field sites. However, the same pattern could not be clearly seen for the dry weight, as PRA and TUB mosquitoes had in average 0.32 mg, while it was 0.38–0.39 mg in average in VILA and CUR, respectively (Supplementary Table S5).

Male Ae. aegypti originated from 160 larval habitats. The univariate LLMs indicated a significant association between average male dry weight per container with dissolved O2, conductivity, phosphorous and DOC (Supplementary Table S6). Therefore, GLMMs were fitted using these four independent variables and were compared to the null model. The full model was the one with the highest support, with male dry weight increasing with water conductivity (Table 8), which was in median 322.50 μS/cm (ranging from 2.25 to 1,374 μS/cm), with no variation between neighborhoods or container categories, as shown by RF and PCA analysis. The R²c was 0.09 for this LLM. For male wing length, only water temperature had a significant effect, suggesting an increase of wing length with temperature. As seen for female Ae. aegypti, male tended to have larger wings in TUB (2.06 mm in average against 1.90–1.98 mm in average in the other neighborhoods, Supplementary Table S5). However, the R²c was quite low (0.06), as shown for the other LMMs fitted for body size measurements. Neighborhood variance estimates were near zero, suggesting there was virtually no among-neighborhood variance left to be explained by the random intercepts (Table 8).

The two species were simultaneously observed in 12 out of 204 larval habitats (5.9%), six from TUB and in two containers from each of the other field sites. A total of 30 male and 22 female Ae. aegypti mosquitoes were collected from shared larval habitats with Ae. albopictus, while 545 males and 301 Ae. aegypti females originated from larval habitats where this species occurred exclusively. Due to disparities higher than 10-fold between the number of Ae. aegypti originating from shared and non-shared containers, a random sample of 30 female and 30 male mosquitoes was taken from the second group in order to balance sample sizes for statistical comparison. Median dry weight and wing length tended to be lower for both male and female when co-occurring with Ae. albopictus, although with no statistical significance (Figures 2A–D: Female wing length: KW χ² = 1.88, df = 1, p = 0.17; Female dry weight: KW χ² = 1.93, df = 1, p = 0.16; Male wing length: KW χ² = 1.99, df = 1, p = 0.16; Male dry weight: KW χ² = 2.76, df = 1, p = 0.10).