---
authors:
  - givenNames:
      - Marie
      - C
    familyNames:
      - Russell
    type: Person
    emails:
      - marie.clare.russell@gmail.com
    affiliations:
      - name: >-
          Department of Life Sciences, Imperial College London, Silwood Park
          Campus
        address:
          addressCountry: United Kingdom
          addressLocality: Ascot
          type: PostalAddress
        type: Organization
  - givenNames:
      - Catherine
      - M
    familyNames:
      - Herzog
    type: Person
    affiliations:
      - name: Center for Infectious Disease Dynamics, Pennsylvania State University
        address:
          addressCountry: United States
          addressLocality: University Park
          type: PostalAddress
        type: Organization
  - givenNames:
      - Zachary
    familyNames:
      - Gajewski
    type: Person
    affiliations:
      - name: >-
          Department of Biological Sciences, Virginia Polytechnic Institute and
          State University
        address:
          addressCountry: United States
          addressLocality: Blacksburg
          type: PostalAddress
        type: Organization
  - givenNames:
      - Chloe
    familyNames:
      - Ramsay
    type: Person
    affiliations:
      - name: Department of Biological Sciences, University of Notre Dame
        address:
          addressCountry: United States
          addressLocality: Notre Dame
          type: PostalAddress
        type: Organization
  - givenNames:
      - Fadoua
    familyNames:
      - El
      - Moustaid
    type: Person
    affiliations:
      - name: >-
          Department of Biological Sciences, Virginia Polytechnic Institute and
          State University
        address:
          addressCountry: United States
          addressLocality: Blacksburg
          type: PostalAddress
        type: Organization
  - givenNames:
      - Michelle
      - V
    familyNames:
      - Evans
    type: Person
    affiliations:
      - name: >-
          Odum School of Ecology & Center for Ecology of Infectious Diseases,
          University of Georgia
        address:
          addressCountry: United States
          addressLocality: Athens
          type: PostalAddress
        type: Organization
      - name: MIVEGEC, IRD, CNRS, Université Montpellier
        address:
          addressCountry: France
          addressLocality: Montpellier
          type: PostalAddress
        type: Organization
  - givenNames:
      - Trishna
    familyNames:
      - Desai
    type: Person
    affiliations:
      - name: Nuffield Department of Population Health, University of Oxford
        address:
          addressCountry: United Kingdom
          addressLocality: Oxford
          type: PostalAddress
        type: Organization
  - givenNames:
      - Nicole
      - L
    familyNames:
      - Gottdenker
    type: Person
    affiliations:
      - name: >-
          Odum School of Ecology & Center for Ecology of Infectious Diseases,
          University of Georgia
        address:
          addressCountry: United States
          addressLocality: Athens
          type: PostalAddress
        type: Organization
      - name: >-
          Department of Veterinary Pathology, University of Georgia College of
          Veterinary Medicine
        address:
          addressCountry: United States
          addressLocality: Athens
          type: PostalAddress
        type: Organization
  - givenNames:
      - Sara
      - L
    familyNames:
      - Hermann
    type: Person
    affiliations:
      - name: Department of Entomology, Pennsylvania State University
        address:
          addressCountry: United States
          addressLocality: University Park
          type: PostalAddress
        type: Organization
  - givenNames:
      - Alison
      - G
    familyNames:
      - Power
    type: Person
    affiliations:
      - name: Department of Ecology & Evolutionary Biology, Cornell University
        address:
          addressCountry: United States
          addressLocality: Ithaca
          type: PostalAddress
        type: Organization
  - givenNames:
      - Andrew
      - C
    familyNames:
      - McCall
    type: Person
    affiliations:
      - name: Biology Department, Denison University
        address:
          addressCountry: United States
          addressLocality: Granville
          type: PostalAddress
        type: Organization
editors:
  - givenNames:
      - Sarah
      - E
    familyNames:
      - Cobey
    type: Person
    affiliations:
      - name: University of Chicago
        address:
          addressCountry: United States
          type: PostalAddress
        type: Organization
datePublished:
  value: '2022-01-19'
  type: Date
dateReceived:
  value: '2021-06-21'
  type: Date
dateAccepted:
  value: '2021-12-01'
  type: Date
title: >-
  Both consumptive and non-consumptive effects of predators impact mosquito
  populations and have implications for disease transmission
description:
  - >-
    Predator-prey interactions influence prey traits through both consumptive
    and non-consumptive effects, and variation in these traits can shape
    vector-borne disease dynamics. Meta-analysis methods were employed to
    generate predation effect sizes by different categories of predators and
    mosquito prey. This analysis showed that multiple families of aquatic
    predators are effective in consumptively reducing mosquito survival, and
    that the survival of 
  - type: Emphasis
    content:
      - Aedes
  - ', '
  - type: Emphasis
    content:
      - Anopheles
  - ', and '
  - type: Emphasis
    content:
      - Culex
  - ' mosquitoes is negatively impacted by consumptive effects of predators. Mosquito larval size was found to play a more important role in explaining the heterogeneity of consumptive effects from predators than mosquito genus. Mosquito survival and body size were reduced by non-consumptive effects of predators, but development time was not significantly impacted. In addition, '
  - type: Emphasis
    content:
      - Culex
  - ' vectors demonstrated predator avoidance behavior during oviposition. The results of this meta-analysis suggest that predators limit disease transmission by reducing both vector survival and vector size, and that associations between drought and human West Nile virus cases could be driven by the vector behavior of predator avoidance during oviposition. These findings are likely to be useful to infectious disease modelers who rely on vector traits as predictors of transmission.'
isPartOf:
  volumeNumber: 11
  isPartOf:
    title: eLife
    issns:
      - 2050-084X
    identifiers:
      - name: nlm-ta
        propertyID: https://registry.identifiers.org/registry/nlm-ta
        value: elife
        type: PropertyValue
      - name: publisher-id
        propertyID: https://registry.identifiers.org/registry/publisher-id
        value: eLife
        type: PropertyValue
    publisher:
      name: eLife Sciences Publications, Ltd
      type: Organization
    type: Periodical
  type: PublicationVolume
licenses:
  - url: http://creativecommons.org/licenses/by/4.0/
    content:
      - content:
          - 'This article is distributed under the terms of the '
          - content:
              - Creative Commons Attribution License
            target: http://creativecommons.org/licenses/by/4.0/
            type: Link
          - >-
            , which permits unrestricted use and redistribution provided that
            the original author and source are credited.
        type: Paragraph
    type: CreativeWork
keywords:
  - vector ecology
  - predation
  - disease dynamics
  - meta-analysis
  - Mosquito
identifiers:
  - name: publisher-id
    propertyID: https://registry.identifiers.org/registry/publisher-id
    value: 71503
    type: PropertyValue
  - name: doi
    propertyID: https://registry.identifiers.org/registry/doi
    value: 10.7554/eLife.71503
    type: PropertyValue
  - name: elocation-id
    propertyID: https://registry.identifiers.org/registry/elocation-id
    value: e71503
    type: PropertyValue
fundedBy:
  - identifiers:
      - value: 1R01AI122284-01
        type: PropertyValue
    funders:
      - name: National Institutes of Health
        type: Organization
    type: MonetaryGrant
  - identifiers:
      - value: BB/N013573/1
        type: PropertyValue
    funders:
      - name: Biotechnology and Biological Sciences Research Council
        type: Organization
    type: MonetaryGrant
  - identifiers:
      - value: President's PhD Scholarship
        type: PropertyValue
    funders:
      - name: Imperial College London
        type: Organization
    type: MonetaryGrant
about:
  - name: Ecology
    type: DefinedTerm
  - name: Epidemiology and Global Health
    type: DefinedTerm
genre:
  - Research Article
bibliography: article-71503.references.bib
---

# Introduction

While it is well known that predation reduces vector populations through consumptive effects, non-consumptive effects of predators can also greatly impact prey demographics [@bib106]. Mosquitoes are vectors of a variety of debilitating and deadly diseases, including malaria, lymphatic filariasis, and arboviruses, such as chikungunya, Zika, and dengue [@bib147; @bib148]. Consequently, there is motivation from a public health perspective to better understand the different drivers of variation in mosquito traits that can ultimately impact vector population growth and disease transmission. In addition, recent work has suggested that incorporation of vector trait variation into disease models can improve the reliability of their predictions [@bib23]. In this study, systematic review and meta-analysis methods are used to synthesize a clearer understanding of the consumptive and non-consumptive effects of predators on mosquito traits, including survival, oviposition, development, and size.

Mosquito insecticide resistance is recognized as a growing problem [@bib54; @bib57; @bib82] leading some to suggest that control efforts should rely more heavily on ‘non-insecticide based strategies’ [@bib14]. The consumptive effects of predators on mosquitoes have previously been harnessed for biocontrol purposes. Past biocontrol efforts have used predators such as cyclopoid copepods [@bib69; @bib87; @bib120; @bib142] and mosquitofish [@bib107; @bib121] to target the mosquito’s aquatic larval stage. The strength of the consumptive effects of these predators on mosquitoes can be influenced by multiple factors, including predator-prey size ratio and temperature. Predator-prey body size ratios tend to be higher in freshwater habitats than other types of habitats [@bib18], and attack rate tends to increase with temperature [@bib68; @bib29], although other studies suggest a unimodal response to temperature [@bib139; @bib41].

