Analysis of natural female post-mating responses of Anopheles gambiae and Anopheles coluzzii unravels similarities and differences in their reproductive ecology

Anopheles gambiae and An. coluzzii, the two most important malaria vectors in sub-Saharan Africa, are recently radiated sibling species that are reproductively isolated even in areas of sympatry. In females from these species, sexual transfer of male accessory gland products, including the steroid hormone 20-hydroxyecdysone (20E), induces vast behavioral, physiological, and transcriptional changes that profoundly shape their post-mating ecology, and that may have contributed to the insurgence of post-mating, prezygotic reproductive barriers. As these barriers can be detected by studying transcriptional changes induced by mating, we set out to analyze the post-mating response of An. gambiae and An. coluzzii females captured in natural mating swarms in Burkina Faso. While the molecular pathways shaping short- and long-term mating-induced changes are largely conserved in females from the two species, we unravel significant inter-specific differences that suggest divergent regulation of key reproductive processes such as egg development, processing of seminal secretion, and mating behavior, that may have played a role in reproductive isolation. Interestingly, a number of these changes occur in genes previously shown to be regulated by the sexual transfer of 20E and may be due to divergent utilization of this steroid hormone in the two species.


Results
Collection of mating couples from natural swarms. In order to analyze the natural post-mating response of females from the two anopheline sibling species, we collected 91 An. gambiae and 75 An. coluzzii mating couples from different swarms in the villages of Soumousso and Vallèe du Kou (Burkina Faso) (Fig. 1). Females of each couple were then dissected at either 1 day or 4 days post mating (PM), to capture the short-term as well as the lasting, long-term response to copulation. We dissected the lower reproductive tract (LRT) comprising atrium and spermatheca, and the rest of the body (carcass).
Virgin females were instead produced by collecting larvae from natural breeding sites, and LRTs and carcasses were dissected from resulting adult females at 2 and 5 days post emergence. Because the age of mated females could not be determined as they were caught in natural mating swarms, we chose these time points for tissue collection in virgins to approximately age-match these samples to the ones dissected from mated females, given that it is generally believed that females mate on the second night after emergence 44 . Post-mating transcriptional response in the lower reproductive tract (LRT) of field An. gambiae and An. coluzzii females. In our analysis of the post-mating response in the LRT, we focused on ten genes shown to be strongly up-or down-regulated after mating in laboratory experiments and thus likely to be involved in the reproductive processes triggered by copulation 18,22 (Tables 1, 2; Fig. 2). These included 9 genes whose function in An. gambiae has not been determined yet -i.e. one ABC transporter (AGAP011518), three serine proteases (AGAP005194, AGAP005195, AGAP005196), one amino protease (AGAP000885), two metallopeptidases (AGAP001791 and AGAP009791), a protease inhibitor (AGAP009766), and a putative anti-microbial Andropin-like gene (AGAP009429) 45 . The last gene was the mating induced stimulator of oogenesis (MISO, AGAP002620), which is induced by the sexual transfer of the steroid hormone 20E and is implicated in the increase in egg development experienced by mated An. gambiae females after blood feeding 17 .
With the exception of AGAP009766 and Andropin-like, which were previously shown to be upregulated after mating, results were consistent with those obtained in the laboratory, showing that the transcriptional response to mating is mostly conserved after colonization 18,22 . Post-mating transcriptional response of genes related to reproduction, mating behavior and immunity in the carcass of field An. gambiae and An. coluzzii females. We next analyzed the expression of genes in the female carcass, initially focusing on seven factors that may be related to reproductive success or mating behavior (Table 1 and Table 2; Fig. 3). Our analysis included Vitellogenin (Vg, AGAP004203), which encodes a yolk protein that is needed for egg development 46 and which in Aedes aegypti mosquitoes is strongly upregulated by 20E synthetized after blood feeding 47 , and other six genes shown to be differentially expressed in An. gambiae and An. coluzzii virgin females 48 that we reasoned might be associated with assortative mating behavior. These included: the sex determining gene doublesex (dsx, AGAP004050), the antennal carrier protein AP-1, (AGAP004799), the odorant binding protein 25 (OBP25, AGAP012320), the cuticular protein CPF3 (AGAP004690), the glutathione S transferases -epsilon class 2 (GST-E2, AGAP009194); and lingerer (AGAP004817). dsx regulates the terminal sexual differentiation in most insects 49 , and specifically it determines the differentiation of neurons that control male courtship behavior 50 and female sexual receptivity [51][52][53] . In Drosophila, lingerer is involved in the control of male copulatory organs during courtship 54 . CPF3 is a non-canonical cuticular protein with no chitin-binding capacity, which may be part of the epicuticle 55 where it could bind to sex pheromones such as cuticular hydrocarbons (CHCs) 48 . Odorant binding proteins such as AP-1 and OBP25 transfer odorants to specific receptors 56 and might play a role in female mate choice by helping the identification of co-specific males. GST-E2 might instead be involved in the metabolism of chemical stimuli from antennae and other sensory organs 57 , thus regulating the availability of stimulants such as CHCs.
We finally tested whether mating induces a differential immune response in the two species, possibly driven by diverging sexually transmitted pathogens 58,59 . To this aim, we evaluated the expression levels in the female carcass of five immunity-related genes: the thioester containing protein 1 (TEP1), which is a complement-like factor, homologous to the human C3, that binds and mediates killing of pathogens including Plasmodium parasites 60 ; the leucine-rich immune protein 1 (LRIM1), which circulates in the hemolymph as a disulphide-bounded complex with the leucine-rich protein APL1C and interacts with TEP1 controlling its activity 61,62 ; and the antimicrobial peptides cecropin 1 (CEC1), CEC3, and gambicin (GAMB) 63,64 . Only TEP1 showed to be upregulated in An. gambiae females at 1 day PM (ANOVA P = 0.0297; post-hoc test with FDR correction P = 0.0044).

