Mosquito rearing and maintenance. Ae. aegypti mosquitoes used in all experiments were derived from the Ae. aegypti Liverpool strain, which was the source strain for the reference genome sequence 43. Mosquito rearing and maintenance were performed following previously established procedures 17.
Guide RNA design and assessment. To induce ♀-specific lethality, doublesex (dsx, AAEL009114) and Intersex (ix, AAEL010217) genes were targeted for CRISPR/Cas9-mediated disruption. For each target gene, the DNA sequences were initially identified using reference genome assembly, and then genomic target sites were validated using PCR amplification and Sanger sequencing (see Supplementary Table 1 for primer sequences). CHOPCHOP V3.0.0 (https://chopchop.cbu.uib.no) was used to select the gRNA target sites. In total, we selected 4 dsx gRNAs and 2 ix gRNAs, and assembled the gRNAdsx and gRNAix constructs harboring the corresponding gRNAs. To confirm the activity of the gRNAs in vivo, we generated the transgenic strains harboring either the gRNAdsx or gRNAix constructs, and crossed each gRNA strain to the Cas9 strain17. We sequenced DNA regions targeted by each selected gRNA in F1 trans-hemizygous (aka. gRNA/+; Cas9/+) mosquitoes to assess gRNA/Cas9-induced mutagenesis.
Genetic construct design and plasmid assembly. All plasmids were cloned using Gibson enzymatic assembly44. DNA fragments were either amplified from available plasmids or synthesized gBlocks™ (Integrated DNA Technologies) using Q5 Hotstart Start High-Fidelity 2X Master Mix (New England Biolabs, Cat. #M0494S). The generated plasmids were transformed into Zymo JM109 chemically competent E. coli (Zymo Research, Cat # T3005), amplified, isolated (Zymo Research, Zyppy plasmid miniprep kit, Cat. #D4036), and Sanger sequenced. Final selected plasmids were maxi-prepped (Zymo Research, ZymoPURE II Plasmid Maxiprep kit, Cat. #D4202) and sequenced thoroughly using Oxford Nanopore Sequencing at Primordium Labs (https://www.primordiumlabs.com). All plasmids and annotated DNA sequence maps are available at www.addgene.com under accession numbers: 200252 (gRNAdsx, 1067B) 200251 (gRNAix, 1055J), 200253 (gRNAdsx,ix,βTub, 1067L) and 164846 (Nup50-Cas9 or Cas9, 874PA).
Engineering transgenic strains. Transgenic strains were generated by microinjecting preblastoderm stage embryos (0.5–1 hr old) with a mixture of the piggyBac plasmid (100 ng/ul) and a transposase helper plasmid (pHsp-piggyBac, 100 ng/ul). Embryonic collection, microinjections, transgenic lines generation and rearing were performed following previously established procedures 17,45,46.
Genetic assessment of gRNA strains. To evaluate the activity of gRNAdsx and gRNAix strains, we reciprocally crossed hemizygous mosquitoes of each generated strains to the homozygous Nup-Cas9 strain17, and thoroughly examined F1 trans-hemizygous ♀ progeny (Supplementary Tables 2–3). Each of three established gRNAdsx,ix,βTub strains was subjected to single-pair sibling matings over 5–7 generations to generate the homozygous stocks. Zygosity was confirmed genetically by repeated test crosses to WT. To generate trans-hemizygous progeny with a maternal Cas9 (pgSIT♀Cas9) and paternal Cas9 (pgSIT♂Cas9), homozygous gRNAdsx,ix,βTub ♂’s and ♀’s were mated to homozygous Cas9 ♀’s and ♂’s, respectively. Then, the generated trans-hemizygous mosquitoes were mated to WT mosquitoes of the opposite sex, and numbers of laid and hatched eggs were scored to measure the fecundity. Multiple genetic control crosses were also performed for comparisons: WT ♂ 𝗑 WT ♀; WT ♂ 𝗑 Cas9 ♀; Cas9 ♂ 𝗑 WT ♀; gRNA ♂ 𝗑 Cas9 ♀; gRNA ♀ 𝗑 Cas9 ♂; gRNA ♀ 𝗑 WT ♂; and gRNA ♂ 𝗑 WT ♀. For the genetic cross, adult mosquitoes were allowed to mate in a cage for 4–5 days, then blood meals were provided, and eggs were collected and hatched. To calculate the normalized percentage of laid eggs per ♀ (e.i. fecundity), the number of laid eggs per ♀ for each genotype replicate was divided by the maximum number of laid eggs per ♀ (Supplementary Table 4). The percentage of egg hatching (i.e. fertility) was estimated by dividing the number of laid eggs by the number of hatched eggs. Egg-to-adult survival ratios were the number of laid eggs divided by the total number of eclosed adults. Larvae-to-adult survival rates were calculated by dividing the number of eclosed adults by the number of larvae. Pupae-adult survival rates were calculated by dividing the number of pupae by the number of eclosed adults. Blood feeding rates were calculated by dividing the number of blood-fed ♀’s or ⚥’s by the total number of ♀ or ⚥’s. Trans-hemizygous ⚥’s had many morphological features visibly different from those of WT ♀’s. To assess differences in external and internal anatomical structures, we dissected twenty trans-hemizygous and control mosquitoes for each genotype and sex (♂, ♀ and ⚥) in 1% PBS buffer and imaged structures on the Leica M165FC fluorescent stereomicroscope equipped with a Leica DMC4500 color camera. The measurements were performed using Leica Application Suite X (LAS X by Leica Microsystems) on the acquired images (Supplementary Table 2).
Validation of induced mutations. To confirm that the ♀ lethality or transformation correlated with specific mutations at the target loci, dsx and ix gRNA targets were PCR amplified from the genomic DNA extracted from trans-hemizygous and control mosquitoes. PCR amplicons were purified directly or gel purified using Zymoclean Gel DNA Recovery Kit (Zymo Research Cat. #D4007) and then Sanger sequenced. Presence of induced mutations were inferred as precipitous drops in sequence read quality by analyzing base peak chromatograms (Supplementary Fig. 2). The primers used for PCR and sequencing, including gRNA target sequences, are listed in Supplementary Table 1. In addition, induced mutations and their precise localization to six target sites in ix, dsx, and βTub genes were confirmed by sequencing genomes of gRNAdsx,ix,βTub#1/Cas9 ♂’s and ⚥’s using Oxford Nanopore (see below) and aligning them to the AaegL5 genome (GCF_002204515.2). Mutations in the DNA using Nanopore data and mutations in the RNA using RNAseq data were visualized using an integrated genome browser (Supplementary Figs. 7–10).
Ear Johnston’s Organ microanatomy and Immunohistochemistry (IHC). The heads of two-day old mosquitoes were removed, fixed in 4% PFA and 0.25% Triton X-100 PBS, and sectioned using a Leica VT1200S vibratome. The anti-SYNORF1 3C11 primary monoclonal antibody (AB_528479, 1:30, Developmental Studies Hybridoma Bank (DSHB), University of Iowa) and two secondary antibodies, Alexa Fluor 488-conjugated anti-mouse IgG (#A-11029, 1:300, ThermoFisher) and anti-horseradish peroxidase 555 (anti-HRP, AB_2338959, 1:100, Jackson ImmunoResearch), were used for IHC. All imaging was conducted using a laser-scanning confocal microscope (FV3000, Olympus). Images were taken with a silicone-oil immersion 60× Plan-Apochromat objective lens (UPlanSApo, NA = 1.3) as previously described 47.
Wing Beat Frequency (WBF) measurements. Groups of 30 mosquitoes were transferred to 15cm3 cages containing a microphone array (Knowles NR-23158-000) and provided with a source of 10% glucose water. After two days of LD entrainment, audio recordings of mosquito flight behavior were made using a Picoscope recording device (Picoscope 2408B, Pico Technology) at a sampling rate of 50 kHz. These recordings lasted 30 minutes during periods of peak activity during dusk (ZT12.5-13). Audio recordings were then processed and analyzed using the simbaR package in R 48,49. Four separate cages of mosquitoes were tested for each genotype. Statistical comparisons were made using ANOVA on Ranks tests, followed by pairwise Wilcoxon tests with a Bonferroni correction applied for multiple comparisons.
Phonotaxis assays. Individual ♂’s were provided with a 475 Hz pure tone stimulus calibrated to 80 dB SPL (V ≈ 5 x 10 − 4 ms− 1) provided for 30 seconds by a speaker. Prior to tone playback, a gentle puff of air was provided to each ♂ to stimulate flight initiation. Mosquitoes were considered to demonstrate phonotaxis behaviors if they orientated to the speaker and then demonstrated abdominal bending behaviors. A video camera (GoPro) was used to record all responses to stimulation, with responses being scored as either 0 (no approach/no bending) or 1 (approach and bending). Mosquitoes were tested on two consecutive days, with mosquitoes showing behavior considered as a phonotactic response to stimuli on at least one test day being regarded as responders. Three repeats were conducted for each genotype (across 3 generations), with at least 10 ♂’s tested per repeat. Statistical comparisons were made using Chi-squared tests.
Transgene copy number and genomic integration site. To infer the transgene copy number(s) and insertion site(s), we performed Oxford Nanopore sequencing of the trans-hemizygous gRNAdsx,ixβTub#1/Cas9 mosquitoes. Genomic DNA was extracted from 10 ♂’s and ♀’s using Blood & Cell Culture DNA Midi Kit (Qiagen, Cat# 13343) and the sequencing library was prepared using the Oxford Nanopore SQK-LSK110 genomic DNA by ligation kit. The library was sequenced for 72 hrs on MinION flowcell (R9.4.1) followed by base calling with Guppy v6.3.2 and produced 16.18 Gb of filtered data with Q score cutoff of 10. To determine transgene copy number(s), reads were mapped to the AaegL5 genome supplemented with transgene sequences gRNAdsx,ix,βTub (1067L plasmid) and Cas9 (874PA plasmid) using minimap2 and sequencing depth was calculated using samtools. To identify transgene integration sites, nanopore reads were mapped to the gRNAdsx,ix,βTub and Cas9 plasmids and the sequences that aligned to the expected regions between the piggyBac sites were parsed. These validated reads were then mapped to the AaegL5 genome (GCF_002204515.2) and sites of insertions were identified by examining the alignments in the Interactive Genomics Viewer (IGV) browser. The gRNAdsx,ix,βTub#1 strain harbors a single copy of gRNAdsx,ix,βTub inserted on chromosome 3 ( +/+ orientation) between 2993 and 16524 at NC_035109.1:362,233,930 TTAA site (Supplementary Fig. 5). We also confirmed the previously identified insertion site of Cas917: chromosome 3 (+/+ orientation) between 2933 and 14573 at NC_035109.1:33,210,107 TTAA.
Transcriptional profiling and expression analysis. To quantify target gene reduction and expression from transgenes as well as to assess global expression patterns, we performed Illumina RNA sequencing. We extracted total RNA using miRNeasy Mini Kit (Qiagen, Cat# 217004) from 20 sexed pupae: WT ♀’s, WT ♂’s, trans-hemizygous pgSIT♀Cas9 ♂’s, and pgSIT♀Cas9 ⚥’s in biological triplicate (12 samples total), following the manufacturer’s protocol. DNase treatment was conducted using DNase I, RNase-free (ThermoFisher Scientific, Cat# EN0521), following total RNA extraction. RNA integrity was assessed using the RNA 6000 Pico Kit for Bioanalyzer (Agilent Technologies #5067 − 1513), and mRNA was isolated from ~ 1 µg of total RNA using NEBNext Poly(A) mRNA Magnetic Isolation Module (NEB #E7490). RNA-seq libraries were constructed using the NEBNext Ultra II RNA Library Prep Kit for Illumina (NEB #E7770) following the manufacturer’s protocols. Libraries were quantified using a Qubit dsDNA HS Kit (ThermoFisher Scientific #Q32854), and the size distribution was confirmed using a High Sensitivity DNA Kit for Bioanalyzer (Agilent Technologies #5067 − 4626). Libraries were sequenced on an Illumina NextSeq2000 in paired end mode with the read length of 50 nt and sequencing depth of 20 million reads per library. Base calling and FASTQ conversion were performed with DRAGEN 3.8.4. The reads were mapped to the AaegL5.0 (GCF_002204515.2) genome supplemented with the Cas9 and gRNAdsx,ix,βTub plasmid sequences using STAR aligner. On average, ~ 97.8% of the reads were mapped. Gene expression was then quantified using featureCounts against the annotation release 101 GTF downloaded from NCBI (GCF_002204515.2_AaegL5.0_genomic.gtf). TPM values were calculated from counts produced by featureCounts. PCA and hierarchical clustering were performed in R and plotted using the ggplot2 package. Differential expression analyses between pgSIT vs WT samples within each sex were performed with DESeq2. Gene Ontology overrepresentation analyses were done with the topGO package.
Blood feeding assays. Ten mature WT ♀’s or pgSIT ⚥’s (i.e gRNAdsx+ix+βTub#1/Cas9) generated with Cas9 inherited from either the mother (pgSIT♀Cas9) or father (pgSIT♂Cas9) were caged with 100 mature WT ♂’s for two days to insure that every WT ♀ was mated before WT ♂’s were removed from all cages. Then, one anesthetized mouse (similar size and weight) was put into each cage for 20 minutes (from 9:00 am to 9:20 am) daily for five days, and two times were recorded for each mosquito ♀ or ⚥. First, the ‘time to bite’ is the time interval between a mouse placed in a cage and ♀/⚥ landing on the mouse. Second, the ‘feeding time’ is the duration of blood feeding (Supplementary Table 5). When individual ♀/⚥’s finished blood feeding they were removed from the cage, while unfed ♀/⚥’s were kept in the cage and the blood feeding was repeated the next day. In total, each ♀/⚥ had 6000 seconds to initiate bite and feed during 5 days. To assess whether distance affects blood feeding of pgSIT ⚥’s, we conducted the experiment in two sizes of cages: small (24.5 X 24.5 X 24.5 cm) and large (60 X 60 X 60 cm). All experiments were repeated 3 times independently.
Flight activity assays. We followed the protocol in our previous study17 to assess the flight activity of trans-hemizygous pgSIT♀Cas9 (i.e., gRNAdsx+ix+βTub#1/Cas9) mosquitoes (Supplementary Table 6).
Mating assays. We followed a previously described protocol17, to assess the mating ability of both types of pgSIT males, pgSIT♀Cas9 and Cas9 pgSIT♂Cas9, and confirm that prior matings with pgSIT ♂’s suppressed WT ♀ fertility (Supplementary Table 8).
Fitness parameters and longevity. To assess fitness of pgSIT mosquitoes, we measured developmental times, fertility, and longevity of WT, homozygous gRNAdsx+ix+βTub#1, homozygous Cas9, and pgSIT♀Cas9 and pgSIT♂Cas9 mosquitoes. The previously published protocol was used to assess fitness parameters and longevity17. Briefly, all fitness parameters and longevity were assessed in three groups (3 replicates) of twenty individual mosquitoes (Supplementary Table 7).
Multigenerational population cage trials. To assess the competitiveness of pgSIT ♂’s, we performed discrete non-overlapping multigenerational population cage trials as described in our previous study 17. Briefly, 3/4-day-old mature WT adult ♂’s were released along with mature (3/4 days old) pgSIT adult ♂’s at release ratios (pgSIT: WT): 1:1 (50:50), 5:1 (250:50), 10:1 (500:50), 20:1 (1000:50), and 40:1 (2000:50), with three biological replicates for each release ratio (15 cages total). One hour later, 50 mature (3/4 days old) WT adult ♀’s were released into each cage. All adults were allowed to mate for 2 days. ♀’s were then blood fed and eggs were collected. Eggs were counted and stored for four days to allow full embryonic development. Then, 100 eggs were randomly selected, hatched, and reared to the pupal stage, and the pupae were separated into ♂ and ♀ groups and transferred to separate cages. Three days post eclosion, ratios (pgSIT: WT) of 50 (1:1), 250 (5:1), 500 (10:1), 1000 (20:1), and 2000 (40:1) age-matched pgSIT mature ♂ adults were caged with these mature ♂’s from 100 selected eggs. One hour later, mature ♀’s from 100 selected eggs were transferred into each cage. All adults were allowed to mate for 2 days. ♀’s were blood fed, and eggs were collected. Eggs were counted and stored for 4 days to allow full embryonic development. The remaining eggs were hatched to measure hatching rates and to screen for the possible presence of transformation markers. The hatching rate was estimated by dividing the number of hatched eggs by the total number of eggs. This procedure continued for all subsequent generations.
Mathematical modeling. To model the expected performance of pgSIT at suppressing and eliminating local Ae. aegypti populations, we used the MGDrivE simulation framework50. This framework models the egg, larval, pupal, and adult mosquito life stages with overlapping generations, larval mortality increasing with larval density, and a mating structure in which ♀’s retain the genetic material of the adult ♂ with whom they mate for the duration of their adult lifespan. The inheritance pattern of the pgSIT system was modeled within the inheritance module of MGDrivE. This module includes parameters for ♀ pupatory success (as a proxy for viability), and ♂ fertility and mating competitiveness. We simulated a randomly-mixing Ae. aegypti population of 10,000 adults and implemented the stochastic version of MGDrivE to capture chance effects that are relevant to suppression and elimination scenarios. Density-independent mortality rates for juvenile life stages were calculated for consistency with the population growth rate in the absence of density-dependent mortality, and density-dependent mortality was applied to the larval stage following Deredec et al.51. Sensitivity analyses were conducted with “window of protection” (mean duration for which the Ae. aegypti population is suppressed by at least 90%), probability of Ae. aegypti elimination, and reduction in cumulative potential for transmission as the outcomes, and model parameters sampled describing the construct (gRNA cutting rate and maternal deposition frequency), release scheme (number, size and frequency of releases), ♀ viability, and ♂ fertility and mating competitiveness. Model parameters were sampled according to a latin-hypercube scheme with 216 parameter sets and 20 repetitions per parameter set. Regression models were fitted to the sensitivity analysis data using multi-layered perceptrons as emulators of the main simulation in predicting the window of protection outcome. Subsequent simulations focused on the most sensitive parameters from the sensitivity analysis (those describing the release scheme, ♀ viability, and ♂ fertility) and were explored via factorial simulation design. Complete model and intervention parameters are listed in Supplementary Table 10–11.
Transgene detection method in mosquitoes. SENSR assays were performed as described previously52 using a two-step nucleic acid detection protocol. Target sequences were first amplified in a 30 min isothermal preamplification reaction using recombinase polymerase amplification (RPA). RPA primers were designed to amplify 30 bp gRNA spacer complement regions flanked by 30 bp priming regions from the transgenic elements (Cas9 coding sequence or gRNA scaffold for Cas9 or gRNAdsx,ix,βTub construct, respectively) while simultaneously incorporating a T7 promoter sequence and “GGG'' into the 5′ end of the dsDNA gene fragment to increase transcription efficiency53. The RPA primer sequences and synthetic target sequences can be found in Supplementary Table 12. The RPA reaction was performed at 42°C for 30 min by using TwistAmp Basic (TwistDx #TABAS03KIT). The final conditions for the optimized (28.5% sample input) RPA protocol in a 10-µL reaction are as follows: 5.9 µL rehydration buffer (TwistDx), 0.35 µL primer mix (10 µM each), 0.5 µL MgOAc (280 mM), and 50 ng of genomic DNA. The RPA reaction was then transferred to a second reaction, termed the Cas Cleavage Reaction (CCR), which contained a T7 polymerase and the CasRx ribonucleoprotein. In the second reaction, in vitro transcription was coupled with the cleavage assay and a fluorescence readout using 6-carboxyfluorescein (6-FAM). A 6-nt poly-U probe (FRU) conjugated to a 5′ 6-FAM and a 3′ IABlkFQ (Iowa Black Fluorescence Quencher) was designed and custom-ordered from IDT54. Detection experiments were performed in 20-µL reactions by adding the following reagents to the 10 µL RPA preamplification reaction: 2.82 µL water, 0.4 µL HEPES pH 7.2 (1 M), 0.18 µL MgCl2 (1 M), 1 µL rNTPs (25 mM each), 2 µL CasRx (55.4 ng/µL), 1 µL RNase inhibitor (40 U/µL), 0.6 µL T7 Polymerase (50 U/µL), 1 µL gRNA (10 ng/µL), and 1 µL FRU probe (2 µM). Experiments were immediately run on a LightCycler 96 (Roche #05815916001) at 37°C for 60 min with the following acquisition protocol: 5 sec acquisition followed by 5 sec incubation for the first 15 min, followed by 5 sec acquisition and 55 sec incubation for up to 45 min. Fluorescence readouts were analyzed by calculating the background-subtracted fluorescence of each timepoint and by subtracting the initial fluorescence value from the final value. Statistical significance was calculated using a one-way ANOVA followed by specified multiple comparison tests (n = 4). The detection speed of each gRNA was calculated from the Half-Maximum Fluorescence (HMF) 54. Raw SENSR assays data can be found in Supplementary Table 13.
We generated 12 gRNAs, six target regions within Cas9 or gRNA transgenes, to assess the detection capabilities in a SENSR reaction for each selected gRNA using Cas9 and gRNAdsx,ix,βTub plasmids for detection of Cas9 and gRNAdsx,ix,βTub transgenic sequences, respectively. After identifying the best gRNA for detection of Cas9 and gRNA sequences, we run SENSR of DNA preps made from Ae. aegypti pgSIT♀Cas9 and WT mosquitoes harboring both Cas9 and gRNAdsx,ix,βTub constructs at different ratios of pgSIT♀Cas9 to WT mosquitoes: 1:0; 1:1; 1:5; 1:10; 1:25; 1:50; 0:1; and No Template Control (NTC). We challenged the best gRNA (A4 and B3) in a SENSR assay and found that both gRNA detected transgenic DNA even at the lowest dilution and no signal was observed in both controls, only WT mosquitoes and NTC.
Statistical analysis. Statistical analysis was performed in JMP8.0.2 by SAS Institute Inc and Prism9 for macOS by GraphPad Software, LLC. We used at least three biological replicates to generate statistical means for comparisons. P values were calculated using a two-sided Student’s t test with equal variance. The departure significance for survival curves was assessed with the Log-rank (Mantel–Cox) texts. Statistical comparisons of WBF for each sex were made using ANOVA on Ranks tests, followed by pairwise Wilcoxon tests with a Bonferroni correction applied for multiple comparisons. Chi-squared test was used to assess significance of departure in phonotaxis assay. All plots were constructed using Prism 9.1 for macOS by GraphPad Software and then organized for the final presentation using Adobe Illustrator.