The urban-adapted underground mosquito Culex pipiens form molestus maintains exogenously influenceable circadian rhythms

Examining the ways in which distinct life-history traits are linked to one another can provide a better understanding of the mechanisms through which organisms adapt to, or are constrained by, environmental conditions (Blows and Hoffmann, 2005; Mauro and Ghalambor, 2020). In insects, two important life-history traits that are tightly linked to environmental variation are circadian rhythms and seasonal dormancy in response to prolonged, adverse climatic conditions (Tauber and Tauber, 1981; Danks, 1987; Kyriacou et al., 2008). Individual variation in both of these traits can have substantial fitness consequences, with mis-timed behavioral responses leading to reduced or absent mating opportunities, reduced access to food resources or increased probability of mortality (Danks, 1987; Yerushalmi and Green, 2009).

The perception of external light cues plays an important role in both the maintenance of circadian rhythms and the initiation of seasonal dormancy (Beck, 1980; Denlinger et al., 2017). Insect circadian rhythms are usually exogenously entrainable, meaning that their timing is determined (‘entrained’) over a 24 h period by the daily succession of light and dark. This direct influence (i.e. entrainability) allows individuals to adjust their behavior to match seasonal or geographic variation in day length (Meireles-Filho and Kyriacou, 2013). A dormancy response to predictable seasonal changes, such as colder periods in temperate climates (i.e. ‘winter’), is also often strongly influenced by changes in day length (Adkisson, 1966). In insects, one of the most common responses to oncoming inclement weather that is seasonally consistent is a prolonged form of dormancy known as diapause. The induction of diapause typically involves many specific metabolic changes, and during the time an insect is in the diapause state, growth, development and reproduction are generally suspended (Tauber and Tauber, 1976). Diapause is common in insects living in temperate regions where the summers are warm with longer days and the winters are cold with shorter days (Denlinger, 2002).

The northern house mosquito, Culex pipiens, is an important disease vector of many arboviruses including West Nile virus and Saint Louis encephalitis (Turell, 2012). For this reason, understanding the nature and manifestation of circadian rhythms and seasonal dormancy in this species has implications for vector population control and disease mitigation (Ewing et al., 2019). In the pipiens form of this mosquito, females enter diapause with the onset of winter, spending this time remaining relatively stationary in sheltered locations (reviewed in Vinogradova, 2000). The ovarian follicles of females grow relatively slowly while they are in a diapause state (Sanburg and Larsen, 1973), and they generally will not seek out a blood meal until spring and the termination of diapause (Eldridge and Bailey, 1979; Mitchell, 1983; Mitchell and Briegel, 1989; Bowen, 1992). Like most diapausing insects, they rely on the lipid reserves that are stored prior to entering diapause as part of its initiation process (Zhou and Miesfeld, 2009). Male Culex pipiens f. pipiens do not diapause, and die with the onset of winter. Culex pipiens f. pipiens populations go through multiple generations each year and, as such, diapause is facultative, being induced predominantly by the shorter photoperiods and lower temperatures experienced by fourth instar larvae and early pupae in late summer and autumn (Eldridge, 1968; Spielman and Wong, 1973).

As both circadian rhythms and the initiation of diapause are closely tied to the perception of external light cues, it is perhaps not surprising that they are correspondingly linked molecularly in Cx. pipiens f. pipiens. Specifically, many well-known circadian clock genes have been shown to also influence the initiation of diapause in this mosquito (Meuti and Denlinger, 2013; Meuti et al., 2015; Chang and Meuti, 2020). For example, using RNA interference (RNAi), Meuti and colleagues (2015) suppressed the expression of the negative circadian regulators period (per), timeless (tim) and cryptochrome2 (cry2), causing females reared in diapause-inducing conditions to continue ovarian follicle development and fail to accumulate lipid reserves, two characteristics indicating the absence of diapause initiation. In a reciprocal experiment, suppressed expression of the gene pigment dispersing factor (pdf) caused females reared in diapause-averting conditions to enter diapause, with corresponding lipid reserve accumulation and the arrested development of ovarian follicles (Meuti et al., 2015). This work convincingly shows the importance of functioning circadian clock genes for the initiation of diapause.

In temperate climates, a second Cx. pipiens, known as the molestus form, is an urban-adapted mosquito that lives predominantly in subterranean locations such as flooded basements, sewers and subway tunnels (Byrne and Nichols, 1999; reviewed in Vinogradova, 2000). The taxonomic status of Cx. pipiens f. molestus remains a challenging biological problem, and many authors treat it either as a subspecies or ecological form of Cx. pipiens (e.g. Vinogradova, 2000; Harbach, 2012). Correspondingly, it has been suggested that Cx. pipiens f. molestus is recently derived (Kothera et al., 2010; Harbach, 2012; Becker et al., 2012; Asgharian et al., 2015), potentially arising with the emergence of human settlements (Fonseca et al., 2004; Gomes et al., 2009; Price and Fonseca, 2015).

Unlike Cx. pipiens f. pipiens, Cx. pipiens f. molestus will readily mate in enclosed spaces (‘stenogamy’; Roubaud, 1930; as cited in Vinogradova, 2000), and has the ability to lay one batch of eggs before requiring a blood meal from a vertebrate host (‘autogeny’; Roubaud, 1929; Spielman, 1971). Additionally, the exogenous entrainment of circadian rhythms via light cues may be absent or diminished in this mosquito, as several previous studies have shown that Cx. pipiens f. molestus exhibits less sensitivity to changes in photoperiod than other Culex taxa (e.g. Chiba et al., 1981; Karpova, 2009; Fritz et al., 2014). Finally, both observational and experimental studies indicate that Cx. pipiens f. molestus lacks the ability to enter a diapause state in response to a shortened photoperiod, and instead remains reproductively active year round (Richards, 1941; Spielman and Wong, 1973; Vinogradova, 2000; Kassim et al., 2013; Nelms et al., 2013; Bajwa and Zuzworsky, 2016).