Discussion
Our results on the transcriptional response to mating in An. gambiae and An. coluzzii females collected from natural mating swarms largely corroborate previous data obtained under laboratory conditions 17,18,22 , demonstrating the opportunity of studying complex phenomena such as mating and post-mating behavior in laboratory colonies. This result is remarkable when considering that gene expression is age-dependent [65][66][67] and that in our study it was not possible to precisely age-match mated females to virgin ones. For this reason -as well as for the limited number of samples we analyzed due to intrinsic difficulties in collecting couples from natural mating swarms -we observed some variability in our results that probably limited our power to detect subtler, age-dependent changes.
Despite these limitations, some interesting differences were detected in the post-mating responses of the two species samples. Although field An. gambiae and An. coluzzii males and females from the same geographic areas studied here share largely overlapping reproductive microbiomes 42 , we detected a mating-induced regulation of TEP1, a key immune gene, in the carcass of An. gambiae females. This species-specific upregulation may be due to sexual transfer of microorganisms populating the An. gambiae male reproductive tract, similarly to what observed in D. melanogaster where mating anticipates immune reactions to sexually transmitted pathogens possibly as a mechanism to enhance fecundity 68 .
Perhaps more interestingly, our data also highlight differential mating-induced changes in genes involved in oogenesis, which may reflect inter-specific differences in the physiological processes leading to egg development. First, we show that MISO -an atrial gene strongly induced by sexual transfer of 20E that regulates the number of eggs developed by females after mating and blood feeding 17,19 -was significantly upregulated only in An. coluzzii at 1 day PM, although a trend towards an increase was also observed in An. gambiae at the same time point (Tables 1 and 2, Fig. 2). Given that MISO interacts in the atrium with 20E transferred during mating 17 , the differential transcriptional dynamics of this gene in the two species suggests that the timing of release of the steroid hormone from the mating plug may be regulated in a species-specific fashion. Second, we reveal that another 20E-induced gene important for oogenesis, Vg, is differentially regulated in the female carcass of the two anophelines. This yolk protein precursor, produced in the fat body and incorporated in the developing eggs  69 , was strongly upregulated in An. gambiae at 1 day PM. This difference may reflect a reduced reliance of An. coluzzii females on mating for oogenesis, possibly due to an increased ability to store nutritional reserves during larval development 70 , and may provide some cues on why females of this species are competent to start egg development as virgins, while An. gambiae females generally need a mating-induced boost to promote the same process 44,70,71 . Even if the two species have a similar competence for Plasmodium transmission in the laboratory 72,73 , the fact that An. gambiae females often require multiple blood feedings to complete oogenesis 70 may have important implication for malaria transmission in field settings, as it may increase its chances to become infected with Plasmodium parasites earlier in adult life and be associated with higher infection prevalence, as observed in some regions 74,75 .
We also detected differences in the regulation of the atrial proteolytic machinery which may be involved in the digestion of the mating plug and other seminal secretions. While the protease AGAP009791 was significantly repressed in both species at both time points analyzed, other proteases were downregulated in a time-and species-specific manner. Specifically, AGAP001791 was repressed at both time points in An. gambiae but only at 4 days PM in An. coluzzii; at 1 days PM AGAP000885, AGAP005195 and AGAP005196 were downregulated in An. gambiae, with AGAP000885 and AGAP005195 repressed also at 4 days PM in this species, while AGAP005194 levels were downregulated in An. coluzzii at 1 day PM only. As postulated for MISO, the differential expression of proteolytic enzymes is consistent with the occurrence of species-specific timing of digestion of seminal secretions, which is associated with fertility in Drosophila as well as An. gambiae 20,25,26 and, when perturbed, can lead to speciation 23,24 . Intriguingly, several codons -including those close to the catalytic portion -of the genes encoding the atrial proteases AGAP005194, AGAP005195 and AGAP005196 are evolving under long-term and episodic positive selection in the An. gambiae complex 76 , supporting the hypothesis that timely and proper mating plug digestion might drive the emergence of post-mating pre-zygotic barriers in species of this complex. Similar to the activation of MISO and Vg, the post-mating downregulation of the proteolytic machinery appears to depend on the sexual transfer of 20E, as all six proteases analyzed here were repressed in the atrium of virgin laboratory females following 20E injection 18 .
Finally, the expression of four genes encoding for factors possibly involved in mating behavior (dsx, CPF3, AP-1 and OBP25) was regulated by mating in An. gambiae females only. dsx is a key gene in the sexual differentiation cascade, and is produced as sex-specific isoforms 77,78 that in Drosophila govern multiple aspects of reproductive biology, including the female receptivity to mating and the development and the activity of neural circuit that regulate sex-specific sexual behavior [51][52][53] . Furthermore, in the fruit fly dsx controls the expression of genes that synthetize female-specific long-chain cuticular hydrocarbons (CHC), notably the desaturase DESAT-F, that are potent pheromones for male courtship behavior 79 . It is therefore possible that dsx may affect the synthesis of    CHC pheromones also in An. gambiae females, consistent with the observation that An. gambiae and An. coluzzii have, indeed, slightly different CHC profiles that are altered after mating 36,38 . Interestingly, CPF3, another gene related to CHC, is also downregulated in An. gambiae females. This cuticular protein is likely expressed in the epicuticle, where it is postulated to bind to cuticular pheromones such as CHCs 48 . Because post-mating changes in the CHCs profiles affect female attractiveness in many monandrous insect species 80,81 , the An. gambiae-specific downregulation of genes related to chemical contact cues might reflect the occurrence of different post-mating signals in the two species. Interestingly, both cuticular proteins (CPs) and CHCs have been linked to 20E function, as this ecdysteroid reduces CP expression levels during development 82 and is involved in CHC production in adult Drosophila 83 . Although a link between expression of the genes studied here and male-transferred 20E has yet to be confirmed in field setting, the different post-mating regulation of genes shown in laboratory conditions to be controlled by this steroid hormone 17,18,21 supports the hypothesis of divergent male 20E effects in the two species, consistently with the finding of differential 20E levels in the MAGs of An. coluzzii and An. gambiae males in Burkina Faso 84 . It is intriguing to speculate that the transcriptional differences observed here could represent signatures of a divergent evolutionary arms race between the sexes, which in turn may have led to changes in key reproductive processes and possibly to the development of mechanisms of sexual isolation 85,86 .  Virgin females were obtained from larvae collected in natural breeding sites (3 in Soumousso and 6 in Vallèe du Kou) and reared to the adult stage in natural climatic conditions and photoperiod using cages placed in the outdoor space available at the IRSS laboratory in Bobo-Dioulasso. In order to ensure females would not mate, pupae were individually transferred in single cups and their sex determined at emergence. Using this method, adult females and males were never in contact with each other.