Predators can also have non-consumptive effects on prey [@bib103], and these effects are thought to be more pronounced in aquatic ecosystems than in terrestrial ecosystems [@bib106]. Non-consumptive effects of predators are the result of the prey initiating anti-predator behavioral and/or physiological trait changes that can aid in predator avoidance [@bib59; @bib80]. Such plasticity in certain prey traits may also result in energetic costs [@bib81]. Predator detection is key for these trait changes to occur and can be mediated by chemical, tactile, and visual cues [@bib58]. In mosquitoes, exposure to predators is known to affect a variety of traits including behavior, size, development, and survival [@bib10; @bib17; @bib114; @bib119; @bib151]. Experimental observations of predator effects on mosquito size and development are inconsistent and results sometimes vary by mosquito sex. For example, exposure to predation was found to increase the size of _Culex pipiens_ mosquitoes [@bib2] but decrease the size of _Culiseta longiareolata_!number(0)[@bib128]. In addition, female _Aedes triseriatus_ exhibited shorter development times when exposed to predation at high nutrient availability [@bib100], but male _C. longiareolata_ had longer development times in the presence of predators [@bib128]. In some cases, a shared evolutionary history between predator and prey organisms can strengthen the non-consumptive effects of predators on mosquitoes [@bib21; @bib124].

This investigation assesses the consumptive and non-consumptive effects of predators on mosquito traits and describes how these effects could impact disease transmission. The roles of vector genus, predator family, mosquito larval instar (an indicator of prey size), and temperature are also examined as potential moderators of predator effects. Non-consumptive effects of predators are expected to cause a smaller reduction in mosquito survival than consumptive effects because, in practice, measures of consumptive effects always include both consumptive and non-consumptive effects. Based on previous findings, larger predators are more likely to consumptively reduce mosquito survival [@bib74]. In addition, _Aedes_ mosquito larvae may be more vulnerable to consumption than other genera because of the high degree of motility observed in this genus [@bib35; @bib89; @bib126]. The oviposition response to predation is expected to be weakest among _Aedes_ species that oviposit above the water line, due in part to their delayed-hatching eggs [@bib146]. Predation is predicted to reduce mosquito size and lengthen development time, consistent with the reduced growth response observed in other insect systems [@bib59]. Certain non-consumptive effects of predation, particularly oviposition site selection and decreased vector size, are likely to play important roles in the dynamics of mosquito-borne disease.

# Materials and methods

## Literature screening