The underground locations where temperate populations of Cx. pipiens f. molestus predominantly live are unlikely to experience natural fluctuations in photoperiod. It is therefore possible that selection pressures to maintain this mosquito's ability to respond to external photoperiod cues have been reduced or eliminated. This diminished purifying selection could have allowed novel, derived genetic changes to accumulate in genes influencing circadian rhythms. If ancestral populations of the molestus form had the ability to enter diapause, then such derived changes could also have led to the observed absence (loss) of diapause ability in contemporary urban populations of Cx. pipiens f. molestus.

In this study, we experimentally examined photoperiodically entrainable circadian rhythms and the potential for induction of ecologically relevant reproductive dormancy in a Cx. pipiens f. molestus population originating from an underground location in New York City, USA. Our null hypotheses were that these mosquitoes do not alter circadian behaviors in response to variation in photoperiod and that exposure to a shortened photoperiod during larval and pupal development would not induce any detectable changes in adult reproduction. We also investigated the presence of Cx. pipiens f. molestus-specific derived genetic changes that may be associated with an inability to express a photoperiodically induced dormant state. Specifically, we postulated that an inactivating mutation or other major structural variant in one or more clock genes could simultaneously account for its inability to enter diapause and the hypothesized reduction or loss of entrainable circadian rhythms. Finally, we looked for non-synonymous single nucleotide polymorphism (SNPs, aka ‘missense mutations’) in these genes, and used expression data to survey possible regulatory influences in a binary context on diapause and circadian rhythms (i.e. a gene is expressed or not).

We first wanted to assess whether our underground-derived lab colony of Cx. pipiens f. molestus retains the ability to behaviorally respond to exogenous influences, specifically variation in light cues. These mosquitoes originated from adult females collected in a New York City residential basement in December 2010 (Price and Fonseca, 2015). They have been continuously maintained in various labs from that time to the present. Since 2017, we have kept all life stages of this mosquito colony at ambient temperature in a windowless basement room with no fixed light cycle, mimicking their natural habitat in urban environments. Adults are kept in 60 cm³ screened flight cages (BugDorm-2120 Insect Rearing Tent), where they are allowed to mate freely. They are given access to an 8% sucrose solution ad libitum. As Cx. pipiens f. molestus are autogenous, females do not require a blood meal before they will lay eggs. To start a new generation, egg rafts are collected by placing a black disposable food tray (henceforth, ‘oviposition trays’) with 400 ml of dechlorinated water into the cages for 48 h. The rafts are then removed and the eggs are allowed to hatch. For colony maintenance, first and second instar larvae are transferred from the oviposition trays to 600 ml of dechlorinated water in 1 l white plastic disposable food containers. Larvae are fed a diet of TetraMin tropical fish flakes approximately once per week, with the amount varying depending on larvae size and density. Individuals are allowed to pupate in these containers and, upon emergence, adults are transferred to the aforementioned flight cages.