Mosquito sample preparation. Anopheles gambiae and
Mating couples were collected from naturally occurring swarms as previously described 33,39,87 , allowed to complete copulation, transferred to single cups using mouth aspirators, and brought to the laboratory in sealed containers avoiding shifts in temperature and humidity.
Virgin and mated females were maintained in individual cups and DNA was extracted from single legs removed from live specimens for genotyping 88 prior to dissections of reproductive organs. These were carried out under a dissecting stereo-microscope (5x magnification lens) at different time intervals, i.e. 2 and 5 days post emergence in the case of virgin females, and 24 hours and 4 days post-mating in the case of mated ones (to analyze both short-term and long-term response to mating). The time points for virgin females were selected to match as much as possible the age of mated females based on data showing that most females mate on the second night after emergence 44 . The lower reproductive tract (LRT, comprising atrium and spermathecae) and the rest of the body (carcass) of single females were stored separately in RNAlater solution (Ambion) and pools of five individual tissues/species/time interval were obtained for each time point (Table S1).

RNA extraction and cDNA synthesis.
For tissue-specific analysis, total RNA was extracted using TRI Reagent (Helena Biosciences). The amount of RNA for female carcasses was limited to 1 μg. All samples were treated with DNase I (Invitrogen), according to manufacturer's guidelines. cDNAs were synthesized in 100 μl reactions using 1x First Strand buffer, 5 mM DDT, 0.5 mM dNTPs, 2.5 μM random hexamers, 40 units RNaseOut recombinant ribonuclease inhibitor, and 125 units of M-MLV reverse transcriptase (all reagents from Invitrogen).
Quantitative Reverse Transcription PCR with SYBR green detection. Samples were run in 15 μl reaction volume using 1x Fast SYBR Green Master Mix (Applied Biosystems). Gene expression was quantified in duplicates on a StepOnePlus Real-Time thermocycler (Applied Biosystems) using the following program: 95 °C for 15 min, then 40 cycles (95 °C for 15 sec, 60 °C for 60 sec) followed by a dissociation curve analysis. Primers  Table S2. Three technical replicates were used for each biological replicate for each gene. A standard curve against serial dilutions of cDNA templates (mated and virgin) was used for each gene to determine the linear range of the assay. Statistical analysis. Gene expression levels were normalized using deltaCt method against the ribosomal gene RpL19 (AGAP004422), which is expressed at high levels and does not respond to mating 18,22 . To test for mating-induced changes in gene expression, the two species were studied separately. As we do not know if the primers anneal with the same efficiency or if the reference gene is expressed at the same levels in the two species, we did not perform cross-species comparisons. An ANOVA was first used to test whether each gene showed significant changes in the two time points analyzed, including in the analysis both mated and virgin samples. If the global F test gave positive results, pairwise posthoc contrast tests (Tukey-Kramer procedure) have been used to determine differences between mated and virgin females at each time point analyzed. To control for possible Type I error arising through use of multiple ANOVA tests, the P values were corrected by applying a False Discovery Rate procedure (FDR).