|Home | About | Journals | Submit | Contact Us | Français|
Evolution of osmoregulatory systems is a key factor in the transition of species between fresh- and saltwater habitats. Anopheles coluzzii and An. merus are stenohaline and euryhaline malaria vector mosquitoes belonging to a larger group of sibling species, the Anopheles gambiae complex, which radiated in Africa within the last 2 million years. Comparative ecological genomics of these vector species can provide insight into the mechanisms that permitted the rapid radiation of this species complex into habitats of contrasting salinity. Here we use RNA-Seq to investigate gene expression differences between An. coluzzii and An. merus after briefly exposing both young and old larval instars of each species to either saltwater (SW) or freshwater (FW). Our examination aims to identify candidate genes and pathways responsible for the greater SW tolerance of An. merus. Our results are congruent with the ability of gene induction to mediate salinity tolerance, with both species showing increasing amounts of differential gene expression between SW and FW as salt concentrations increase. Besides ion transporters such as AgAE2 that may serve as effectors for osmoregulation, we also find mitogen activated protein kinases that may serve in a phosphorylation signaling pathway responding to salinity, and report potential cross-talk between the mosquito immune response and osmoregulation. This study provides a key step towards applying the growing molecular knowledge of these malaria vectors to improve understanding of their ecological tolerances and habitat occupancy.
Transitions between freshwater (FW) and saltwater (SW) systems can represent dramatic shifts, requiring adaptation to sudden fluxes of water and ions. Salinity has been described as the most important physical environmental influence on species’ distributions (Whitehead et al. 2011). Within fishes, salinity tolerance differs both among closely related species and within populations of the same species, with killifish and sticklebacks representing two well-studied examples (Jones et al. 2012; Kozak et al. 2013). Differences in saltwater tolerance have been implicated in the stable coexistence of sibling species in cryptic species complexes and local adaptation of aquatic invertebrates such as rotifers (Montero-Pau & Serra 2011; Scheuerl & Stelzer 2013). Rapid evolution of enzymes (V-type H+ ATPase and Na+/K+ ATPase) in the copepod Eurytema affinis has been linked with the capacity of this invasive species to move from marine to FW systems (Lee et al. 2011). Ability to cross the osmotic gradient between such FW and SW systems is central to the biogeography, population dynamics, and local adaptation of a variety of organisms.
Freshwater systems are thought to be the ancestral habitat for the Culicidae family of mosquitoes, with ~95% of species inhabiting FW, but the lack of phylogenetic clustering of species occupying SW suggests that there have been multiple acquisitions (and/or losses) of salinity tolerance (Bradley 1987; Clements 1992). Both SW and FW species co-occur repeatedly in different species complexes across genus Anopheles, including 3 non-sister species tolerant of SW (An. merus, An. melas, and An. bwambae) out of a group of 8 morphologically indistinguishable species in the An. gambiae complex (Fontaine et al. 2015; White et al. 2011). Distinct molecular mechanisms for osmoregulation may apply, even among the three SW tolerant sibling species, as An. bwambae requires high mineral content to complete normal development (Harbach et al. 1997; White et al. 2011), while both An. merus and An. melas larvae can develop in either fresh- or brackish waters. Determining the role of salinity tolerance in local adaptation and ecological speciation of anophelines requires detailed phenotypic and molecular analyses of the trait.
Physiological studies suggest a key role of protein localization in anopheline osmoregulation. In FW, mosquitoes absorb ions from the external medium with their anal papillae, and in the rectum further resorption of ions takes place (Bradley 2008). In tolerant anophelines, SW exposure is mediated by osmoregulation and ion excretion (Bradley 2008; Smith et al. 2008). With antibody stains, Linser and colleagues demonstrated that the location of expression of Na+/K+-ATPase of FW species does not change regardless of water salinity, while in SW-tolerant anophelines (An. albimanus, An. merus) it is up-regulated in the dorsal anterior rectum (DAR) and down-regulated in the non-DAR (Smith et al. 2010; Smith et al. 2008; White et al. 2013). The authors propose that this disrupts coordination of Na+/K+-ATPase with V-ATPase, ultimately transforming the role of non-DAR cells from ion resorption to production of hyper-osmotic urine. However, ability of the larvae to shift Na+/K+-ATPase and survive transfers from FW to saline media declines after ~1 d, with >80% of 48 h-old An. merus larvae placed in SW failing to reach pupation (Smith et al. 2010; Smith et al. 2008; White et al. 2013).
Research on the genetic basis of mosquito osmoregulation is in its infancy. Rheault and colleagues (2007) describe two mechanisms mediating ion transport, the first consisting of electroneutral ion exchange (e.g., by Na+/H+-exchangers) energized via ATP hydrolysis, as by Na+/K+-ATPase. Second, Na+/nH+ antiporters can use voltage gradients established by H+-ATPases to influence ion transport (Rheault et al. 2007). Bioinformatic and physiological analyses have suggested the presence of both mechanisms in anophelines, with reports of anopheline electrophoretic antiporters AgNHA1 and AgNHA2 (Rheault et al. 2007; Xiang et al. 2012), as well as anion exchangers such as AgAE2 (Linser et al. 2012). Likely, an integrated network of signaling systems regulates their expression. In fish, osmoregulation is proposed to result from the combined action of osmosensors that detect osmotic stress and signal transduction pathways that respond, along with effectors such as ion transporters that act to restore homeostasis (Kültz 2012). Mitogen-activated protein kinase (MAPK) pathways are associated with osmotic stress signal transduction in taxa ranging from fish species (Whitehead et al. 2013) to plants (Cowan & Storey 2003), and the high-osmolarity glycerol (HOG) pathway is a MAPK signaling pathway associated with osmotic stress in yeast (Li et al. 2012). Thus similar themes in osmoregulatory responses are seen throughout diverse organisms, and may be shared by mosquitoes.
To further understanding of anopheline SW tolerance, we began with QTL mapping from a cross of euryhaline An. merus and stenohaline An. coluzzii, using 2387 markers distributed genome-wide (Smith et al. 2015). We identified six genomic regions associated with osmoregulation: one on chromosome arm 2R, one encompassing the centromere of chromosome 2, two on arm 3R, and two on the X chromosome. Our results suggested the potential for multiple interacting loci, and perhaps different yet fungible combinations of alleles, to contribute to SW survival. However, these six QTL intervals were too large to allow inference of specific genes involved in the trait. Here we expand these findings with RNA-Seq to associate mRNA expression levels of individual genes with SW exposure. We hypothesize that inducible gene expression is associated with osmoregulation, influencing the overall greater tolerance of SW by An. merus compared to An. coluzzii. To address this hypothesis we compare transcriptomic responses of young (18 h post-hatch) An. merus and An. coluzzii transferred for 6 h from FW into one of six solutions, (FW, 10, 20, 30, 40, and 50% the salinity of seawater). We also investigated the response of older (66 h) larvae, which cannot survive a transition to SW (White et al. 2013).
Colonies of the Mali-NIH strain of Anopheles coluzzii, formerly known as An. gambiae “M” form (Coetzee et al. 2013), and the An. merus MAF strain were maintained at the University of Notre Dame insectary as in (White et al. 2013), with minor modifications. Strains were obtained from MR4 (MRA-1156 and MRA-860, MR4, ATCC® Manassas, Virginia). To minimize maternal effects from prior SW exposure, larvae of both the SW-tolerant species An. merus as well as the FW species An. coluzzii were reared in FW (reverse-osmosis water) for ≥2 generations prior to the experiment. Except during the 6 h experimental exposures, larvae were fed daily until pupation with 2:1 w/w ground mixture of Koi fish food to active dry yeast. Saltwater used in experimental assays consisted of sea salts (~98.6% NaCl). In an initial experiment in 2012, species-specific near-lethal salinity levels were chosen for SW concentrations based on the chronic exposure bioassays of White et al. (2013). Survivorship through pupation declines at ~20% the salinity of seawater for An. coluzzii, while for An. merus survivorship declines at ~60% (White et al. 2013). Thus for An. merus the salinity was set slightly below that inducing substantial mortality, at 50% salinity (15.85 g/L sea salts, defining seawater as 31.7 g/L salt), and for An. coluzzii we used 20% salinity (6.34 g/L). Subsequently in 2014, we repeated the experiment exposing each species to a gradient of 0%, 10%, 20%, 30%, 40%, and 50% the salinity of seawater, as detailed below.
In both years, experiments were performed separately but simultaneously on An. merus and An. coluzzii (Fig. 1). Cages of adults received blood meals to stimulate egg production. Synchronous hatching was induced in FW, following (White et al. 2013). This procedure requires a very brief (~6 min.) exposure of eggs to air, and hatching in FW rather than beginning experimental exposure in SW prevents uncontrolled increases in salinity from evaporative water loss. Neonates were placed in small plastic containers (~12.5 cm2) at densities of 50 larvae / 200 mL FW; all 50 larvae were kept together for the remainder of the experiment as a pooled replicate sample. At 18 h post-hatch (approximately 6–7 a.m.), experimental exposures were begun. All larvae in each replicate container were transferred to a new container with 200 mL of one of six solutions: water of 0%, 10%, 20%, 30%, 40%, or 50% salinity. Exposure lasted 6 h, after which time ~35–50 animals from each container were pooled into 1.5 mL eppendorf tubes, snap-frozen in liquid nitrogen, and stored at −80°C until RNA extraction. In 2012 the experiment only tested FW vs. a single saline water concentration of 50% for An. merus, and 20% for An. coluzzii. However in addition to 18 h larvae, the identical protocol was performed on sibling larvae that remained in FW until 66 h post-hatch, and then were transferred to FW or SW for 6 h. With two salinities (FW or SW), two age classes of exposure, and two species, we had eight separate treatments replicated four times apiece, for a total of 32 samples in 2012. In the 2014 experiment we only assayed 18 h larvae, but used the gradient of salinities of 0–50% in 10% increments. Thus with two species, six salinities, and three replicates (12 treatments), we produced 36 samples in 2014.
No mortality was observed for An. coluzzii at salinities <50% during the short 6 h exposure period, nor for An. merus at any salinity. Yet roughly 1/3rd of the animals in the An. coluzzii samples at 50% salinity died. Due to the impossibility of separating the tiny live and dead first instar larvae quickly enough to avoid prolonging exposure beyond 6 h, all An. coluzzii larvae were preserved for library preparation. As discussed below, results for An. coluzzii at salinities ≥30% were treated with caution as some of this gene expression data may reflect impending mortality.
All samples of ~35–50 pooled larvae were used to construct separate mRNA-Seq libraries. Total RNA was extracted and DNA removed with the RNeasy® Mini Kit and RNase-Free DNAse Set following the manufacturer’s instructions (QIAGEN). Construction and barcoding of libraries was performed with the TruSeq® RNA Sample Preparation kit v2 (Illumina) from 600 ng of total RNA per sample. For sequencing we pooled one replicate of each of the eight treatments per lane, with a total of four lanes for the 32 libraries in 2012. Pooled libraries were submitted to the Princeton University Sequencing Core Facility for sequencing on the Illumina HiSeq (V3 chemistry), with 101 base pair single read sequencing. This protocol was repeated in 2014 with all 36 libraries made by the same individual (HAU) using the same eight barcodes. However, HiSeq sequencing was performed by the Beijing Genomics Institute at the University of California Davis facility on five lanes, with 101 base pair paired read sequencing.
Preliminary filtering included application of Illumina’s default Chastity filter, demultiplexing, and visualization using either the FASTX-Toolkit (http://hannonlab.cshl.edu/fastx_toolkit/index.html) or FASTQC (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/) to ensure reads were of high quality. The program Trimmomatic (Lohse et al. 2012) was used for trimming adapters and filtering based on quality scores (default settings). Reads were then mapped to the An. gambiae PEST genome with the software Bowtie2 (Langmead & Salzberg 2012) and TopHat2 (Kim et al. 2013) using the most current geneset downloaded from VectorBase (Giraldo-Calderón et al. 2015), consisting of AgamP3.7 for the 2012 experiment and AgamP4.2 for the 2014 data. Numbers of uniquely mapped reads were counted with HTSeq using the default “union” setting and excluding multi-mapped reads (Anders et al. 2015). We obtained a mean (±1 SE) of ≥6.26 (±0.86) million uniquely mapped reads per treatment in 2012, and a mean of ≥14.0 (±0.6) million uniquely mapped reads per treatment in 2014.
Mapping An. merus reads to the An. gambiae PEST genome could have missed detection of some genes, although these species are very closely related. However we chose this reference genome because the An. merus reference assembly was not well annotated at the time of analysis. Because of the difference in the year of the experiment and the use of paired end data in 2014, statistical analyses of differential gene expression for the 2012 and 2014 data were performed separately. We focus our study on the more comprehensive 2014 salinity experiment, unless otherwise noted. However, 66 h larvae were only assayed in 2012, and hence analysis of ontogenetic effects on SW tolerance is derived from this dataset. Overall patterns and major conclusions for treatments tested in either year are similar for both years. In both 2012 and 2014, An. merus at 50% salinity differentially expresses far more genes than An. coluzzii at 20% salinity. Statistically there is little evidence for an interaction effect of year and water type on gene expression. For An. coluzzii, merely seven genes were significant for a two-way interaction effect of year and water salinity (20% vs. 0%). For An. merus, only four genes were significant for the interaction of year and water salinity (50% vs. 0%).
Finally, raw library counts were normalized with the TMM approach (Robinson & Oshlack 2010) and tested for differential expression using generalized linear model likelihood ratio tests with edgeR (McCarthy et al. 2012). Specific contrasts tested are in Table S1, S2 (Supporting information). We only retained genes with >1 count per million in the same number of libraries as we had replicates, thereby filtering out genes with extremely low expression. This resulted in analysis of 10,905 genes in 2014 and 11,025 in 2012. We refined our list of significantly differentially expressed genes by focusing on those with FDR <0.05 that also showed a log2 fold change <-2 or >2 between SW and FW. Although small-effect changes missed by applying the fold change threshold may be important, we assumed most of the critical genes whose induced expression influences saltwater tolerance would show larger changes. For An. coluzzii, log2 fold change values ranged from −10.47 to 9.44, while for An. merus values ranged from −6.96 to 7.98.
We used DAVID (Huang et al. 2008) to cluster genes in similar functional categories (default settings), focusing on clusters of ≥5 genes that were significantly enriched (enrichment score >1.3, corresponding to a P-value <0.05). This tool incorporates information from multiple functional annotations to group similar genes, and determine enrichment relative to a background set. Here, DAVID was used to identify functional enrichment of differentially expressed genes, relative to a background of all genes passing our filters for low expression. Manual scans of VectorBase gene descriptions for the members of each significant DAVID cluster were used to confirm functional identity of the cluster.
Based on DAVID clusters, five key functional groups emerged. We assessed the norm of reaction for each functional category by computing the number of genes up- and down-regulated across salinities. We used Gene Ontology (GO) membership as reported in VectorBase to identify membership in these five groups. First we designated all genes with GO terms of transport, carrier, channel, antiporter, symporter, or exchanger as transport genes. Second, a list generated by Rob Waterhouse was used to identify immune-related genes, as in (Cassone et al. 2014; Waterhouse et al. 2013). Third, a combined list of potential cytoskeleton and cuticle components was made of GO terms including cuticle, microtubule, cytoskeleton, cytoskeletal, chitin, and chitinase. Fourth, a list of cell cycle genes was made of GO terms including cell cycle, meiosis, meiotic, mitotic, mitosis, cell division, anaphase, checkpoint, cytokinesis, and cyclin (metaphase, telophase, and interphase either were not found, or were combined with other terms). Fifth, a list of sensory genes was made of genes with GO terms of odorant, olfactory, and sensory. Although assignment of gene functional categories from homology approaches such as this is inevitably somewhat subjective, it permitted a coarse assessment of how not only numbers but also functional categories of genes expressed changes across a salinity gradient.
Construction of co-expression modules for network analysis was performed with WGCNA version 1.34 (Langfelder & Horvath 2008). This analysis grouped genes into modules based on their showing similar patterns of increasing or decreasing expression across treatments (species identity as An. coluzzii or An. merus, by water salinity as FW or SW). Data was pooled from both 2012 and 2014 experiments to increase the number of samples, which can increase reliability of such network analyses. Prior to implementing WGCNA the TMM normalized data was log transformed, and effects of year of the experiment were removed using the combat algorithm of the SVA package v. 3.12.0 (Leek et al. 2012). Only the young 18 h data from 2012 was used for network analysis, to avoid the influence of age on gene expression overwhelming differences from salinity and species identity. Large transcriptomic changes are known to occur during larval development (Cassone et al. 2014; Koutsos et al. 2007). Prior work has shown that roughly half of An. coluzzii larvae do not survive to pupation at salinities of 20% and most die at ≥30% salinity; in contrast An. merus survives salinities from 0–50% (White et al. 2013). Thus we questioned whether data from higher salinities for An. coluzzii would be representative of typical gene expression networks, as it is not known how long after exposure to high salinities the expression data reflects impending death. Hence salinities of 30–50% were excluded for An. coluzzii.
Co-expression signed networks were made of positively correlated genes, using data for all genes except those removed due to low expression (genes with ≤1 count per million in the same number of libraries as we had replicates). Thus networks include differentially expressed genes as well as those that do not change across salinities. Construction of networks was based on non-parametric Spearman rank correlations of gene expression levels; other parameters were kept at default settings. We used a soft-thresholding power of 13, yielding an r2 of 0.829 in plots of scale-free topology model fit. Log2 fold change values based on the 2014 dataset alone were plotted for all genes in each module using the “matplot” function in R. These plots provided visual confirmation that the identified modules did represent groups of genes with similar expression shifts across salinities.
Additionally, we assigned regulatory impact factors (RIF) to genes in a differential co-expression analysis, implementing the method of (Hudson et al. 2012) with a custom R script adapting it to RNA-Seq data (Script S1, Supporting information). The RIF statistic is reported as a z-score that ranks genes with potential regulatory roles, based on the difference between two conditions in the regulatory gene’s correlation to differentially expressed genes. This correlation is weighted by the average abundance of the differentially expressed genes, and by the difference in their level of expression (here, TMM normalized mRNA counts) between the two conditions. Because we anticipate transcriptomic networks vary between species, we focus on the difference in Spearman correlation between them, with the two species identities as our two conditions. Similar to our concern in the WGCNA networks, we did not want any subversion of typical gene interactions due to the mortality of An. coluzzii at high salinity. Yet the RIF analysis focuses on differential correlations between species, so we omitted data from salinities ≥30% for both An. merus and An. coluzzii. Expression levels are taken from the 2014 dataset, averaged across 0%, 10%, and 20% salinity.
The RIF analysis was designed to test for regulatory roles of genes such as transcription factors, which may correlate with and influence differentially expressed genes even if they themselves are not differentially expressed. Here we employed RIF twice. First in following the original design of Hudson et al. (2012), we generated a list of potential transcription factors using all genes in the AgamP4.2 geneset with a Gene Ontology category including the word “transcription.” This yielded 416 genes, and is notably a liberal list of possible transcription factors and other genes that may influence transcription. Second, we applied the RIF analysis to test for a relationship between differentially expressed genes, and genes we previously localized to QTL peak regions (Smith et al. 2015). These QTL peaks were identified by mapping genomic loci associated with saltwater tolerance in backcross larvae, which were generated by mating parental An. merus with F1 hybrids of An. merus and An. coluzzii. In this second formulation of the RIF analysis, the transcription factor gene list was replaced with the genes within all six QTL peaks, based on the genomic coordinates reported in (Smith et al. 2015).
Data support our overall hypothesis that inducible gene expression occurs in response to transfer from FW to SW (Fig. 2, Table S1–S2, Supporting information). We assessed the reaction norm of gene expression for young larvae of the SW tolerant An. merus and the obligate FW species An. coluzzii transferred from FW (0%), back into FW or to one of five salt levels: 10%, 20%, 30%, 40%, and 50% the salinity of seawater. Comparing exposure of 10% salinity to FW, An. merus significantly down-regulates 10 genes in saltwater and up-regulates two genes, while An. coluzzii only differentially expresses 3 genes, all down-regulated in SW. At all other salinities, An. coluzzii shows more differential expression than An. merus (Fig. 2).
The greater induction by An. merus than its sibling species at 10% salinity could be important, revealing greater sensitivity to small increases in ion concentrations. Yet the species only exhibited differential expression of nine more genes than An. coluzzii. Moreover in a two-way test of genes with differential induction between species (interaction of species and water salinity), only two genes were significant (AGAP002331, a sodium-independent sulfate anion transporter, and AGAP013746, a gene of unknown function with no GO terms or description in VectorBase at the time of writing). Each species induces a distinct set of genes, sometimes with little overlap (Fig. 2). Indeed at 10% salinity there is no overlap in genes differentially expressed by the two species. Anopheles merus differentially expresses some genes between SW and FW that do not change in An. coluzzii, and vice versa, for all salinities tested.
Analysis of genes uniquely induced by one species reveals a few interesting genes. Allantoinase (AGAP000239) shows inconsistent results in An. coluzzii with up-regulation at 30% and down-regulation at 50%, whereas in An. merus it is up-regulated in all salt levels ≥10%. Allantoinase is involved in nitrogenous waste processing in mosquitoes (Esquivel et al. 2014). This response may reflect increased nitrogen processing from catabolism, as energy is consumed by the SW response. Another gene of interest is the Solute carrier family 4 member 1 (AGAP006968), designated as AgAE2 in the Slc-4-like anion transport family (Linser et al. 2012). We find that it is up-regulated by An. merus at 20–40% SW, and also up-regulated at 10% and 50% although it did not meet our fold change threshold (FDR <0.05, but log2 fold changes of ~1.98 and ~1.67, respectively). The gene is neither up- nor down-regulated by An. coluzzii in any SW vs. FW comparison. Considering this induction by An. merus and lack of induction in the FW species An. coluzzii, AgAE2 represents a strong candidate for an ion transport effector gene that contributes to SW tolerance in An. merus.
For both species, the number of up- and down-regulated genes tends to increase as salinity increases. Neither species shows strong induction until higher salinities, ≥30% for An. coluzzii and 50% for An. merus (Fig. 2, Fig. 3). Rather than linear, the reaction norm of gene expression vs. salinity stress appears closer to exponential. For An. merus the most striking change is seen between the 158 genes differentially expressed at 40% salinity, almost tripling to reach 464 genes differentially expressed at 50% salinity. Down-regulation consistently exceeds up-regulation of expression across species and salinities.
Spearman rank correlations were computed to relate gene expression (averaged across replicate libraries) to salinity levels. Correlations were computed separately for each species (Table S3, Supporting information). Of those genes with very strong correlations (rho >0.8, P <0.05), 15 are correlated positively to salinity for An. merus and negatively for An. coluzzii, while another 15 genes show the opposite relations to salinity (positive for An. coluzzii, negative for An. merus) (Table 1). Of particular interest is GTPase Kras (AGAP002219), a gene correlated positively to salinity for An. merus and negatively for An. coluzzii. This gene is an upstream part of the ERK signaling pathway (Horton et al. 2011), which is one of the three primary MAPK families (together with the JNK and p38 pathways) responsible for transmitting external signals and stressors into cellular responses (Cowan & Storey 2003). The ERK pathway is known for roles in development as well as immunity, and in mosquitoes can be induced by the Plasmodium parasite (Cowan & Storey 2003; Horton et al. 2011). Conversely, the eukaryotic translation initiation factor 3 subunit E (AGAP006944) is correlated negatively with salinity in An. merus and positively in An. coluzzii. Decreasing expression of the eukaryotic translation initiation factor 3 subunit E gene in An. merus, could reflect a tradeoff towards down-regulating protein production for growth, until osmotic homeostasis is restored.
To gain additional insight into transcriptional responses, we focused on the induction of genes potentially involved in similar processes, aided by the DAVID clustering algorithm. Clusters identified by DAVID were summarized into five functional categories: immunity, transport, sensory, cell cycle, and cuticle/cytoskeleton (Fig. 4, Table 2, Table S4–S5, Supporting information). At 10% salinity, insufficient numbers of genes are induced to recover functional clusters in either species. At 20%, both species induce genes associated with immune stress, such as caspase, prophenoloxidase, and fibrinogen. For An. merus, transport represents the category most strongly differentially expressed at low salinity. More transport genes are up- than down-regulated at and above 20% salinity, in contrast to the overall trend of greater gene down-regulation in SW. For sensory genes such as gustatory and odorant receptors, both species show more up- than down-regulation, although sensory genes are not induced by An. merus until ≥40% salinity. An. merus cell cycle genes are not induced until 50% salinity, a sub-lethal stress level for An. merus. In An. coluzzii, cell cycle genes are not differentially expressed until lethal salinity concentrations are reached (≥30%). The cuticle/cytoskeleton genes are not robustly induced in An. merus until salinities reach 40–50%, with a slightly earlier induction at ~30% for An. coluzzii.
Overall, induction of these functions is not unexpected. Transporters are expected to show induction, with a role in mobilizing ions for excretion, and potentially mobilizing energy reserves. Cuticle composition may be altered to combat osmoregulatory stress by decreasing permeability in SW, although for anophelines this strategy has been better studied in the adult stage (Reidenbach et al. 2014). Development is slowed upon saltwater exposure in both species (White et al. 2013), potentially indicating a life-history tradeoff between growth and stress tolerance and explaining impacts on the cell cycle. It is possible that sensory receptors are utilized in detecting changes in ionic composition or heightening sensitivity in the face of environmental stress. Osmosensing is thought to be an important part of the salinity stress response in fishes (Kültz 2012).
Another intriguing category induced by salt stress is the immune system. Prior evidence exists for cross-talk between immune and other stress-related genes in dipterans (Davies et al. 2012). In mosquitoes, the innate immune system responds to infection and can decrease survival of Plasmodium parasites within the insect (Clayton et al. 2014; Garver et al. 2013). Our results suggest the possibility for salinity stress to inhibit the immune system by down-regulating immunity genes, which would be consistent with a diversion of resources from combatting internal pathogens to maintaining homeostasis under environmental stress. It is not known whether the apparent down-regulation of immunity genes would continue to the adult form, when Plasmodium infection can occur. If it persists, larval stress exposure could have important implications for malaria transmission to humans.
We also performed a weighted correlation network analysis (WGCNA) on all genes tested for differential expression, using both the 2014 and 2012 datasets. This allowed clustering of genes that showed similar expression level changes across species and salinity concentrations (Fig. 5, Table S6, Supporting information). The sixteen modules of co-expressed genes identified by WGCNA were assigned a unique color and assessed for the significance of the relationship to salinity and species. Most modules (N =12) were significantly related to salinity (P <0.05), with 9 remaining significant if a Bonferroni-corrected α of 0.003 (=0.05/16) was imposed. For this analysis, modules with positive correlations for species contain genes with higher expression in An. merus; those with positive correlations for salinity contain genes whose expression increases with increasing salinity (Fig. 5).
The WGCNA analysis revealed several co-expression modules of interest. The brown module is even more strongly correlated to species identity (rho =0.85, P <0.001) than it is to salinity (rho =0.35, P <0.02). This module contains the Na+/H+ antiporter AgNHA1 (AGAP002093). Another module of particular interest is the yellow module (correlation to species: rho =0.4, P =0.008; salinity: rho =0.91, P <0.001). This module (Fig. 6) contains several osmoregulatory candidates, including AgAE2 and three potential members of a MAPK phosphorylation signaling pathway: AgJNKa (AGAP009460), Jun (AGAP006386), Puc (AGAP004353). This module also includes Na+/K+-ATPase (AGAP002858), and another Na+/H+ antiporter, AgNHA2 (AGAP002324). The yellow module is particularly intriguing considering the proposed role of electrophoretic antiporters and a shift in Na+/K+-ATPase localization in anopheline osmoregulation (Rheault et al. 2007; Smith et al. 2008; Xiang et al. 2012). Co-occurrence of several candidates, and the module’s overall correlation to salinity and species identity, suggests that genes in this module may function in an overall network of signals (e.g., JNK phosphorylation) and effectors (e.g., AgNHA2) facilitating the survival of An. merus in saltwater. The particularly strong induction of allantoinase (AGAP000293) by An. merus could reflect additional catabolism to provide energy for osmoregulatory processes. Alternatively induction of allantoinase may be a byproduct of an overall up-regulation of excretory processes, including nitrogen excretion, to maintain osmotic homeostasis in saline waters. While postprandial nitrogen excretion has been extensively studied in adult females (Coast et al. 2005; Esquivel et al. 2014), further research is needed on the larval stage of anophelines.
Intramodular connectivity scores from WGCNA gradually declined, rather than being divided into clear groups of either highly connected “hub” genes or non-hub genes with low connectivity. Although our data do not provide strong support for the existence of hubs, one of the most highly connected genes (27th of 1720 genes) in the yellow module is of interest: AGAP010837, identified as MAP4K1 (Horton et al. 2011). This gene is an ortholog of the Drosophila melanogaster gene happyhour, and based on its known role in other systems, Horton et al. (2011) suggest that this kinase may ultimately activate JNK and p38 MAPK signaling. The JNK pathway is our strongest candidate for a phosphorylation signaling cascade induced by SW, as discussed below.
Although little is known of the function of AgJNKa, one can speculate on potential pathways based on those predicted for its neighboring paralog AgJNKb (AGAP009461, not differentially expressed). JNK protein phosphorylates the transcription factor Jun (encoded by AGAP006386, an ortholog of the D. melanogaster Jun-related antigen or Jra) (Garver et al. 2013; Ríos-Barrera & Riesgo-Escovar 2013). However this activity can be halted by the phosphatase Puc, which dephosphorylates JNK in a negative feedback loop (encoded by AGAP004353, an ortholog of the D. melanogaster gene Puckered) (Garver et al. 2013; Martín-Blanco et al. 1998; Ríos-Barrera & Riesgo-Escovar 2013). Up-regulation of all three genes by An. merus in SW genes suggests the JNK signaling pathway may be induced by osmotic stress by An. merus.
In mammals, JNK contributes to regulating mRNA expression of Na+/K+-ATPase subunits (Capasso et al. 2004; Vadász et al. 2012; Wang et al. 2007). Considering the importance of Na+/K+-ATPase protein localization to mosquito SW tolerance, a similar regulatory role of the anopheline JNK is plausible. In general c-Jun N-terminal kinases interact with a variety of transcription factors such as p53, Jun, and NFAT molecules, and were first described as stress-activated protein kinases that mediate responses to several environmental stressors (e.g., ultraviolet radiation and oxidative stress in Drosophila) (Bogoyevitch & Kobe 2006; Ríos-Barrera & Riesgo-Escovar 2013). Thus it is possible that up-regulation of AgJNKa and the downstream transcription factor Jun contributes to regulation of gene expression by An. merus in response to salinity.
The RIF statistic assigns a score to a list of input genes, based on the difference between mosquito species in the correlation of those genes to differentially expressed genes. Even if a given input gene does not itself change expression in SW, a strong correlation with genes that are differentially expressed suggests the potential for it to serve as a regulator of the transcriptomic response to salinity. When using a predicted list of transcription factors as input genes and analyzing salinities of 0–20%, 16 received high RIF z-scores >2 (Table 3, Table S7A, Supporting information). Among these is POU domain transcription factor, class 3 (AGAP005878). This gene has the 10th-highest RIF score (Table S7A, Supporting information), with RIF =2.99. The POU homeobox family contains transcription factors, at least some of which are expressed in osmoregulatory tissues and proposed to regulate osmoregulatory genes (Bürglin & Ruvkun 2001; Wang et al. 2012).
Besides ranking transcription factors with the RIF, we also were interested in exploring the relationship between genomic regions with sequence differences associated with SW tolerance, and gene expression changes. Previously we backcrossed hybrids of An. merus and An. coluzzii to parental An. merus, and related SW tolerance of the backcross progeny to genomic sequences at thousands of loci across the genome (Smith et al. 2015). Based on this QTL mapping experiment, we identified six genomic locations associated with SW tolerance, which altogether encompass 2271 genes. Thus we ran the RIF analysis a second time, using as the input gene list all 2271 genes located within these six QTL. With this implementation, we found 125 genes in the QTL regions with high RIF z-scores of >2 (or <-2). Among these are mitogen-activated protein kinase kinase kinase 7-interacting protein 1 or MAP3K7IP1 (AGAP002953) in the 2Rop inversion with RIF=4.01 and rank 20 (Table S7B, Supporting information), and mitogen-activated protein kinase kinase 3 or MAP2K3 (AGAP000310) on the X chromosome with RIF=3.22 and rank 38. According to the kinome depicted by (Horton et al. 2011), MAP3K7 is part of the JNK pathway, while MAP2K3 is upstream in the p38 pathway.
Older larvae of An. merus alter expression of more genes than younger larvae, in the single-factor tests of response to 50% SW (423 vs. 351 genes, respectively) (Table S1, Supporting information). Functional clustering of differentially expressed genes reveals similar functional classes between age groups (Table S5, S8, Supporting information), including sensory genes (e.g., odorant and gustatory receptors), and immune genes (e.g., fibrinogens).
The knowledge that larval tolerance of a transfer to SW diminishes at ages of approximately ≥24 h (White et al. 2013) prompted us to examine candidate genes (Table S9, Supporting information) for patterns of greater expression at the younger vs. older age class. For example, expression of AgJNKa and the Jun transcription factor is high in both age classes. Thus, while the JNK pathway may play a role in the osmotic stress response, it likely is not responsible for the unique ability of young An. merus to survive high salinity. However, comparison of age classes reveals AgAE2 to be a strong candidate for the SW tolerance of young larvae. Not only is expression of AgAE2 consistently higher in SW between An. merus and An. coluzzii, but also the expression of this gene barely changes between FW and SW in older An. merus (a mere log2 fold of ~0.14) (Table S9, Supporting information).
By assaying differential expression of larval An. coluzzii and An. merus in SW and FW, we supported our hypothesis that salinity exposure induces differential gene expression, which may in turn influence their different SW tolerances and habitat occupancy. One of the greatest shifts in gene expression was seen for AgJNKa, a mosquito-specific MAPK, while AgAE2 represents one of the most promising ion transport candidate genes up-regulated by young but not old An. merus. Future knockdown or gene silencing assays such as CRISPR/Cas (Sander & Joung 2014, Bono et al. 2015) should be employed to confirm the roles of candidate genes in osmoregulation. It also will be important to investigate tissue-specific changes in gene expression; our use of whole bodies may miss genes that change expression only in one tissue, or change expression in opposite directions among different tissues. Field research is needed to explore how transcriptomic shifts during salinity transitions differ in a natural setting where other biotic (e.g., predation, competition, food resources) and abiotic (e.g., thermal fluctuations, hypoxia, pollution) forces may modulate the response. Bradley (2008) suggests that osmoregulation is an energetically costly strategy for tolerating salinity, but viable in environments with abundant food resources. Once successful laboratory colonies are established for additional members of the An. gambiae complex, it will be interesting to compare transcriptomic responses to a salinity gradient in other SW-tolerant species and explore the potential for adaptive plasticity vs. constitutive expression. Indeed the effect of age on the ability of An. merus to survive SW has been interpreted as a form of developmental plasticity by White et al. (2013).
Despite its limitations, the controlled laboratory system employed here allows for a crucial first step in developing hypotheses of candidate osmoregulatory genes and pathways. Overlap between genes responding to salinity with those associated with immunity reveals that such studies also may provide important insights regarding vectorial capacity, and perhaps differences in malaria transmission between species or habitat types. With climate change and extreme weather events creating the potential for saltwater intrusion and salinization of inland waters (Ramasamy & Surendran 2011), improved understanding of SW-tolerant disease vectors may be particularly timely.
We gratefully acknowledge the help of Luke Xiao and Marcy Kern in the insectary, the Notre Dame Genomics Core for quality tests during library construction, Paul Howell and MR4 for anopheline eggs (Anopheles gambiae Mali-NIH M form, MRA-860, deposited by T. Lehmann NJB, and Anopheles merus MAF, MRA-1156, deposited by M. Coetzee), and Dr. Bradley White for discussions that aided project conceptualization. We appreciate the insight of Dr. Antonio Reverter regarding RIFs, and comments by reviewers that improved this manuscript. Funding was provided by National Institutes of Health award 1R21AI101459 to NJB.
The views and opinions expressed by the authors do not necessarily reflect those of the agencies or institutions with which they are affiliated.
Expression levels (raw counts from HTSeq) for each gene identified by their VectorBase gene ID (AGAP + 6 digits for gene identifier), along with links to fastq files, are available from NCBI GEO reference series GSE74820.
Author ContributionsThis project was designed by NJB and HAU. Larval water exposure, library construction, and statistical analyses were conducted by HAU; RIF analyses were performed by HAU and CC. The manuscript was drafted by HAU, and edited and approved by NJB, CC, and HAU.