A systematic search was conducted for studies on predation of mosquitoes that were published between 1970 and July 1, 2019 using both PubMed and Web of Science search engines, according to the PRISMA protocol [@bib93]. Mosquito vectors of the _Anopheles_ and _Aedes_ genera were specifically highlighted in our search terms because these genera contain the vector species that transmit malaria, yellow fever, and dengue – the three most deadly mosquito-borne diseases worldwide [@bib62]. Searches included 18 combinations of three vector predation terms (mosquito predat\*, _Anopheles_ predat\*, _Aedes_ predat\*) and six trait terms (survival, mortality, development, fecundity, dispers\*, host preference). Abstracts from the 1136 studies were each screened by two different co-authors, using the ‘metagear’ package in R [@bib76; @bib108]. If either screener thought the study had information relevant to predation of mosquitoes, or both screeners thought the abstract was ambiguous, the study was read in full. This resulted in 306 studies that were fully reviewed to determine if any predation data could be extracted ([Figure 1](#fig1)).

figure: Figure 1.
:::
![](article-71503.rmd.media/fig1.jpg)

### Flowchart demonstrating the literature search, screening process, data exclusions, and the resulting seven different vector trait data subsets.
:::
{#fig1}

## Study exclusion criteria

Data were extracted from studies that collected data on non-consumptive and/or consumptive effects of predators on mosquitoes. Studies were required to have a mean, error measurement, and at least two replicates for both control and predator treatments. The control treatment was required to have all the same conditions as the predator treatment, such as prey density and type of water, without the predators. Studies that were not published in English and studies that did not differentiate between predators of multiple families were excluded. Studies were also excluded if oviposition by free-flying female mosquitoes could have interfered with observing the consumptive effects of predators on vector survival. The final database comprised data extracted from 60 studies ([Supplementary file 1](#supp1)). The data included observations from laboratory experiments, as well as semi-field experiments, in which mesocosms of different treatments were observed in outdoor settings.

## Data extraction

Variables related to the publication, the vector, the predator, and the effect size ([Table 1](#table1)) were extracted from each study. Data from tables and text were recorded as they were published, and data from figures were extracted using WebPlotDigitizer [@bib116]. Error measurements that were not originally presented as standard deviations were converted to standard deviations prior to the effect size calculation.

table: Table 1.
:::
### Variables extracted from included studies.

| Variable                                     | Description                                                                                                                                          |
| -------------------------------------------- | ---------------------------------------------------------------------------------------------------------------------------------------------------- |
| Publication data:                            |                                                                                                                                                      |
| Title                                        | Full study title                                                                                                                                     |
| Journal                                      | Name of journal that published the study                                                                                                             |
| Year                                         | Year of publication                                                                                                                                  |
| Study environment                            | Environment where the experiment took place: lab or semi-field                                                                                       |
|                                              |                                                                                                                                                      |
| Vector data:                                 |                                                                                                                                                      |
| Order, Family, Genus, Species                | Taxonomic identification                                                                                                                             |
| Trait                                        | Outcome that was measured (e.g. survival, development, etc.)                                                                                         |
| Stage                                        | Life stage: egg, larva, pupa, or adult                                                                                                               |
| Larval instar                                | Early (1st and 2nd instars), late (3rd and 4th instars), both, or NA (eggs, pupae, or adults)                                                        |
| Sex                                          | Male or female                                                                                                                                       |
|                                              |                                                                                                                                                      |
| Predator data:                               |                                                                                                                                                      |
| Phylum, Class, Order, Family, Genus, Species | Taxonomic identification                                                                                                                             |
| Starved                                      | Whether the predator was starved: yes or no                                                                                                          |
| Time starved                                 | Amount of time that the predator was starved (in minutes)                                                                                            |
| Predation effect                             | Consumptive or non-consumptive                                                                                                                       |
|                                              |                                                                                                                                                      |
| Effect size data:                            |                                                                                                                                                      |
| Units                                        | Units of extracted data                                                                                                                              |
| Control mean                                 | Average of the outcome measured among the controls                                                                                                   |
| Control standard deviation                   | Standard deviation of the outcome measured in the controls                                                                                           |
| Control number of replicates                 | Number of control replicates                                                                                                                         |
| Predation mean                               | Average of the outcome measured in the predator treatment                                                                                            |
| Predation standard deviation                 | Standard deviation of the outcome measured in the predator treatment                                                                                 |
| Predation number of replicates               | Number of predation replicates                                                                                                                       |
| Experiment ID                                | Alphabetic assignment to mark observations sharing a control group or representing the same prey individuals as originating from the same experiment |
|                                              |                                                                                                                                                      |
| Additional data:                             |                                                                                                                                                      |
| Experiment time (days)                       | Duration of the experiment in days                                                                                                                   |
| Data source                                  | Graph or text                                                                                                                                        |
| Number of predators                          | Number of predators with access to prey, or ‘cue’ if there are no predators with direct access to prey                                               |
| Number of prey (vectors)                     | Number of mosquito prey that are exposed to predation                                                                                                |
| Arena volume (mL)                            | Volume of the arena where prey encounter predators                                                                                                   |
| Time exposed to predator(s)                  | Amount of time (in days) when the predator has direct access to the mosquito prey                                                                    |
| Temperature (°C)                             | Temperature during the predation interaction                                                                                                         |
| Type of predator cue                         | Predator cues, or cues from both predator(s) and dying conspecifics; NA for observations with a consumptive predation effect                         |
:::
{#table1}

## Data exclusions

A PRISMA plot of literature inclusion and exclusion is provided in [Figure 1](#fig1). Observations where insecticide was used were excluded because insecticides are known to interfere with consumptive and non-consumptive effects of predators [@bib32; @bib63]. In addition, observations from experiments with mosquito prey of two or more species were excluded because it was not possible to account for effects from apparent competition or prey-switching. Observations of vector fecundity, vector competence, behavioral traits other than oviposition, as well as observations where the vector trait was marked as ‘other’ were not analyzed because each of these traits were only recorded from three or fewer studies.

Due to protandry, the earlier emergence of males to maximize their reproductive success, mosquitoes respond to sex-specific selective forces that influence their development time and body size [@bib71]. Under low resource conditions, female mosquitoes are likely to maximize body mass by extending their development time, whereas males tend to minimize their development time at the expense of lower body mass [@bib71]. Observations of mosquito development time and body size in our database that were not sex-specific were excluded so that these vector traits could be analyzed while controlling for sex. In addition, among the observations of development time and body size, some predator means did not necessarily represent an evenly weighted average of the replicates. For example, if a total of 20 mosquitoes from three different predator replicates survived to adulthood, the mean development time and size of those 20 individuals may have been reported. To represent an evenly weighted average of the replicates, it is necessary to first calculate summary statistics among multiple individuals that emerge from the same replicate, and then report the average of the replicate-specific means. Observations that might have been influenced by uneven representation of replicates were excluded to prevent pseudo-replication from altering later meta-analyses.

For consumptive observations where life stage-specific survival was reported after more than 10 days of predator exposure, only data on survival marked by adult emergence were included for analysis. Effects observed among immature vector stages after such a long period of predator exposure were not analyzed because they could have resulted from a combination of non-consumptive effects on development, and consumptive effects on survival. Development time observations that were reported as the inverse of development time (units of days^–1^) were excluded because although their means could be converted to units of days, their standard deviations could not be converted to match units of days. In cases where multiple body sections of the same mosquitoes were measured to produce multiple size observations, only the wing measurement was included in the analysis to prevent pseudo-replication. Observations in which both the control and the predator treatments had standard deviations of zero were excluded because the meta-analysis methods did not support non-positive sampling variances.

## Exclusions and data substitutions for predator treatment means of zero

One study that was included in our database reported egg survival data as the hatch rate of field collected _Culex pervigilans_ rafts [@bib151]. However, mosquitoes have been shown to lay eggs independent of mating [@bib101], and hatch rates of zero have previously been observed in rafts laid by _Culex_ females that were held separately from males [@bib131]. Thus, hatch rates of zero were excluded from further analysis because these values may represent unfertilized egg rafts, rather than a strong impact of predators on survival. Twenty of the 187 consumptive survival observations had a predation mean of zero, and each of these zeros resulted from experiments that began with a specified number of live larvae. Consumptive survival zeros were each replaced with 0.5% of the starting number of mosquito prey to avoid undefined effect sizes. In addition, there was one zero out of the 36 oviposition predation means; this value had units of ‘number of egg rafts laid’ and was replaced with 0.5 rafts. Similar methods for replacing zero values in the treatment mean with small non-zero values have previously been employed [@bib134].

The final analysis dataset included seven subsets: consumptive effects on survival, non-consumptive effects on survival, oviposition, development (female and male), and size (female and male). The data included 187 observations from 34 studies of consumptive survival, 24 observations from seven studies of non-consumptive survival, 36 observations from 12 studies of oviposition, 14 observations from seven studies of female development, 14 observations from seven studies of male development, 27 observations from 10 studies of female size, and 18 observations from nine studies of male size ([Figure 1](#fig1)). These observations covered seven different classes of predator families ([Figure 2](#fig2)).

figure: Figure 2.
:::
![](article-71503.rmd.media/fig2.jpg)

### Mosquito predator classes (bold font) and families (italicized font) included in the database and the vector traits that they may influence (in parentheses); predator images not to scale, and placed randomly with respect to the different mosquito life stages.

Image sources: [phylopic.org](https://phylopic.org/) (CC BY 3.0 or public domain): Actinopterygii (creator: Milton Tan), Arachnida (creators: Sidney Frederic Harmer & Arthur Everett Shipley, vectorized by Maxime Dahirel), Branchiopoda (creator: Africa Gomez), and Insecta (creator: Marie Russell). [BioRender.com](https://biorender.com/): Amphibia, Hexanauplia, and Malacostraca class silhouettes; mosquito larval instars, pupa, and blood-feeding adult. Trishna Desai: mosquito egg raft.
:::
{#fig2}

## Data analysis

```{r warning=FALSE}
options(stringsAsFactors = F, scipen = 999)

library("metafor")
library("tidyverse")
library("broom")

## Full database:
database = read.csv("Predator_effects_on_mosquitoes_database.csv")

length(unique(database$Study_ID))

## Data exclusions & zero substitutions:
clean_data = database %>%
  #Insecticide used with predators
  filter(Pesticide!="YES") %>%
  #Multiple prey species
  filter(Mixed_prey!="YES") %>%
  #Traits and behavior categories that didn't have enough data for meta-analysis
  filter(Vec_trait!="OTHER" & Vec_trait!="FECUNDITY" & Vec_trait!="COMPETENCE") %>%
  filter(!(grepl("BEHAVIOR", Vec_trait) & !grepl("oviposition", Vec_trait))) %>%
  #Size or development data that were not specific to sex
  filter(!(Vec_trait=="DEVELOPMENT" & (Sex=="BOTH" | is.na(Sex)))) %>%
  filter(!(Vec_trait=="SIZE" & (Sex=="BOTH" | is.na(Sex)))) %>%
  #Pseudo-replication due to uneven representation of reps
  filter(!(grepl("individual means", Units))) %>%
  #Immature life stage-specific survival reported after more than ten days
  filter(!(Pred_effect=="CONSUMPTIVE" & Study_ID=="U106" & grepl("larval", Units))) %>%
  filter(!(Pred_effect=="CONSUMPTIVE" & Study_ID==1175 & is.na(Sex))) %>%
  filter(!(Pred_effect=="CONSUMPTIVE" & Study_ID=="U105" & grepl("larvae", Units))) %>%
  filter(!(Pred_effect=="CONSUMPTIVE" & Study_ID=="U105" & grepl("pupae", Units))) %>%
  #Inverse dev time
  filter(!(Vec_trait=="DEVELOPMENT" & grepl("1/median", Units))) %>%
  #Pseudo-replication due to measuring different body parts of the same mosquitoes
  filter(!(Study_ID=="1224" & grepl("thorax", Units))) %>%
  filter(!(Study_ID=="1226" & grepl("thorax", Units))) %>%
  filter(!(Study_ID=="U166" & grepl("femur", Other_trait))) %>%
  filter(!(Study_ID=="U166" & grepl("thorax", Other_trait))) %>%
  #Both the control sd and the predation sd are zeros. This will result in non-positive sampling variances for the ROM measure of effect.
  filter(!(Con_sd==0 & Pred_sd==0)) %>%
  #For oviposition observations with a predation treatment mean of zero, replace zero egg rafts with 0.5 egg rafts.
  mutate(Pred_mean = ifelse(Pred_mean==0 & (grepl("oviposition", Vec_trait)), 0.5, Pred_mean)) %>%
  #Delete the egg stage survival data if the predation treatment has a mean of zero. These could represent unfertilized egg rafts.
  filter(!(Vec_trait=="SURVIVAL" & Vec_stage=="EGG" & Pred_mean==0)) %>%
  #Replace larval survival observations that have zero predator means with 0.5% of the starting larvae.
  mutate(Pred_mean = ifelse(Pred_mean==0 & (grepl("%", Units) | grepl("percent", Units)), 0.5,
                            ifelse(Pred_mean==0 & (grepl("number", Units) | grepl("#", Units) | grepl("abundance", Units)), (0.005*Num_prey),
                                   ifelse(Pred_mean==0 & grepl("proportion", Units), 0.005, Pred_mean)))) %>%
  #Add the Effect ID column.
  mutate(Effect_ID = 1: n())

###############################################################################
#########################  CONSUMPTIVE SURVIVAL  ##############################
###############################################################################

consumptive_survival = clean_data %>%
  filter(Pred_effect == "CONSUMPTIVE")

length(unique(consumptive_survival$Study_ID))

#Calculate the Ratio of Means (ROM) effect size estimate.
ROM_consumptive_survival = escalc(measure = "ROM",
                                  n1i = Pred_n,
                                  n2i = Con_n,
                                  m1i = Pred_mean,
                                  m2i = Con_mean,
                                  sd1i = Pred_sd,
                                  sd2i = Con_sd,
                                  data = consumptive_survival,
                                  var.names = c("ROM_est", "ROM_sv"),
                                  append=T)

#Random effects model:
ROM.uni_consumptive_survival <- rma.uni(yi = ROM_est,
                                        vi = ROM_sv,
                                        data = ROM_consumptive_survival,
                                        method = "REML")
ROM.uni_consumptive_survival[["version"]]=NULL

coef(summary(ROM.uni_consumptive_survival))

#Assess publication bias visually with a funnel plot.
funnel(ROM.uni_consumptive_survival)

#Use Egger's regression test to quantitatively assess publication bias.
regtest(ROM.uni_consumptive_survival, model = "rma", predictor="sei")
```

```{r warnings=FALSE}
#Multilevel mixed effects models:
#No moderators
ROM.no_mods <- rma.mv(yi = ROM_est,
                      V = ROM_sv,
                      random = list(~1|Effect_ID, ~1|Study_ID/Experimental_ID),
                      data = ROM_consumptive_survival,
                      method = "REML")
ROM.no_mods[["version"]]=NULL
coef(summary(ROM.no_mods))
```

```{r warning=FALSE}
#Vector genus as a moderator
ROM.vec <- rma.mv(yi = ROM_est,
                  V = ROM_sv,
                  mods = ~Vec_genus-1,
                  random = list(~1|Effect_ID, ~1|Study_ID/Experimental_ID),
                  data = ROM_consumptive_survival,
                  method = "REML")
ROM.vec[["version"]]=NULL

coef(summary(ROM.vec))
```

```{r warning=FALSE}

#Predator family as a moderator
ROM.pred_fam <- rma.mv(yi = ROM_est,
                       V = ROM_sv,
                       mods = ~Pred_family-1,
                       random = list(~1|Effect_ID, ~1|Study_ID/Experimental_ID),
                       data = ROM_consumptive_survival,
                       method = "REML")
ROM.pred_fam[["version"]]=NULL
coef(summary(ROM.pred_fam))
```

```{r warning=FALSE}
#Predator family and vector genus as moderators
ROM.pred_and_vec <- rma.mv(yi = ROM_est,
                           V = ROM_sv,
                           mods = ~Pred_family+Vec_genus -1,
                           random = list(~1|Effect_ID, ~1|Study_ID/Experimental_ID),
                           data = ROM_consumptive_survival,
                           method = "REML")
ROM.pred_and_vec[["version"]]=NULL
coef(summary(ROM.pred_and_vec))
```

```{r warning=FALSE}
#Predator family and vector genus interacting as moderators
ROM.pred_vec_interact <- rma.mv(yi = ROM_est,
                                V = ROM_sv,
                                mods = ~Pred_family*Vec_genus -1,
                                random = list(~1|Effect_ID, ~1|Study_ID/Experimental_ID),
                                data = ROM_consumptive_survival,
                                method = "REML")
ROM.pred_vec_interact[["version"]]=NULL
coef(summary(ROM.pred_vec_interact))
```

```{r message=FALSE, warning=FALSE}
#182 out of 187 consumptive observations, 97%, are of the larval stage.
#Testing moderators on observations of larval survival only:

larval_stage = consumptive_survival %>%
  filter(!(is.na(Larval_instar)))

length(unique(larval_stage$Study_ID))

#Calculate the Ratio of Means (ROM).
ROM_larval_stage = escalc(measure = "ROM",
                          n1i = Pred_n,
                          n2i = Con_n,
                          m1i = Pred_mean,
                          m2i = Con_mean,
                          sd1i = Pred_sd,
                          sd2i = Con_sd,
                          data = larval_stage,
                          var.names = c("ROM_est", "ROM_sv"),
                          append=T)

#Random effects model:
ROM.uni_larval_stage <- rma.uni(yi = ROM_est,
                                vi = ROM_sv,
                                data = ROM_larval_stage,
                                method = "REML")
ROM.uni_larval_stage[["version"]]=NULL
coef(summary(ROM.uni_larval_stage))

#Assess publication bias visually with a funnel plot.
funnel(ROM.uni_larval_stage)

#Use Egger's regression test to quantitatively assess publication bias.
regtest(ROM.uni_larval_stage, model = "rma", predictor="sei")

#Multilevel mixed effects models:
#No moderators
ROM.larval_no_mods <- rma.mv(yi = ROM_est,
                             V = ROM_sv,
                             random = list(~1|Effect_ID, ~1|Study_ID/Experimental_ID),
                             data = ROM_larval_stage,
                             method = "REML")
coef(summary(ROM.larval_no_mods))
```

```{r message=FALSE, warning=FALSE}
#Larval instar as a moderator
ROM.larval_instar <- rma.mv(yi = ROM_est,
                            V = ROM_sv,
                            mods = ~Larval_instar-1,
                            random = list(~1|Effect_ID, ~1|Study_ID/Experimental_ID),
                            data = ROM_larval_stage,
                            method = "REML")
ROM.larval_instar[["version"]]=NULL
coef(summary(ROM.larval_instar))
```

```{r message=FALSE, warning=FALSE}
#Vector genus as a moderator
ROM.larval_vec <- rma.mv(yi = ROM_est,
                         V = ROM_sv,
                         mods = ~Vec_genus-1,
                         random = list(~1|Effect_ID, ~1|Study_ID/Experimental_ID),
                         data = ROM_larval_stage,
                         method = "REML")
ROM.larval_vec[["version"]]=NULL
coef(summary(ROM.larval_vec))
```

```{r message=FALSE, warning=FALSE}
#Predator family as a moderator
ROM.larval_pred <- rma.mv(yi = ROM_est,
                          V = ROM_sv,
                          mods = ~Pred_family-1,
                          random = list(~1|Effect_ID, ~1|Study_ID/Experimental_ID),
                          data = ROM_larval_stage,
                          method = "REML")
ROM.larval_pred[["version"]]=NULL
coef(summary(ROM.larval_pred))
```

```{r message=FALSE, warning=FALSE}
#Larval instar and predator family as moderators
ROM.larval_instar_and_pred <- rma.mv(yi = ROM_est,
                                     V = ROM_sv,
                                     mods = ~Larval_instar+Pred_family-1,
                                     random = list(~1|Effect_ID, ~1|Study_ID/Experimental_ID),
                                     data = ROM_larval_stage,
                                     method = "REML")
ROM.larval_instar_and_pred[["version"]]=NULL
coef(summary(ROM.larval_instar_and_pred))
```

```{r message=FALSE, warning=FALSE}
#Vector genus and predator family as moderators
ROM.larval_pred_and_vec <- rma.mv(yi = ROM_est,
                                  V = ROM_sv,
                                  mods = ~Vec_genus+Pred_family-1,
                                  random = list(~1|Effect_ID, ~1|Study_ID/Experimental_ID),
                                  data = ROM_larval_stage,
                                  method = "REML")
ROM.larval_pred_and_vec[["version"]]=NULL
coef(summary(ROM.larval_pred_and_vec))
```

```{r message=FALSE, warning=FALSE}
#Larval instar and predator family interacting as moderators
ROM.larval_instar_pred_interact <- rma.mv(yi = ROM_est,
                                          V = ROM_sv,
                                          mods = ~Larval_instar*Pred_family-1,
                                          random = list(~1|Effect_ID, ~1|Study_ID/Experimental_ID),
                                          data = ROM_larval_stage,
                                          method = "REML")
ROM.larval_instar_pred_interact[["version"]]=NULL
coef(summary(ROM.larval_instar_pred_interact))
```

```{r message=FALSE, warning=FALSE}
#Vector genus and predator family interacting as moderators
ROM.larval_pred_vec_interact <- rma.mv(yi = ROM_est,
                                       V = ROM_sv,
                                       mods = ~Vec_genus*Pred_family-1,
                                       random = list(~1|Effect_ID, ~1|Study_ID/Experimental_ID),
                                       data = ROM_larval_stage,
                                       method = "REML")
ROM.larval_pred_vec_interact[["version"]]=NULL
coef(summary(ROM.larval_pred_vec_interact))
```

```{r message=FALSE, warning=FALSE}
#Test moderators only on observations that included temperature.
#The temperature range is from 16 to 28.25 degrees C.

temp = consumptive_survival %>%
  filter(!(is.na(Temp))) %>%
  mutate(Temp_squared = Temp*Temp)

length(unique(temp$Study_ID))

#Calculate the Ratio of Means (ROM).
ROM_temp = escalc(measure = "ROM",
                  n1i = Pred_n,
                  n2i = Con_n,
                  m1i = Pred_mean,
                  m2i = Con_mean,
                  sd1i = Pred_sd,
                  sd2i = Con_sd,
                  data = temp,
                  var.names = c("ROM_est", "ROM_sv"),
                  append=T)

#Random effects model:
ROM.uni_temp <- rma.uni(yi = ROM_est,
                        vi = ROM_sv,
                        data = ROM_temp,
                        method = "REML")
ROM.uni_temp[["version"]]=NULL
coef(summary(ROM.uni_temp))
```

```{r message=FALSE, warning=FALSE}
#Assess publication bias visually with a funnel plot.
funnel(ROM.uni_temp)

#Use Egger's regression test to quantitatively assess publication bias.
regtest(ROM.uni_temp, model = "rma", predictor="sei")

#Multilevel mixed effects models:
#No moderators
ROM.temp_no_mods <- rma.mv(yi = ROM_est,
                           V = ROM_sv,
                           random = list(~1|Effect_ID, ~1|Study_ID/Experimental_ID),
                           data = ROM_temp,
                           method = "REML")
ROM.temp_no_mods[["version"]]=NULL
coef(summary(ROM.temp_no_mods))
```

```{r message=FALSE, warning=FALSE}
#Temp as a moderator
ROM.temp <- rma.mv(yi = ROM_est,
                   V = ROM_sv,
                   mods = ~Temp,
                   random = list(~1|Effect_ID, ~1|Study_ID/Experimental_ID),
                   data = ROM_temp,
                   method = "REML")
ROM.temp[["version"]]=NULL
coef(summary(ROM.temp))
```

```{r message=FALSE, warning=FALSE}
#Temp_squared as a moderator
ROM.temp_squared <- rma.mv(yi = ROM_est,
                           V = ROM_sv,
                           mods = ~Temp_squared,
                           random = list(~1|Effect_ID, ~1|Study_ID/Experimental_ID),
                           data = ROM_temp,
                           method = "REML")
ROM.temp_squared[["version"]]=NULL
coef(summary(ROM.temp_squared))
```

```{r message=FALSE, warning=FALSE}
###############################################################################
######################  NON-CONSUMPTIVE SURVIVAL  #############################
###############################################################################

non_con_survival = clean_data %>%
  filter(Pred_effect=="NONCONSUMPTIVE") %>%
  filter(Vec_trait=="SURVIVAL")

length(unique(non_con_survival$Study_ID))

#Calculate the Ratio of Means (ROM).
ROM_non_con_survival = escalc(measure = "ROM",
                              n1i = Pred_n,
                              n2i = Con_n,
                              m1i = Pred_mean,
                              m2i = Con_mean,
                              sd1i = Pred_sd,
                              sd2i = Con_sd,
                              data = non_con_survival,
                              var.names = c("ROM_est", "ROM_sv"),
                              append=T)

#Random effects model:
ROM.uni_non_con_survival <- rma.uni(yi = ROM_est,
                                    vi = ROM_sv,
                                    data = ROM_non_con_survival,
                                    method = "REML")
ROM.uni_non_con_survival[["version"]]=NULL
coef(summary(ROM.uni_non_con_survival))

#Assess publication bias visually with a funnel plot.
funnel(ROM.uni_non_con_survival)

#Use Egger's regression test to quantitatively assess publication bias.
regtest(ROM.uni_non_con_survival, model = "rma", predictor="sei")
```

```{r message=FALSE, warning=FALSE}
#There is significant publication bias, so the "trimfill" function is used to impute values from missing studies. 
trimfill(ROM.uni_non_con_survival, estimator = "L0", side = "left")
trimfill(ROM.uni_non_con_survival, estimator = "L0", side = "right")

#Visualization of the three imputed values on the left side of the funnel plot:
funnel(trimfill(ROM.uni_non_con_survival, estimator = "L0", side = "left"))
```

```{r message=FALSE, warning=FALSE}
###############################################################################
##############################  OVIPOSITION  ##################################
###############################################################################

oviposition = clean_data %>%
  filter(Vec_trait=="BEHAVIOR: oviposition")

length(unique(oviposition$Study_ID))

#Calculate the Ratio of Means (ROM).
ROM_oviposition = escalc(measure = "ROM",
                         n1i = Pred_n,
                         n2i = Con_n,
                         m1i = Pred_mean,
                         m2i = Con_mean,
                         sd1i = Pred_sd,
                         sd2i = Con_sd,
                         data = oviposition,
                         var.names = c("ROM_est", "ROM_sv"),
                         append=T)

#Random effects model:
ROM.uni_oviposition <- rma.uni(yi = ROM_est,
                               vi = ROM_sv,
                               data = ROM_oviposition,
                               method = "REML")
ROM.uni_oviposition[["version"]]=NULL
coef(summary(ROM.uni_oviposition))
```

```{r message=FALSE, warning=FALSE}
#Assess publication bias visually with a funnel plot.
funnel(ROM.uni_oviposition)

#Use Egger's regression test to quantitatively assess publication bias.
regtest(ROM.uni_oviposition, model = "rma", predictor="sei")

#Multilevel mixed effects models:
#No moderators
ROM.oviposition_no_mods <- rma.mv(yi = ROM_est,
                                  V = ROM_sv,
                                  random = list(~1|Effect_ID, ~1|Study_ID/Experimental_ID),
                                  data = ROM_oviposition,
                                  method = "REML")
ROM.oviposition_no_mods[["version"]]=NULL
coef(summary(ROM.oviposition_no_mods))
```

```{r message=FALSE, warning=FALSE}
#Vector genus as a moderator
ROM.oviposition_vec <- rma.mv(yi = ROM_est,
                              V = ROM_sv,
                              mods = ~Vec_genus-1,
                              random = list(~1|Effect_ID, ~1|Study_ID/Experimental_ID),
                              data = ROM_oviposition,
                              method = "REML")
ROM.oviposition_vec[["version"]]=NULL
coef(summary(ROM.oviposition_vec))
```

```{r message=FALSE, warning=FALSE}
#Predator family as a moderator
ROM.oviposition_pred <- rma.mv(yi = ROM_est,
                               V = ROM_sv,
                               mods = ~Pred_family-1,
                               random = list(~1|Effect_ID, ~1|Study_ID/Experimental_ID),
                               data = ROM_oviposition,
                               method = "REML")
ROM.oviposition_pred[["version"]]=NULL
coef(summary(ROM.oviposition_pred))
```

```{r message=FALSE, warning=FALSE}
###############################################################################
##########################  FEMALE DEVELOPMENT  ###############################
###############################################################################

development_females = clean_data %>%
  filter(Vec_trait=="DEVELOPMENT") %>%
  filter(Sex=="FEMALE") 

length(unique(development_females$Study_ID))

#Calculate the Ratio of Means (ROM).
ROM_development_females = escalc(measure = "ROM",
                                 n1i = Pred_n,
                                 n2i = Con_n,
                                 m1i = Pred_mean,
                                 m2i = Con_mean,
                                 sd1i = Pred_sd,
                                 sd2i = Con_sd,
                                 data = development_females,
                                 var.names = c("ROM_est", "ROM_sv"),
                                 append=T)

#Random effects model:
ROM.uni_development_females <- rma.uni(yi = ROM_est,
                                       vi = ROM_sv,
                                       data = ROM_development_females,
                                       method = "REML")

ROM.uni_development_females[["version"]]=NULL
coef(summary(ROM.uni_development_females))

#Assess publication bias visually with a funnel plot.
funnel(ROM.uni_development_females)

#Use Egger's regression test to quantitatively assess publication bias.
regtest(ROM.uni_development_females, model = "rma", predictor="sei")

#There is significant publication bias, so the "trim-fill" function is used to impute values from missing studies. 
trimfill(ROM.uni_development_females, estimator = "L0", side = "left")
trimfill(ROM.uni_development_females, estimator = "L0", side = "right")
```


```{r message=FALSE, warning=FALSE}
###############################################################################
###########################  MALE DEVELOPMENT  ################################
###############################################################################

development_males = clean_data %>%
  filter(Vec_trait=="DEVELOPMENT") %>%
  filter(Sex=="MALE") 

length(unique(development_males$Study_ID))

#Calculate the Ratio of Means (ROM).
ROM_development_males = escalc(measure = "ROM",
                               n1i = Pred_n,
                               n2i = Con_n,
                               m1i = Pred_mean,
                               m2i = Con_mean,
                               sd1i = Pred_sd,
                               sd2i = Con_sd,
                               data = development_males,
                               var.names = c("ROM_est", "ROM_sv"),
                               append=T)

#Random effects model:
ROM.uni_development_males <- rma.uni(yi = ROM_est,
                                     vi = ROM_sv,
                                     data = ROM_development_males,
                                     method = "REML")
ROM.uni_development_males[["version"]]=NULL
coef(summary(ROM.uni_development_males))

#Assess publication bias visually with a funnel plot.
funnel(ROM.uni_development_males)

#Use Egger's regression test to quantitatively assess publication bias.
regtest(ROM.uni_development_males, model = "rma", predictor="sei")

#There is significant publication bias, so the "trimfill" function is used to impute values from missing studies. 
trimfill(ROM.uni_development_males, estimator = "L0", side = "left")
trimfill(ROM.uni_development_males, estimator = "L0", side = "right")
```

```{r message=FALSE, warning=FALSE}
###############################################################################
##############################  FEMALE SIZE  ##################################
###############################################################################

size_females = clean_data %>%
  filter(Vec_trait=="SIZE") %>%
  filter(Sex=="FEMALE") 

length(unique(size_females$Study_ID))

#Calculate the Ratio of Means (ROM).
ROM_size_females = escalc(measure = "ROM",
                          n1i = Pred_n,
                          n2i = Con_n,
                          m1i = Pred_mean,
                          m2i = Con_mean,
                          sd1i = Pred_sd,
                          sd2i = Con_sd,
                          data = size_females,
                          var.names = c("ROM_est", "ROM_sv"),
                          append=T)

#Random effects model:
ROM.uni_size_females <- rma.uni(yi = ROM_est,
                                vi = ROM_sv,
                                data = ROM_size_females,
                                method = "REML")

ROM.uni_development_males[["version"]]=NULL
coef(summary(ROM.uni_size_females))

#Assess publication bias visually with a funnel plot.
funnel(ROM.uni_size_females)

#Use Egger's regression test to quantitatively assess publication bias.
regtest(ROM.uni_size_females, model = "rma", predictor="sei")

#There is significant publication bias, so the "trimfill" function is used to impute values from missing studies. 
trimfill(ROM.uni_size_females, estimator = "L0", side = "left")
trimfill(ROM.uni_size_females, estimator = "L0", side = "right")
```

```{r message=FALSE, warning=FALSE}
#Visualization of the two imputed values on the right side of the funnel plot:
funnel(trimfill(ROM.uni_size_females, estimator = "L0", side = "right"))
```

```{r warning=FALSE}
###############################################################################
###############################  MALE SIZE  ###################################
###############################################################################

size_males = clean_data %>%
  filter(Vec_trait=="SIZE") %>%
  filter(Sex=="MALE") 

length(unique(size_males$Study_ID))

#Calculate the Ratio of Means (ROM).
ROM_size_males = escalc(measure = "ROM",
                        n1i = Pred_n,
                        n2i = Con_n,
                        m1i = Pred_mean,
                        m2i = Con_mean,
                        sd1i = Pred_sd,
                        sd2i = Con_sd,
                        data = size_males,
                        var.names = c("ROM_est", "ROM_sv"),
                        append=T)

#Random effects model:
ROM.uni_size_males <- rma.uni(yi = ROM_est,
                              vi = ROM_sv,
                              data = ROM_size_males,
                              method = "REML")
ROM.uni_size_males[["version"]]=NULL
coef(summary(ROM.uni_size_males))

#Assess publication bias visually with a funnel plot.
funnel(ROM.uni_size_males)

#Use Egger's regression test to quantitatively assess publication bias.
regtest(ROM.uni_size_males, model = "rma", predictor="sei")
```

### Measuring effect sizes and heterogeneity

All analyses were conducted in R version 4.0.2 [@bib108]. For each subset of trait data ([Figure 1](#fig1)), the ratio of means (ROM) measure of effect size was calculated using the ‘escalc’ function from the ‘metafor’ package; this effect measure is equal to a log-transformed fraction, where predation mean is the numerator and control mean is the denominator [@bib143]. Random effects models, using the ‘rma.uni’ function, were run with the ROM effect sizes as response variables; each model had a normal error distribution and a restricted maximum likelihood (REML) estimator for τ^2^, the variance of the distribution of true effect sizes [@bib143]. Although these random effects models could not account for multiple random effects or moderators, they provided overall estimates of the ROM effect sizes and estimates of the _I_^2^ statistics. Each _I_^2^ statistic represented the percentage of total variation across studies due to heterogeneity [@bib60]. If the _I_^2^ statistic was equal to or greater than 75%, the heterogeneity was considered to be high [@bib60], and high heterogeneity has previously motivated further testing of moderators [@bib144].

### Assessing publication bias

Publication bias was assessed by visually inspecting funnel plots and conducting Egger’s regression test (‘regtest’ function) with standard error as the predictor [@bib129; @bib143]. If the Egger’s regression test showed significant evidence of publication bias based on funnel plot asymmetry, the ‘trim and fill’ method (‘trimfill’ function) was used to estimate how the predation effect size might change after imputing values from missing studies [@bib37; @bib38; @bib143]. The trim and fill method has previously been recommended for testing the robustness of conclusions related to topics in ecology and evolution [@bib64]. Of the two trim and fill estimators, R~null~ and L~null~, that were originally recommended [@bib37; @bib38], the L~null~ estimator was used in this study because it is more appropriate for smaller datasets [@bib123].

### Testing moderators

Data subsets that had high heterogeneity, observations from at least 10 studies, and no evidence of publication bias according to Egger’s regression results were analyzed further using multilevel mixed effects models with the ‘rma.mv’ function [@bib143; @bib61]. All multilevel mixed effects models had normal error distributions, REML estimators for τ^2^, and accounted for two random factors: effect size ID, and experiment ID nested within study ID. Moderators, such as predator family, vector genus, larval instar (directly correlated to prey size), and temperature, were tested within each data subset to determine if they affected the observed heterogeneity in ROM effect sizes. For categorical moderators, the intercept of the multilevel mixed effects model was removed, allowing an analysis of variance (ANOVA) referred to as the ‘test of moderators’ to indicate if any of the categories had an effect size different than zero. For data subsets with observations from 10 to 29 studies, only one moderator was tested at a time to account for sample size constraints. For subsets with observations from a higher number of studies (30 or more), up to two moderators were tested at once, and interaction between moderators was also tested. The small sample corrected Akaike Information Criterion (AICc) was used to compare multilevel mixed effects models and to select the model of best fit within each data subset; differences in AICc greater than two were considered meaningful [@bib22].

# Results

## Random effects models

Each data subset ([Figure 1](#fig1)) had an _I_^2^ statistic of greater than 75%, indicating high heterogeneity [@bib60]. Random effects model results showed that predators consumptively decreased mosquito survival with an effect size of –1.23 (95% CI −1.43,–1.03), p-value &lt; 0.0001, and non-consumptively reduced survival with a smaller effect size of –0.11 (95% CI −0.17,–0.04), p-value = 0.0016. In addition, predators non-consumptively reduced oviposition behavior with an effect size of –0.87 (95% CI −1.31,–0.42), p-value = 0.0001, and mosquito body size was non-consumptively reduced by predators in both males and females; the female effect size was –0.13 (95% CI −0.19,–0.06), p-value = 0.0002, and the male effect size was –0.03 (95% CI −0.06,–0.01), p-value = 0.0184. There was not a significant non-consumptive effect of predators on either male or female development time; the female effect size was –0.01 (95% CI –0.09, 0.07), p-value = 0.7901, and the male effect size was –0.04 (95% CI –0.12, 0.04), p-value = 0.3273.

The Egger’s regression test results showed that the non-consumptive survival subset, both development time subsets (male and female), and the female size subset exhibited funnel plot asymmetry indicative of publication bias. The ‘trim and fill’ procedure identified missing studies in the non-consumptive survival subset and the female size subset, but the procedure did not identify any missing studies in either of the development time subsets. Three studies were estimated to be missing from the non-consumptive survival data, and accounting for imputed values from missing studies resulted in a shift in the predation effect size from –0.11 (95% CI −0.17,–0.04), p-value = 0.0016, to -0.13 (95% CI −0.20,–0.07), p-value &lt; 0.0001. Two studies were estimated to be missing from the female size data, and accounting for imputed values from these missing studies shifted the predation effect size from –0.13 (95% CI −0.19,–0.06), p-value = 0.0002, to -0.10 (95% CI −0.17,–0.03), p-value = 0.0083. Shifts in effect size estimates due to the trim and fill procedure were minor and did not cause any of the observed effects of predators to change direction or become insignificant.

## Multilevel mixed effects models

The consumptive survival and oviposition data subsets met the criteria of high heterogeneity, observations from at least 10 studies, and no evidence of publication bias. Therefore, these data subsets were tested for moderators using multilevel mixed effects models. Predator families that decreased mosquito survival included Cyprinidae: –3.44 (95% CI −5.79,–1.09), p-value = 0.0042; Poeciliidae: –1.42 (95% CI −2.67,–0.16), p-value = 0.0270; Ambystomatidae: –5.18 (95% CI −7.94,–2.42), p-value = 0.0002; Aeshnidae: –2.93 (95% CI −4.80,–1.07), p-value = 0.0020; and Notonectidae: –2.14 (95% CI −3.07,–1.21), p-value &lt; 0.0001 ([Figure 3a](#fig3)). Vector genera that experienced significant decreases in survival due to consumptive effects of predators included _Aedes_: –1.23 (95% CI −1.81,–0.65), p-value &lt; 0.0001; _Anopheles_: –1.34 (95% CI −2.01,–0.66), p-value = 0.0001; and _Culex_: –1.41 (95% CI −1.96,–0.86), p-value &lt; 0.0001 ([Figure 3b](#fig3)). Among all 187 consumptive survival observations from 34 studies, the best model fit, according to AICc value, was achieved when an interaction between predator family and vector genus was included in the model ([Table 2](#table2)). However, among the 163 larval stage consumptive survival observations from 30 studies, adding an interactive term between larval instar (an indicator of prey size) and predator family had a greater improvement on model fit than adding an interactive term between vector genus and predator family ([Figure 3c](#fig3), [Table 3](#table3)). Temperature did not affect the heterogeneity of consumptive survival data, either as a linear moderator: –0.01 (95% CI –0.10, 0.07), p-value = 0.7559, or a quadratic moderator: 0.00 (95% CI 0.00, 0.00), p-value = 0.8184. The best oviposition model fit, according to AICc value, was achieved when vector genus was added as a moderator ([Table 4](#table4)). The mean oviposition effect size was not significantly different than zero for _Aedes_: 0.32 (95% CI –2.14, 2.79), p-value = 0.7970, or _Culiseta_: –0.61 (95% CI –1.83, 0.62), p-value = 0.3329, but for _Culex_ mosquitoes, oviposition was significantly decreased by predator presence: –1.69 (95% CI −2.82,–0.56), p-value = 0.0033 ([Figure 4](#fig4)).

table: Table 2.
:::
### Candidate multilevel mixed effects models of consumptive effects from predators on mosquito survival, fitted to dataset of effect sizes (n = 187 from 34 studies), and ranked by corrected Akaike’s information criterion (AICc).

| Moderator(s)                   | Test of moderators(degrees of freedom, p-value) | AICc  | ΔAICc |
| ------------------------------ | ----------------------------------------------- | ----- | ----- |
| Predator family x vector genus | 28, &lt; 0.0001                                 | 500.5 | !null |
| Predator family                | 19, &lt; 0.0001                                 | 507   | 6.5   |
| Predator family + vector genus | 23, &lt; 0.0001                                 | 508.1 | 7.6   |
| Vector genus                   | 5, &lt; 0.0001                                  | 573   | 72.5  |
| None                           | ----                                            | 576.5 | 76    |
:::
{#table2}

table: Table 3.
:::
### Candidate multilevel mixed effects models of consumptive effects from predators, fitted to dataset of effect sizes where larval instar is not missing (n = 163 from 30 studies), and ranked by corrected Akaike’s information criterion (AICc).

| Moderator(s)                    | Test of moderators(degrees of freedom, p-value) | AICc  | ΔAICc |
| ------------------------------- | ----------------------------------------------- | ----- | ----- |
| Predator family x larval instar | 25, &lt; 0.0001                                 | 429.2 | !null |
| Predator family + larval instar | 19, &lt; 0.0001                                 | 443.5 | 14.3  |
| Predator family x vector genus  | 25, &lt; 0.0001                                 | 455   | 25.8  |
| Predator family                 | 17, &lt; 0.0001                                 | 456.8 | 27.6  |
| Predator family + vector genus  | 21, &lt; 0.0001                                 | 458.4 | 29.2  |
| Larval instar                   | 3, &lt; 0.0001                                  | 503.1 | 73.9  |
| Vector genus                    | 5, &lt; 0.0001                                  | 504.7 | 75.5  |
| None                            | ----                                            | 508.5 | 79.3  |
:::
{#table3}

table: Table 4.
:::
### Candidate multilevel mixed effects models of non-consumptive effects of predators on mosquito oviposition behavior, fitted to dataset of effect sizes (n = 36 from 12 studies), and ranked by corrected Akaike’s information criterion (AICc).

| Moderator(s)    | Test of moderators(degrees of freedom, p-value) | AICc  | ΔAICc |
| --------------- | ----------------------------------------------- | ----- | ----- |
| Vector genus    | 3, 0.0149                                       | 122.1 | !null |
| None            | ----                                            | 125.2 | 3.1   |
| Predator family | 12, 0.8855                                      | 167.9 | 45.8  |
:::
{#table4}

figure: Figure 3.
:::
![](article-71503.rmd.media/fig3.jpg)

### Effect sizes and 95 % confidence intervals for consumptive effects of predators, for different categories of moderators (with number of studies in parentheses).

(**a**) predator family with predator class in the right-hand column, (**b**) vector genus, and (**c**) larval instar.
:::
{#fig3}

figure: Figure 4.
:::
![](article-71503.rmd.media/fig4.jpg)

### Oviposition effect sizes and 95 % confidence intervals for different categories of vector genus (with number of studies in parentheses).
:::
{#fig4}

# Discussion

In this study, laboratory and semi-field empirical data were obtained through a systematic literature review and used to conduct a meta-analysis that assessed consumptive and non-consumptive effects of predators on mosquito prey. Some results agree with previously observed trends, such as greater consumptive effects from larger predators [@bib74; @bib104] and no oviposition response to predator cues among container-breeding _Aedes_ mosquitoes [@bib146]. However, this meta-analysis revealed additional trends. Mosquito larval instar had an important role in moderating consumptive effects of predators, likely because of its direct correlation to prey size. Furthermore, a small, but significant, decrease in mosquito survival due to non-consumptive effects of predators was observed, suggesting that mosquitoes can be ‘scared to death’ by predators [@bib106]. Both male and female body sizes were also reduced among mosquitoes that had been exposed to predators, and predator avoidance during oviposition was observed among female _Culex_ mosquitoes. Effects of predators on different vector traits, particularly survival, body size, and oviposition behavior, have the potential to influence infectious disease dynamics.

## Consumptive effects of predators on survival

Several larger predators reduced mosquito survival, including freshwater fish (Cyprinidae and Poeciliidae), salamander larvae (Ambystomatidae), dragonfly larvae (Aeshnidae), and backswimmers (Notonectidae) ([Figure 3a](#fig3)). This finding is consistent with a previous analysis which showed a positive linear relationship between predator body mass and ingestion rate across taxa [@bib104]. In addition, more effect size heterogeneity in the consumptive survival data was explained by an interaction between predator family and larval instar than was explained by an interaction between predator family and vector genus ([Table 3](#table3)). This result suggests that the relative sizes of predator and prey groups could play a more important role in determining consumptive mosquito survival than variations in predator responses to different behaviors of prey genera, which are likely to be shaped by the degree of shared evolutionary history between trophic levels [@bib21]. Larval instar is an indicator of mosquito size, and previous modeling work has provided evidence of prey size selection by predators to maximize energetic gain [@bib92]. While smaller cyclopoid copepods are more effective against early instar mosquito larvae [@bib34], larger predators including tadpoles, giant water bugs, dragonfly larvae, fish, and backswimmers are more effective against late instar larvae [@bib75].

## Non-consumptive effects of predators on survival

Exposure to predation cues significantly lowered mosquito survival, and this non-consumptive effect has also been observed in dragonfly larvae prey (_Leucorrhinia intacta_) that were exposed to caged predators [@bib90]. The reduction in mosquito survival from non-consumptive effects of predators was significantly smaller than the reduction that was observed from consumptive effects. This is partially due to the practical constraints of most experimental designs, which cause consumptive and non-consumptive effects of predators on survival to be grouped together and reported as consumptive effects. The greater impact of combined consumptive and non-consumptive effects, in comparison to only non-consumptive effects, has previously been observed in pea aphids (_Acyrthosiphon pisum_) [@bib96].

## Non-consumptive effects of predators on body size

While predators did not significantly impact mosquito development time through non-consumptive effects in either sex, mosquito body size was decreased by the non-consumptive effects of predators in both sexes. Smaller body size is associated with lower reproductive success in mosquitoes because smaller females lay fewer eggs [@bib16; @bib84; @bib98; @bib130; @bib138], and smaller males produce less sperm [@bib55; @bib105]. These effects suggest that predation could non-consumptively reduce mosquito population growth. The smaller size of mosquitoes exposed to predators could also limit disease transmission. Vector lifespan contributes disproportionately to disease transmission because older vectors are more likely to have been exposed to pathogens, more likely to already be infectious after having survived the extrinsic incubation period, and more likely to survive long enough to bite subsequent hosts [@bib23]. It is well-established that smaller mosquito body size is associated with shorter mosquito lifespan [@bib9; @bib56; @bib111; @bib112; @bib149]. Therefore, non-consumptive effects of predators may limit the transmission of mosquito-borne diseases.

## Non-consumptive effects of predators on oviposition behavior

Predator presence also non-consumptively reduced oviposition behavior in adult female mosquitoes. Meta-regression results showed that _Culex_ females significantly avoid oviposition sites that contain predators or predator cues, but _Aedes_ and _Culiseta_ females do not avoid these sites, despite a slight non-significant trend toward predator avoidance in _Culiseta_ ([Figure 4](#fig4)). Both _Culex_ and _Culiseta_ mosquitoes have an ‘all-or-none’ oviposition strategy [@bib66], in which they lay hundreds of rapidly hatching eggs in rafts on the water’s surface [@bib30]. Such an oviposition strategy is conducive to evolving predator avoidance behaviors, and a previous meta-analysis showed significant predator avoidance in both _Culex_ and _Culiseta_ during oviposition [@bib146]. Conversely, it is likely that an oviposition response to predation is not particularly advantageous for _Aedes_ because the delayed hatching of their eggs [@bib30] can prevent the level of predation risk at the time of oviposition from matching the level of predation risk present in the eventual larval environment [@bib146]. The predator avoidance response in _Aedes_ species that lay their eggs above the water’s edge in containers has previously been described as ‘non-existent’ [@bib146]. Both _Aedes_ species included in this study’s oviposition data subset, _Ae. albopictus_ and _Ae. aegypti_, meet the criterion of ovipositing above water in containers [@bib67]. Predator avoidance during oviposition has previously been found to increase the mosquito population size at equilibrium [@bib127]. However, this study’s results and those of a previous meta-analysis [@bib146] suggest that models of oviposition site selection, such as those using parameters from Notonectidae predators and _Culiseta_ prey [@bib70], are not generalizable to _Aedes_ vectors.

## Implications for West Nile Virus disease dynamics

Predator avoidance during oviposition by _Culex_ mosquitoes ([Figure 4](#fig4)) may be of particular importance to West Nile virus (WNV) disease dynamics. Previous work has shown that _Cx. pipiens_, _Cx. restuans_, and _Cx. tarsalis_ all avoid predator habitats [@bib146], and that _Cx. pipiens_ is the primary bridge vector of WNV responsible for spill-over transmission from avian reservoir hosts to humans [@bib44; @bib51; @bib72; @bib8]. _Cx. pipiens_ mosquitoes can live in permanent aquatic environments, such as ground pools [@bib3; @bib11; @bib33; @bib133], ponds [@bib83], stream edges [@bib3], and lake edges [@bib145] that are more common in rural areas, but _Cx. pipiens_ are also found in urban and suburban residential areas, where they typically breed in artificial containers [@bib133], including tires [@bib83; @bib97; @bib140], rainwater tanks [@bib136], and catch basins [@bib45]. Small artificial containers, such as discarded tires, are generally unlikely to harbor larger predators, including freshwater fish (Cyprinidae and Poeciliidae), salamander larvae (Ambystomatidae), dragonfly larvae (Aeshnidae), and backswimmers (Notonectidae), because temporary aquatic environments cannot support the relatively long development times of these organisms. The mean dispersal distance of adult _Culex_ mosquitoes is greater than one kilometer [@bib27; @bib53], and female _Cx. pipiens_ have exhibited longer dispersal distances after developing in the presence of a fish predator [@bib2]. Therefore, predator avoidance during oviposition may cause _Cx. pipiens_ populations to disperse from permanent aquatic environments in more rural areas to artificial container environments in urbanized areas, where the risk of human WNV infection is higher [@bib19].

Predator cue levels may be altered by climate conditions, and these changes in cue levels can impact WNV transmission to humans. Drought has previously been associated with human WNV cases [@bib65; @bib86; @bib115; @bib122; @bib42; @bib102], but the association has thus far lacked a clear underlying mechanism. Under drought conditions, the density of aquatic organisms increases and predation pressures can intensify due to compressed space and high encounter rates [@bib4]. A previous study of a stream ecosystem found that impacts of fish predation are more severe during the dry season [@bib36]. In addition, reductions in water volume can facilitate consumption of mosquito larvae by crane fly larvae (Tipulidae), whereas mosquito consumption by tipulids was not observed at a higher water level [@bib4]. Laboratory and semi-field studies have shown that mosquitoes respond to a gradient of predator cues [@bib118; @bib125]. The frequency of larval anti-predator behavior is correlated with the concentration of predator cues [@bib118], and adult female mosquitoes prefer oviposition sites with lower predator densities [@bib125]. Therefore, as predator cue levels increase due to drought, permanent aquatic habitats are likely to transition from suitable oviposition sites for one generation of female mosquitoes, to unsuitable oviposition sites for the next generation.

When suitable oviposition sites are absent, females retain their eggs until sites become available [@bib15]. _Cx. pipiens_ females can retain their eggs for up to five weeks, allowing them enough time to find container sites with low predation risk, often located in residential areas [@bib66]. The movement of gravid female _Cx. pipiens_ to residential areas increases the risk of WNV spill-over to humans because these vectors are likely to have already blood-fed at least once [@bib28], suggesting that they have a higher risk of WNV infection, relative to non-gravid mosquitoes. This is consistent with studies that have reported associations between drought and WNV-infected mosquitoes in urban and residential areas [@bib65; @bib102]. In addition, vertical transmission of WNV from gravid females to their progeny may occur during oviposition [@bib117], when the virus is transmitted by an accessory gland fluid that attaches eggs to one another [@bib95]. Because the rate of vertical transmission in _Cx. pipiens_ increases with the number of days following WNV infection [@bib7], extended searches for oviposition sites due to drought could increase the frequency of vertical transmission. However, the impact of vertical transmission on WNV epidemics is thought to be minimal because when transmission to an egg raft did occur, only 4.7% of the progeny were found to be infected as adults [@bib7], and only about half of those infected adults are estimated to be female. In summary, the movement of _Cx. pipiens_ females toward more residential areas, combined with potential limited WNV amplification from increased vertical transmission, suggests that the vector trait of predator avoidance during oviposition can serve as a plausible explanation for associations between drought and human WNV cases.

Another theory for the association between drought and human WNV cases is based on the hypothesis that increased contact between mosquito vectors and passerine reservoir hosts occurs during drought conditions [@bib102; @bib122]. The proposed aggregation of bird and mosquito populations during drought was originally thought to occur in humid, densely vegetated hammocks – a type of habitat that is specific to southern Florida [@bib122], but WNV incidence is more consistently clustered in other regions of the US, particularly the Northern Great Plains [@bib24; @bib132]. Northern cardinals (_Cardinalis cardinalis_), American robins (_Turdus migratorius_), and house sparrows (_Passer domesticus_) were among the bird species that most frequently tested seropositive for WNV antibodies in 2005 and 2006 in Chicago, where high numbers of human cases were reported [@bib52], and these passerine species are more abundant in residential areas, regardless of precipitation patterns [@bib6; @bib13; @bib79]. Apart from drought, landowners’ participation in supplemental bird feeding, providing bird houses, gardening, and maintaining vegetation can strongly influence passerine abundance in residential areas [@bib78]. Furthermore, as terrestrial foragers that can obtain hydration from their diet of insects, fruits, and other plant material [@bib5; @bib20; @bib85; @bib113], passerine reservoir hosts of WNV are less likely to move in response to drought than the mosquito vectors of WNV, which have obligate aquatic life stages.

While hatch-year birds are more vulnerable to mosquito biting, and thus contribute to the amplification of WNV [@bib52], it is illogical to expect an increased abundance of hatch-year birds during drought conditions. However, some have argued that in cases where drought decreases the abundance of juvenile birds, the ratio of mosquitoes to birds increases, and this could lead to higher WNV prevalence in the mosquito population [@bib102]. Although reductions in both hatching success [@bib46] and survival of recently fledged birds [@bib150] have been observed during drought conditions, the impact of drought on avian abundance varies widely by species [@bib141]. In particular, synanthropic species, such as those likely to harbor WNV, are less negatively affected by drought [@bib1]. Additionally, the droughts that impact avian abundance often occur over much longer periods of time than the seasonal droughts that predict WNV transmission to humans. For example, avian abundance has been modeled based on precipitation metrics spanning 32 weeks, and house wren (_Troglodytes aedon_) abundance has been predicted by precipitation averages spanning four years [@bib141]. Finally, birds with higher levels of stress hormones are more likely to be fed on by mosquitoes, and certain factors associated with residential areas, such as road noise, light pollution, and pesticide exposure, can cause avian stress [@bib47]. Therefore, elevated avian stress hormones in these habitats may contribute to WNV prevalence in the mosquito population, independent of drought conditions.

## Implications for mosquito-borne disease modeling

Although the aquatic phase of the mosquito life cycle is often overlooked in mathematical models of mosquito-borne pathogen transmission [@bib110], vector survival at immature stages plays an important role in determining mosquito population abundance, which is an essential factor for predicting disease transmission [@bib12]. The results of this study show that mosquito survival decreases among the _Aedes_, _Anopheles_, and _Culex_ genera due to consumptive effects of predators ([Figure 3b](#fig3)), and that there is also a reduction in mosquito survival due to non-consumptive effects. Other studies have demonstrated that aquatic predators dramatically impact mosquito survival and abundance. For example, a biocontrol intervention relying on the application of copepod predators eliminated _Aedes albopictus_ from three communes in Nam Dinh, Vietnam, where dengue transmission was previously detected, and reduced vector abundance by 86–98% in three other communes [@bib69]. Conversely, the annual abundance of _Culex_ and _Anopheles_ mosquitoes was observed to increase 15-fold in semi-permanent wetlands in the year following a drought, likely because the drought eliminated aquatic predators from wetlands that dried completely, and mosquitoes were able to re-colonize newly formed aquatic habitats more quickly than their most effective predators [@bib26].

While relationships between temperature and different vector traits, such as fecundity and lifespan, have been incorporated into models of temperature effects on mosquito population density [@bib40], models of predator effects on vector borne disease transmission have focused primarily on the impacts of predation on vector survival. Previous models have shown that predators of vector species can decrease or eliminate pathogen infection in host populations as vector fecundity increases [@bib94]. The findings of this meta-analysis suggest that predators also decrease vector fecundity through non-consumptive effects on vector body size. In addition, the entomological inoculation rate (EIR) is likely to be reduced by effects of predators on mosquito fecundity and lifespan, as well as effects of predators on mosquito survival. The EIR has been defined as the product of three variables: (_m_) the number of mosquitoes per host, (_a_) the daily rate of mosquito biting, and (_s_) the proportion of mosquitoes that are infectious [@bib12]. Based on this study’s findings, predators are likely to decrease the number of mosquitoes per host by reducing mosquito survival through both consumptive and non-consumptive effects, and by reducing mosquito fecundity through non-consumptive effects on body size. In addition, predators are likely to decrease the proportion of mosquitoes that are infectious by shortening the vector lifespan through non-consumptive effects on body size. The relationship between mosquito body size and biting rate is unclear, with some studies showing higher biting rates among larger mosquitoes [@bib9; @bib50], and others reporting higher biting rates among smaller mosquitoes [@bib43; @bib77]. The links between factors that influence the EIR and observed effects of predators on mosquito prey demonstrate the necessity of including both consumptive and non-consumptive effects of predators in models of mosquito-borne disease.

## Conclusion

This meta-analysis on mosquito predation demonstrates that predators not only play an important role in directly reducing mosquito populations, but also have non-consumptive effects on surviving mosquitoes that may ultimately reduce further population growth and decrease disease transmission. While families of larger sized predators were effective in reducing mosquito survival, other factors, such as impacts on native species, as well as the economic cost of mass-rearing and field applications [@bib73; @bib107], should be carefully considered before selecting a predator as a suitable biocontrol agent. Predictive disease models are likely to be more reliable when the non-consumptive effects of predation are incorporated. Although exposure of mosquito larvae to predators is commonplace in outdoor field settings, it remains rare in most laboratory-based assessments of vector traits. Therefore, mosquitoes observed in nature are likely to have smaller body sizes than those observed under optimal laboratory conditions. It is important for disease modelers to recognize these impacts of predation on vector traits as they can reduce mosquito population growth and limit disease transmission due to shorter vector lifespans. Within the WNV disease system, consideration of the oviposition behavioral response to predation cues by _Culex_ vectors can improve current understanding of the association between drought and human cases. This study provides general estimates of the effects of predators on selected mosquito traits for use in predictive disease models.

## Future directions

Modeling efforts that aim to optimize the application of biocontrol predators should also consider incorporating predator effects on vector survival, fecundity, and lifespan. These additions to predictive models of various biocontrol interventions are likely to help public health officials choose the most cost-effective strategies for limiting disease transmission. In the 60-study database that was compiled, only one study was designed to directly measure the effect of larval-stage predation on vector competence [@bib119]. Therefore, future efforts to assess the impact of predators on mosquito-borne disease transmission should prioritize experimental studies in which infected mosquito larvae are observed throughout an initial period of aquatic exposure to predators, followed by a period of blood-feeding in the adult stage.

Two studies from the compiled database examined the compatibility of predators with _Bacillus thuringiensis_ var. _israelensis_ (_Bti_), a commonly used bacterial biocontrol agent [@bib25; @bib99]. Previous studies have supported the simultaneous application of cyclopoid copepod predators and _Bti_!number(0)[@bib88; @bib135], but additional analyses are needed on the use of _Bti_ with other families of mosquito predators. Populations of other insect pests, such as the southern green stink bug (_Nezara viridula_), are known to be regulated by both predators and parasites [@bib39]. The literature search conducted for this meta-analysis returned studies on water mite parasites [@bib109] and nematode parasitoids [@bib31] of mosquitoes, and ascogregarine parasites have previously been evaluated as biocontrol agents against _Aedes_ mosquitoes [@bib137]. A more thorough review of the impacts of parasites and parasitoids on vector traits, such as survival, fecundity, and lifespan, is needed before incorporating these potential biocontrol agents into integrated vector control plans.

Three studies in the 60-study database included experiments where two mosquito prey species were made available to the predator species [@bib48; @bib49; @bib91]. In these cases, the effect size measurement for each mosquito species could be influenced by interspecific competition, or a preference of the predator species for a certain prey species. Hetero-specific prey observations were excluded from this meta-analysis, but future analyses centered on the concepts of interspecific competition or predator preferences might further evaluate these data. In addition, this meta-analysis investigated consumptive and non-consumptive effects of predators separately. More research is needed to determine how models should combine these different types of predator effects to accurately reflect predation interactions as they occur in natural environments.