Over three separate trials, we obtained 17 egg rafts from our maintenance colony as described above. After hatching, we placed 5–15 larvae from each egg raft (i.e. full siblings) into a 100 ml glass jar with 50 ml of dechlorinated water. These jars were then kept in a TriTech^TM Research Digitherm^® 38 l Heating/Cooling Incubator, with optional light-blocking door coating and seal (https://www.tritechresearch.com/DT2-MP-47L.html) at one of four environmental conditions: 12 h:12 h light:dark (lights on at 06:00 h, lights off at 18:00 h), 12 h:12 h dark:light (lights off at 06:00 h, lights on at 18:00 h), 24 h of light or 24 h of dark. All four incubators were kept at 23°C and larvae were fed as described above with the amount varying based on larval size and density.

After pupation, individual pupae were transferred to polystyrene ‘wide’ Drosophila vials (Genesee Scientific) that contained approximately 25 ml of dechlorinated water. These vials were then moved to a custom-built dark chamber with a light regime of 12 h:12 h light:dark (lights on at 06:00 h, lights off at 18:00 h). During the dark period, a red light would come on once every hour for 1 min in order to help detect any emergence activity. Adult emergence was filmed using a YI 1080p home security camera with night vision capability. Video technology such as this has been used to record mosquito behavior and activity in laboratory settings previously, including under light–dark and constant dark conditions (e.g. Araujo et al., 2020). We focused on the transition event (eclosion) between the pupal and adult states because it is relatively quick and discrete, and is easily noted in video footage.

We reviewed all the video recordings and documented the times of all emergences that occurred during the trials, as evidenced by observing the mosquito emerging from its pupal case or the first point at which an adult was observed treading on the surface of the water. For the purposes of analysis, we converted minutes into decimals (i.e. an adult emergence at 12:30 h was recorded as 12.50 h). The Rayleigh test was used to test for uniformity of the data (i.e. whether adult emergence occurred at a certain time of the day). Our null hypothesis was that there was no pattern displayed in the adult eclosion of Cx. pipiens f. molestus. To test for differences between treatments, we used the Watson–Williams test of homogeneity of means, performing all pairwise comparisons. All statistical analyses were done in R v.4.0.2 (http://www.R-project.org/), utilizing the ‘circular’ package (https://r-forge.r-project.org/projects/circular/), and statistical significance was considered at P<0.05.

A true diapause state in female Culex mosquitoes is typically assayed by removing the ovaries and measuring the length of the follicles (e.g. Eldridge, 1968; Spielman and Wong, 1973; Sim and Denlinger, 2008; Meuti et al., 2015). All North American Cx. pipiens f. molestus populations examined in this manner have been shown to undergo ‘normal’ ovarian development when reared in short-day conditions (Spielman and Wong, 1973; Sanburg and Larsen, 1973; Nelms et al., 2013). However, some members of the Culex pipiens species complex, such as Culexaustralicus, are known to enter a quiescent state during inclement weather, which, while not true diapause, nonetheless appears to be an adaptation to variable climates (Dobrotworsky and Drummond, 1953; Dobrotworsky, 1967). It is possible that a dormant state more akin to quiescence than to diapause (and distinct from ovarian arrest) could be induced in Cx. pipiens f. molestus in response to photoperiodic cues (Diniz et al., 2017). Therefore, independent of our circadian rhythm examinations, we wanted to confirm that our New York City Cx. pipiens f. molestus line lacks the ability to enter an ecologically relevant dormant state (e.g. ‘quiescence’) in response to a short photoperiod.

To do this, we collected egg rafts from 12 different females from our stock colony using oviposition trays in the manner described above. Each raft was isolated individually in a 100 ml glass jar along with 50 ml of dechlorinated water and 100 µl of a solution made from vigorously blending 100 mg of fish flakes into 10 ml of dechlorinated water. These were kept at 18°C with a 12 h:12 h light:dark (07:00–19:00 h) photoperiod. Then, 24 h after hatching, we split each family into approximately equivalent (numerically) groups that were assigned to one of two photoperiod conditions, ‘long-day’ and ‘short-day’. Larvae placed in long-day conditions were kept on an 18 h:6 h light:dark photoperiod (04:00–22:00 h) at 18°C, whereas larvae in short-day conditions were kept on a 6 h:18 h light:dark photoperiod (10:00–16:00 h), also at 18°C. All egg, larval and pupal environmental conditions were maintained utilizing a TriTech^TM Research Digitherm^® 38 l Heating/Cooling Incubator, with optional light-blocking door coating and seal (https://www.tritechresearch.com/DT2-MP-47L.html). First and second instar larvae were given 200 µl of fish flake feeding solution daily, third instar larvae were given 400 µl daily, and fourth instar larvae were given 800 µl daily.

Upon pupation, individual mosquitoes were placed in a 100 ml glass jar with 50 ml of dechlorinated water. These pupae were maintained in the same environmental conditions as the larvae, depending on the treatment group (long-day versus short-day). When they eclosed (emerged as adults), we moved adult females to a 30 cm³ flight cage (https://shop.bugdorm.com/bugdorm-1-insect-rearing-cage-p-1.html), with virgin males from a separate, dedicated mating pool. Males from this mating pool were maintained in the long-day photoperiod at 20°C from hatching to adulthood. The slightly higher temperature was to ensure that males had time to reach sexual maturity before their interaction with experimental females. Males were at least 3 days old prior to the introduction of experimental females to ensure these males were sexually active (Vinogradova, 2000). Furthermore, we maintained a ratio of two virgin males to each experimental female in each flight cage.

To maximize the likelihood of insemination, we kept females with males reared at 20°C in the flight cages for 72 h. Because of space limitations, the flight cages for both short-day and long-day treatments were maintained on a long-day photoperiod (18 h:6 h light:dark), at ambient temperature (∼25±2°C for the duration of the study). After 72 h, females were placed individually into 100 ml glass jars with 50 ml of dechlorinated water to lay eggs. These jars were also kept in long-day (18 h:6 h light:dark) conditions and at ambient temperature.

We examined multiple traits associated with reproductive activity in females. First, we compared the percentage of females from each treatment that laid eggs within 10 days of being placed in an oviposition jar. For individual females, we also recorded the time to lay eggs (checked every 12 h at 09:00 h and 21:00 h), and the number of eggs laid. Female Cx. pipiens f. molestus typically lay eggs within 4–5 days of eclosion when maintained at 25°C and 6–9 days when maintained at 20°C (Vinogradova, 2000). In contrast, after diapause termination, the sister taxon Cx. pipiens f. pipiens generally requires 10 days before ovarian follicle development progresses and egg laying commences (Tate and Vincent, 1936; Sanburg and Larsen, 1973).

To count the number of eggs laid, we used a Nikon stereoscopic microscope (model C-PS) with a Gosky 10× microscope Smartphone Camera Adaptor to photograph each egg raft at 40× magnification. The eggs within each image were then marked and counted digitally using the program ImageJ v. 1.53a (Schneider et al., 2012). Because a female's size can influence the number of eggs she lays (Vinogradova, 2000), we also measured wing length as a proxy for size using a metric miniscale (https://www.bioquip.com/search/DispProduct.asp?pid=4828E). Females were killed with ethyl acetate prior to wing measurement.

Because of the outcome of this experiment (see Results), we wanted to examine potential temperature differences within our incubators between ‘lights on’ and ‘lights off’ conditions. To do this, we placed two sets of 40, 100 ml glass jars filled with 50 ml of water into two incubators (one set of 40 jars per incubator), to simulate the conditions of our dormancy induction experiment. One incubator was set for a 12 h:12 h light:dark cycle (lights on from 02:00 h to 14:00 h) and the other incubator was set for the reverse cycle of 12 h:12 h dark:light (lights on from 14:00 h to 02:00 h). This offset was done to partially account for potential bias that could occur over the course of a day from external factors. Both incubators were set to a constant 18°C. On two non-consecutive days, we measured the internal temperature in each incubator 5 times when the lights were on and 5 times when the lights were off (10 measurements in total, per condition, per incubator). We allowed at least 1 h to elapse between measurements and we made no measurements for at least 1 h after the transition from light to dark (or vice versa). The temperature measurements were made using an OMEGA™ 2 Channel Multilogger digital thermometer with two Omega SC-GG-K-30-36-PP Insulated Type K Thermocouples secured inside each of the incubators (one each). This allowed us to make our measurements without opening the incubator doors. To determine whether there was a statistical difference between ‘lights on’ and ‘lights off’ for either incubator, we compared the measurements for each condition using an unpaired t-test as implemented in R. v.4.0.2 (http://www.R-project.org/), with statistical significance considered at P<0.05.

To complement the experimental examinations of our New York City-derived population of Cx. pipiens f. molestus, we also wanted to identify Cx. pipiens f. molestus-specific genetic variation in genes potentially influencing circadian rhythms and diapause. To do this, we took advantage of the assembled and annotated genome of the closely related mosquito Culex quinquefasciatus (VectorBase assembly CpipJ2; Arensburger et al., 2010). This taxon, known as the southern house mosquito, is a member of the Culex pipiens species complex (as are Cx. pipiens f. molestus and Cx. pipiens f. pipiens), and its annotated genome is frequently used in genetic and genomic comparisons between Cx. pipiens f. molestus and Cx. pipiens f. pipiens (e.g. Asgharian et al., 2015; Price and Fonseca, 2015; Aardema et al., 2020; Yurchenko et al., 2020).

We first needed to determine which annotated Cx. quinquefasciatus genes potentially influence circadian rhythms and/or diapause. To do this, we used the search words ‘circadian’, ‘photoperiod’, ‘dormancy’, ‘light stimulus’ and ‘diapause’ to locate GO categories associated with these terms in Drosophila melanogaster (FlyBase, Dmel Release 6.32; Thurmond et al., 2019). The specific GO terms found are listed in Table S1. Next, we determined all D. melanogaster genes that had a ‘biological process’ associated with one or more of these terms. The D. melanogaster peptide sequences for these genes were compiled into a dedicated FASTA file. Using a local protein–protein BLAST v.2.7.1 (‘blastp’; Altschul et al., 1990) program and these D. melanogaster sequences as our queries, we searched the annotated peptide sequences of Cx. quinquefasciatus, with an e-value of 1e−10 and default values for all other settings. We considered the ‘best’ matches (the longest sequence with the highest percentage similarity) between D. melanogaster and Cx. quinquefasciatus as prospective orthologs, provided they had greater than 50% amino acid similarity across 100 or more amino acids. This list of potential circadian genes in Culex contained 154 unique annotated sequences (Table S2). We also included the four genes related to ‘diapause’ that were recently found to be upregulated in Cx. pipiens f. pipiens but not Cx. pipiens f. molestus (Kang et al., 2020). Where noted, it is these 158 genes that we utilized in subsequent analyses (see below). Henceforth, these genes will be referred to as the ‘putative Culex circadian genes’.

We wanted to determine whether one or more of our putative Culex circadian genes harbored major Cx. pipiens f. molestus-specific structural variants or other potentially ‘inactivating’ mutations that could correlate with its insensitivity to photoperiod in the induction of a dormancy state. To do this, we mapped previously published Cx. pipiens f. molestus genomic short-read Illumina data (NCBI SRA accession numbers: SRR10053379, SRR10053380 and SRR10053386) to the Cx. quinquefasciatus genome (Arensburger et al., 2010) using the ‘MEM’ algorithm implemented in the program BWA v.0.7.15 (Li, 2013 preprint) with default settings. The sample used to produce these genomic reads came from a New York City-derived lab strain of Cx. pipiens f. molestus (Aardema et al., 2020). After mapping, we marked duplicate reads with the MarkDuplicates function in Picard v.1.77 (http://broadinstitute.github.io/picard/), then performed indel realignment using the IndelRealigner function of the Genome Analysis Toolkit v.3.8 (‘GATK’; McKenna et al., 2010). Next, we determined genotype likelihoods using the ‘mpileup’ command in bcftools v.1.9 (Li et al., 2009), then identified sites divergent from the reference using the ‘call’ command, also with bcftools. These divergent sites included both SNPs and insertion/deletions (INDELs). We then used the program SnpEff v.4.3 (Cingolani et al., 2012) with a custom database for the gene annotations of Cx. quinquefasciatus to annotate these variants with default parameters. We cross-referenced the gene summary produced by the SnpEff program with our list of 158 putative Culex circadian genes, identifying those on the list that were determined to have a ‘high’ impact variant. Potential high-impact variants included frameshift mutations, loss of the start codon, gain of a stop codon, loss of the stop codon or disruption of an exonic splice site. This cross-referencing generated 33 putative Culex circadian genes that had one or more potential structural variants or other mutations of great effect.

To test whether these variants were valid (i.e. not an artifact of incorrect annotation), we needed to confirm that the expressed coding gene structure in Cx. pipiens f. molestus matched that of the Cx. quinquefasciatus gene annotation. We also wanted to determine that the variants were found in a broad geographic representation of Cx. pipiens f. molestus samples, and absent in the sister taxon Cx. pipiens f. pipiens. To do this, we de novo assembled two Cx. pipiens f. molestus transcriptomes (from samples deriving from the USA and Germany, respectively; Table S3), and two Cx. pipiens f. pipiens transcriptomes (also from samples deriving from the USA and Germany), using the program Trinity v.2.8.4 with default parameters (Haas et al., 2013). We performed a nucleotide–nucleotide BLAST v.2.7.1 (‘blastn’) search with the target exon containing the potential inactivating mutations as our query sequence and each of these four assembled transcriptomes as our database. Using RNA-seq (transcriptome) data allowed us to simultaneously confirm accurate annotation of the gene and that the focal variant was present in both North American and European Cx. pipiens f. molestus, but not in Cx. pipiens f. pipiens.

If our analysis of these four transcriptomes supported a correct annotation of the focal gene and additionally if the potential ‘large-effect’ variant was present in both Cx. pipiens f. molestus samples and neither Cx. pipiens f. pipiens sample, we then examined its presence/absence in four additional de novo assembled genomes from two Cx. pipiens f. molestus samples (one from Belarus and the aforementioned New York City sample; see Table S3), and two Cx. pipiens f. pipiens samples (one from Belarus and one from NJ, USA; see Table S3). These de novo assemblies were done with the program ABySS v.2.2.3 (Jackman et al., 2017), with K values from 56 to 96 at 10 nucleotide intervals. The ‘best’ assembled genome was chosen as the K value with the highest E-size (a measure of probable gene completeness; Lian et al., 2014), as determined with ‘abyss-fac’. Again, we used a nucleotide–nucleotide BLAST v.2.7.1 (‘blastn’) program with the target exon containing the potential inactivating mutations as our query sequence and each of these four assembled genomes as our database.

In addition to possible structural variants in genes known to influence diapause and circadian rhythms, it is possible that derived amino acid changes could also impact the expression of these traits. To examine this possibility, we took advantage of previous research comparing non-synonymous divergence (K_a) between Cx. pipiens f. molestus from New York City and Cx. pipiens f. pipiens from Baden-Württemberg, Germany (Price and Fonseca, 2015). We first compared our list of putative Culex circadian genes with those which were found to have a K_a value greater than zero (indicating non-synonymous divergence between Cx. pipiens f. molestus and Cx. pipiens f. pipiens). From this comparison, we compiled a list of 49 gene candidates to investigate potential derived amino acid changes in Cx. pipiens f. molestus.

Next, we used a local nucleotide–nucleotide BLAST v.2.7.1 (‘blastn’) program to compare each exon within each of these 49 genes from the annotated Cx. quinquefasciatus genome to each of the four previously described de novo genome assemblies (two Cx. pipiens f. molestus samples and two Cx. pipiens f. pipiens samples; see above). We aligned these four focal sequences with the Cx. quinquefasciatus focal exon and the full gene sequence to maintain the correct reading frame (in relation to the annotation). Sequence gaps, codons that spanned an intron and regions that were not present in all four focal genomes were removed. After all exon regions were located, aligned and trimmed, these sequences were concatenated into a single sequence, then converted to amino acids. Alignment, trimming and amino acid conversion were done with SeaView v.4.6.3 (Gouy et al., 2010). Using a custom Perl script (available from M.L.A. on request), we counted the number of derived amino acid changes observed in each gene (relative to Cx. quinquefasciatus).

To further assess the presence of potential inactivating mutations in putative Culex circadian genes, we generated RNA-seq data for four Cx. pipiens f. molestus life stages: larva, pre-pupa, pupa and adult. For each library, we combined 10 individual heads, then extracted mRNA using Trizol. The goal of this analysis was not to quantify differences in expression levels per se, but rather to look for the complete absence of any gene expression. Because of resource limitations, it was most cost effective to increase the number of reads sequenced for each of the four life-stage libraries, rather than produce independent biological replicates for each life stage. Pooled samples, especially when the number of reads is high, have been shown to improve the power to detect gene expression of transcripts present at low abundance (Takele Assefa et al., 2020). While not the optimal experimental design, our use of pooled samples here, sequenced to produce a relatively high number of reads (>20 million each), could suggest something about binary gene expression patterns (i.e. a gene is expressed or not expressed), even for genes only expressed here at low levels.

After RNA extraction, we next prepared sequencing libraries following the Illumina TruSeq protocol, incorporating Covaris shearing as an alternative to chemical shearing and excluding the cDNA DSN normalization step. Paired-end sequencing was performed on a HiSeq2000 platform. We mapped the resulting reads from each sample to the Cx. quinquefasciatus genome with STAR v.2.5.2 (Dobin et al., 2013; Dobin and Gingeras, 2015) as implemented in RSEM v.1.3.1 (Li and Dewey, 2011). RSEM was then used to calculate the transcripts per kilobase million (TPM) for all genes in each of the four samples. We considered any measure of TPM above zero as evidence of possible gene expression.

The Rayleigh test indicated that all treatments displayed clustering in the data (Fig. 1; 12 h:12 h light:dark: 0.6042, P<0.001; dark:light: 0.3413, P<0.001, 24 h light: 0.2051, P=0.009; 24 h dark: 0.2301, P=0.0045). The mean (±s.d.) emergence time for Cx. pipiens f. molestus adult eclosion for individuals reared in 12 h:12 h light:dark (lights on at 06:00 h, lights off at 18:00 h) was 21:58 h (±1 h), approximately 4 h after the onset of the dark cycle during larval development. For individuals reared in 12 h:12 h dark:light, the mean was 10:08 h (±1 h 28 min), again approximately 4 h after the onset of the dark cycle during larval development. The mean adult eclosion time for individuals reared in constant light (24 h light) was 20:15 h (±1 h 47 min), and for individuals reared in constant dark (24 h dark) it was 17:50 h (±1 h 43 min). In our pairwise comparisons, we observed a statistically significant difference between all treatments except between constant light and constant dark (Table 1).

Of the 109 females set up for oviposition from those reared in long-day conditions, 102 laid eggs (93.6%), and of the 101 females set up for oviposition from those reared in short-day conditions, 94 laid eggs (93.1%). On average, a female reared in long-day conditions took 38.4±33.8 h (mean±s.d.) to lay eggs, and females reared in short-day conditions took 47.9±45.1 h. The variance between the two groups was not statistically different (Levene's test, F_1,192=1.79, P=0.1825), so we performed an unpaired t-test assuming equal variances. The results of this test indicated that there was no statistically significant difference in the time it took for females to lay eggs between treatments (t₁₉₂=−1.6683, P=0.097).

On average, females reared in long-day conditions laid 81.4±18.3 eggs (mean±s.d.), whereas females reared in short-day conditions laid an average of 97.8±15.8 eggs. There was a statistically significant difference between the two treatments (t_181.32=−6.5394, P<0.001). An ANOVA additionally showed that these differences between treatments were significant (F_1,192=42.21, P<0.001). However, there is a strong intraspecific relationship between mosquito size and number of eggs laid (Vinogradova, 2000), and we observed that on average female mosquitoes reared in long-day conditions had a wing length of 4.38±0.21 mm (mean±s.d.), whereas females reared in short-day conditions had an average wing length of 4.68±0.16 mm. This difference was significant (t_174.67=−10.851, P<0.001). When we controlled for these observed differences in size using an ANCOVA, there was no difference in the number of eggs laid between the treatments (F_1,181=0.350, P=0.555). Fig. 2 shows the relationship between wing size and the number of eggs laid for both treatment groups.

Our post-rearing analysis of temperature variation indicated that there was a small but statistically significant difference within each of the two examined incubators between ‘lights on’ and ‘lights off’ conditions, with lights on being 0.19°C warmer on average (data combined across both incubators). See Table S4 for specific values and statistical comparison results for each incubator.

Our analysis of potential ‘large effect’ variants in putative Culex circadian genes identified 33 candidate genes (Table S5). However, upon further examination, for 28 of these genes the variant was not observed in all Cx. pipiens f. molestus samples surveyed. For three other genes, the variant was also present in one or more Cx. pipiens f. pipiens samples. Lastly, for one gene, the variant was mis-characterized as a result of an apparent inaccuracy in the Cx. quinquefasciatus genome annotation. The single variant that we determined to be correctly annotated and was present in all examined Cx. pipiens f. molestus samples and absent in all examined Cx. pipiens f. pipiens samples, was a nine nucleotide, in-frame deletion in the fifth exon of the Helicase domino (dom) gene (Fig. 3). We observed this variant in the examined Cx. pipiens f. molestus RNA-seq and genomic data from New York City, Germany and Belarus, but not in the examined Cx. pipiens f. pipiens samples from similar, geographically proximate locations, or in the Cx. quinquefasciatus reference genome.

Based on previous analysis of non-synonymous divergence between Cx. pipiens f. molestus and Cx. pipiens f. pipiens (Price and Fonseca, 2015), we examined 49 putative Culex circadian genes for derived amino acid changes relative to Cx. quinquefasciatus. Of these, eight genes had one derived amino acid change (missense mutation) in both the examined New York City and Belarussian Cx. pipiens f. molestus genome samples, but this was absent in the New Jersey and Belarussian Cx. pipiens f. pipiens genome samples, and in the Cx. quinquefasciatus reference (Table 2; Table S6). As annotated in Cx. quinquefasciatus, these genes were: calmodulin binding transcription activator 2, sodium chloride dependent amino acid transporter, dna photolyase, calmodulin-binding protein trpl, ultraviolet-sensitive opsin, glycogen synthase kinase 3, phospholipase c and a conserved hypothetical protein. There were 12 genes that had derived amino acid changes in both examined Cx. pipiens f. pipiens samples, with many genes harboring more than one missense mutation (21 total derived amino acid changes in Cx. pipiens f. pipiens).

Our sequencing of four pooled libraries each constituting a distinct Cx. pipiens f. molestus life stage (larva, pre-pupa, pupa and adult) resulted in over 21 million read pairs per library (larvae: 24.7 million, pre-pupa: 23.6 million, pupa: 21.1 million, adult: 22.4 million). The sequencing data are deposited in the NCBI SRA database under accession numbers SRR13321044–SRR13321047 (BioProject PRJNA688455). Of the 158 identified putative Culex circadian genes, all but four had evidence of expression in at least one Cx. pipiens f. molestus life-history stage (Table 3; Table S6). These four genes are annotated in the Cx. quinquefasciatus genome as AMP dependent ligase, tubulin beta-3 chain, Dual specificity tyrosine-phosphorylation-regulated kinase and a gene encoding an uncharacterized protein. This last gene appears most similar to E3 ubiquitin-protein ligase TRIP12 in D. melanogaster. Six other genes were expressed at very low levels (TPM<1) in only one of the four examined life stages (Table 3; Table S6). Considering the four genes related to ‘diapause’ recently found to be upregulated in Cx. pipiens f. pipiens, but not Cx. pipiens f. molestus (Kang et al., 2020), all showed some level of expression in at least one of the four examined life stages. However, the maximum expression level for one of these genes, a putative juvenile hormone esterase precursor (CPIJ019485), was just 2.78 transcripts per million in the larval data, with lower expression levels in the other three life stages (Table 3).

In contrast to our predictions, our examination of a circadian behavior (time of adult eclosion) in a New York City-derived line of Cx. pipiens f. molestus indicated that it retains the ability to respond to variation in exogenous photoperiod. This is despite it originating from, and being maintained in, a light asynchronous environment. However, our experimental assessment of ecologically relevant, reproductive dormancy in this mosquito showed that in conditions reported to induce diapause in the sister taxon Cx. pipiens f. pipiens, there were no clear reductions in the tendency to lay eggs, time to lay or number of eggs laid.

Although these results strongly suggest that this line of mosquito lacks any tendency towards ecologically relevant dormancy in response to a shortened photoperiod, we note here that as part of our experimental design, females reared in short-day conditions during their larval and pupal development were subsequently moved to long-day conditions for the duration of their adult life stage. Previous research with Cx. pipiens f. pipiens showed that a sudden shift from a short-day to a long-day photoperiod in adulthood can terminate diapause initiation (Sanburg and Larsen, 1973). This suggests the possibility that our short-day-reared mosquitoes may have been physiologically affected by their rearing conditions (i.e. induced towards some level of dormancy), but the transition to a long-day photoperiod after adult eclosion prevented us from observing this. However, in their experiment, Sanburg and Larsen (1973) observed that after long-day photoperiod exposure, it took approximately 10 days for the ovarian follicles of previously diapause-induced females to reach sizes comparable to those in females reared in long-day conditions throughout their lives. If our short-day-reared females had been induced towards a dormant state, but this was subsequently disrupted by the transition to a long-day photoperiod in adulthood, then we predict evidence of this would manifest itself in the length of time between eclosion and oviposition in these mosquitoes. Although our short-day-reared females did take an average of 9.5 h longer to lay eggs relative to long-day-reared females (47.9 versus 38.4 h), the differences overall were not statistically significant. Furthermore, this time plus the 72 h post-eclosion that these females were housed with males in flight cages means that on average they laid eggs 5 days after eclosion. As noted above, this is half the time previously observed to be necessary for full development of the ovarian follicles in females moved from short-day to long-day conditions (Sanburg and Larsen, 1973).

Interestingly, female Cx. pipiens f. molestus reared in short-day conditions (6 h:18 h light:dark) during the larval stage actually laid more eggs than females reared in long-day conditions (18 h:6 h light:dark). This is the opposite trend we would predict if females reared in short-day conditions were induced towards any degree of reduced reproductive output by the short-day rearing conditions. The proximate explanation for the greater average number of eggs laid by females reared in short-day conditions is that these females were generally larger. There is a strong positive correlation between Culex female size and the number of eggs laid, and lower temperatures during the larval period result in larger adult mosquitoes (Vinogradova, 2000). Our post-experiment assessment of incubator temperature variation between ‘lights on’ and ‘lights off’ conditions indicated that differences, although minimal, were significant. Incubators were on average 0.19°C warmer when the lights were on (data consolidated across incubators). This small but consistent temperature difference may account for the size differences observed (and corresponding variation in fecundity). Although much is still unknown about the mechanistic relationships between temperature, growth rates, size and fecundity in mosquitoes, there are several factors that could play a role in the results we observed. Foremost, higher temperatures speed up larval development but can also increase mortality, potentially as a result of increased stress (Mpho et al., 2002; Ciota et al., 2014). Genetically, temperature influences the expression of many important regulatory genes, the most well known being the heat shock proteins (reviewed in King and MacRae, 2015; Chen et al., 2018). Such gene expression variation may have important consequences for larval growth and fecundity. Temperature variation during mosquito larval development can also influence the gut microbial community, which in turn may influence how efficiently food is converted to biomass and future resources for egg development (Onyango et al., 2020). Although slight but consistent temperature variation is a clear hypothesis to explain the observed differences in size and fecundity between mosquitoes reared in short-day versus long-day conditions, it cannot be discounted that variations in the light environment itself could have influenced growth and development through currently unknown mechanisms. This is an intriguing possibility that will require future assessment.

In our experiment described here, the complete absence of any apparent reduction of reproductive output in response to short-day, diapause-inducing conditions strongly suggests that the examined line of Cx. pipiens f. molestus lacks any ability to exhibit a photoperiodically inducible dormant state. However, the results from our circadian rhythm study indicate that Cx. pipiens f. molestus circadian behavior can be exogenously entrained by light cues (specifically the timing of the photoperiod) during the larval stage. This refutes our hypothesis that genetic changes in one or more Cx. pipiens f. molestus clock gene could simultaneously account for its inability to enter diapause and produce either a loss of entrainable circadian rhythms or a diminished sensitivity to light cues. Our results show that this taxon maintains photoperiodic perception and some degree of circadian entrainment influenced by this perception. However, as we only investigated photoperiodic entrainment within a single Cx. pipiens f. molestus population, it remains unclear whether the degree and strength of this entrainment differ from that in other populations or taxa in the Culex genus.

Despite their ecological differences in habitat preference and diapause induction ability, there was minimal divergence in genes potentially influencing circadian rhythms and/or overwintering behavior among the Cx. pipiens samples examined in this study. This agrees with the close, and often challenging, taxonomic relationship of these mosquitoes (e.g. Smith and Fonseca, 2004; Fonseca et al., 2004; Aardema et al., 2020). Of the 158 putative Culex circadian genes we examined, only one, Helicase domino (dom), harbored a derived, Cx. pipiens f. molestus-specific structural variant, specifically a nine-nucleotide, in-frame deletion. This variant was not observed in either the Cx. pipiens f. pipiens samples or the Cx. quinquefasciatus reference genome. The conserved ATPase domains present in the D. melanogaster ortholog of this gene indicate it is part of the SWI/SNF2 DNA-dependent ATPase family (Ruhf et al., 2001). Genes in this family are important for chromatin remodeling (reviewed in Tang et al., 2010), and in D. melanogaster, the DOM protein has many important functional associations, including wing and dendrite development, cell viability and proliferation, neuroblast maintenance and polarity, and oogenesis (Ruhf et al., 2001; Tea and Luo, 2011; Börner and Becker, 2016; Rust et al., 2018; Liu et al., 2019). Based on our alignment of the Cx. quinquefasciatus sequence with the presumed dom ortholog in D. melanogaster, the variant observed here appears between ATPase domains four (IV) and five (V), and upstream of a PEST sequence (Ruhf et al., 2001). The presence of PEST sequences in this gene suggests that the corresponding protein product is rapidly degraded after translation (Rogers et al., 1986), further indicating that dom is an important regulator of physiological processes. Given its likely multifunctionality in Culex, it is perhaps not surprising that the observed structural variant in the dom gene of Cx. pipiens f. molestus would not radically alter its sequence. However, it is still possible that the observed variant does have an influence on the expression of circadian rhythms and, correspondingly, diapause induction. In D. melanogaster, distinct splice variants of the dom gene dramatically impact circadian behaviors (Liu et al., 2019). Both dom isoforms bind to the promoters of the negative circadian regulators PER and TIM, and act to control their abundance in neural cells. Downregulation of the first dom isoform (encoding DOM-A) via RNAi leads to arrhythmic behavior, whereas downregulation of the second isoform (encoding DOM-B) lengthens the circadian period, while retaining its rhythmicity (Liu et al., 2019). In Cx. pipiens f. pipiens females reared in diapause-inducing conditions, knocking down the expression of both per and tim via RNAi caused these mosquitoes to avert diapause, with corresponding ovarian development and failure to increase lipid reserves (Meuti et al., 2015). The effect that the nine-nucleotide deletion we observed in the Cx. pipiens f. molestus dom gene has on the functioning of the DOM protein, particularly in relation to per and tim expression, will require further investigation.

In addition to the structural variant of dom, we also found eight putative Culex circadian genes that each harbored one derived, non-synonymous amino acid change (missense mutation) in Cx. pipiens f. molestus samples but not in Cx. pipiens f. pipiens samples. One of these genes encodes for the protein PHOSPHOLIPASE C. In D. melanogaster, this gene (norpA, FBgn0262738) is predominantly expressed in the eyes, and mutations in this gene ultimately affect the visual input pathway and circadian entrainment (Collins et al., 2004). This gene may also regulate splicing of the per gene, which could impact both circadian rhythms and diapause. Another gene observed to have a non-synonymous change in Cx. pipiens f. molestus was glycogen synthase kinase 3 (GSK3). In D. melanogaster, variation in expression levels of this gene (sgg, FBgn0003371) correspondingly influence patterns of phosphorylation of TIM and effect nuclear translocation of the heterodimer formed by TIM and PER (Martinek et al., 2001). Intriguingly, expression of another glycogen synthase gene in Cx. pipiens f. pipiens was shown to impact the regulation of glycogen and lipid storage during diapause, and it was deemed essential for survival during winter dormancy (King et al., 2020). While the observed amino acid change in GSK3 could potentially impact the ability of Cx. pipiens f. molestus to enter a dormant state, it may also correlate with the absence of diapause in this mosquito and the lack of a need to accumulate large lipid reserves.

Among the genes for which no expression was detected in any assessed Cx. pipiens f. molestus life stage, AMP dependent ligase is perhaps most interesting. This gene appears to be orthologous to the Drosophila gene Very long-chain-fatty-acid-CoA ligase bubblegum (bgm, FBgn0027348), which is predominantly a metabolic gene that influences fatty acid and lipid metabolism, but which also regulates the circadian sleep/wake cycle (Thimgan et al., 2015). Furthermore, in bumble bees, this gene appears to be upregulated prior to the onset of diapause (when metabolic reserves are being accumulated in the body), and downregulated during the actual diapause period (Amsalem et al., 2015). More broadly, it has been found that many metabolic genes are differentially expressed between diapausing and non-diapausing Cx. pipiens f. pipiens females (Kang et al., 2016), indicating the potential importance of such genes for successful diapause. Yet more evidence of this comes from a follow-up study in which numerous metabolic genes were also found to be upregulated in Cx. pipiens f. pipiens adult females relative to Cx. pipiens f. molestus adult females (Kang et al., 2020). Although our analysis of RNAseq data revealed these intriguing gene candidates pertaining to a possible absence of expression in Cx. pipiens f. molestus, these results must be considered tentative given the limitations of our experimental design.

Fundamentally, this study has only revealed potential candidates with genetic changes specific to Cx. pipiens f. molestus that may correlate with its inability to enter a dormant state in response to photoperiodic cues. We can provide no causative evidence that these changes contribute to this phenotype. More dedicated work with expanded biological replicates will be required to confidently identify genetic variation that may contribute to an absence of photoperiodically induced dormancy in this taxon. Furthermore, the changes we have characterized are all likely to be derived within the Culex pipiens species complex. The justification for focusing on such derived genetic variation is that Cx. pipiens f. molestus is generally presumed to have evolved relatively recently (∼10,000–80,000 years; Fonseca et al., 2004; Shaikevich, 2007), possibly from a Cx. pipiens f. pipiens ancestor (Byrne and Nichols, 1999; Kothera et al., 2010). However, the evidence that Cx. pipiens f. molestus is the derived taxon is limited, and it is interesting to note that our outgroup for comparisons, Cx. quinquefasciatus, also lacks the ability to enter a dormant state (Wilton and Smith, 1985). If an absence of diapause is the ancestral condition in this mosquito group, and if Cx. pipiens f. molestus did not recently derive from facultatively diapausing Cx. pipiens f. pipiens, then it is unlikely that derived genetic variation in contemporary populations is responsible for its inability to enter a dormant state in response to shortened photoperiods. If future work suggests this alternative hypothesis is more probable, then it may be interesting to examine the evolution of diapause as a derived trait in Cx. pipiens f. pipiens.

We have shown that in the urban-adapted underground mosquito Cx. pipiens f. molestus, the ability to entrain circadian rhythms in response to exogenous light cues is decoupled from its inability to enter a photoperiodically induced dormancy state. Because it can mate in enclosed spaces and does not require a blood meal to reproduce for the first time, Cx. pipiens f. molestus is commonly maintained in entomology and disease vector laboratories worldwide. Given the wide usage of this taxon in understanding mosquito biology, our results indicate that it may be a valuable tool for exploring exogenously influenced phenotypes, particularly those that display a circadian rhythm. Greater knowledge of circadian rhythms in mosquitoes is critical for controlling and mitigating mosquito-vectored illnesses (Rund et al., 2016). Specifically, understanding how seasonal changes in photoperiod fine-tune daily mosquito behaviors could greatly improve our knowledge of transmission potential for specific pathogens and vectors. In this study, we uncovered a structural variant in the Helicase domino gene that segregates with the taxa examined here. We postulate that this variant is a contributing factor to the inability of Cx. pipiens f. molestus to enter diapause. Given the substantial influence this gene appears to have on circadian rhythms in Drosophila, and given that circadian genes are known to greatly impact diapause induction in Culex, this gene represents a major target for follow-up studies. Additional genetic variation uncovered here specific to Cx. pipiens f. molestus also offers further opportunities to investigate the genetic underpinnings of the diapause trait. A better understanding of the genetic variation influencing diapause and circadian rhythms in Culex pipiens species complex mosquitoes more broadly may lead to improved vector control and a reduction in disease transmission.

We would like to thank David Epstein for helping to design and build the lighting mechanism used in our investigation of circadian rhythms. Prof. Bridgett vonHoldt graciously provided space for our diapause induction experiment. D. Fonseca and anonymous reviewers provided helpful comments on earlier drafts of this work.