Interestingly, but unexpectedly our study ³¹ also revealed the presence of 537 putative transcripts encoding plant like proteins associated with the sugar fed library, but absent in the blood fed salivary transcriptome database ( Figure 1A; Supplementary material ST1). The surprising discovery of these transcripts, which we labeled as plant like transcripts (PLTs), raised several puzzling, but arguable questions that prompted us to clarify: (i) whether the PLTs are of mosquito origin; (ii) if they are expressed in the mosquito tissues and/or other developmental stages (iii) if expressed in mosquito, what is the possible evolutionary and functional correlation of these transcripts in feeding and (iv) whether these transcripts have any molecular relationship to plant-mosquito-microbe interactions/symbiotic associations. To uncover the molecular nature and possible functions of the putative PLTs, in the present investigation we performed a systematic and comprehensive analysis of PLTs, revealing a unique case of the massive transfer of HGTs from plant to mosquito.

First, to confirm the nature of the PLTs’ origin, we did a deep enquiry with technical staff and confirmed that under standard rearing facilities, mosquitoes are never exposed to any plant material. To further rule out the possibilities of any contamination, we separately maintained the experimental mosquitoes as detailed in the methodology section. For technical validation of the PLTs’ origin, we conducted a series of experiments: (i) in two independent experiments, we examined and verified the RT-PCR based expression of at least 10 selected PLTs ( Figure 1B; Supplementary material S1A), in the salivary glands of adult male and female mosquitoes; (ii) interestingly, we also observed that PLT expression is not only restricted to the mosquito tissues, but is also expressed during the aquatic developmental stages viz. egg, larva, and pupa of the laboratory reared mosquitoes ( Figure 1B). Our relative gene expression analysis revealed that PLTs are more dominantly expressed in the egg, pupa and adult than larval stages ( Figure 1C–E). Although, mosquito egg and pupa stages are metabolically active, and do not take any food material, we suspected that the filter paper being used for mosquito egg laying may be a potential source of environmental contamination carry over. To clarify this doubt, we collected a small piece of moistened filter paper in RNA isolation solution (Trizol) and re-examined PLT expression along with other developmental stages. Absence of any amplification even after 35 PCR cycles, in the filter paper cDNA sample showed no sign of contamination ( Supplementary material S1A); (iii) we also observed positive amplification of selected PLTs through genomic DNA PCR ( Figure 1F); (iv) we further carried out the functional validation of one of the plant homolog PLTs encoding dehydrin protein, by immunoblot analysis as well as immuno-florescence assay ( Figure 2A–G); (v) lastly, from ongoing annotation of another independent transcriptome sequence database originated from non-salivary tissue i.e. olfactory (OLF) tissue of adult female mosquito A. culicifacies (Das De T., Sharma P., Thomas T., Pandey KC., Dixit R. unpublished data), we were able to observe similar PLTs ( Supplementary material S1B); (v) finally to test whether PLT expression is associated with feeding machinery components, we monitored the relative expression of PLTs in four tissues that included salivary glands, midgut, olfactory tissue and hemocytes, collected from 3–4 days old naïve adult female mosquitoes, by real-time PCR ( Figure 5A). Interestingly, we not only observed that PLTs are dominantly expressed in the tissues associated with mosquito feeding machinery (olfactory tissue, salivary gland and midgut), but also noticed a significant down regulation in response to blood meal in the salivary glands ( Figure 5A,B), evidencing that the mosquito genome may code plant like proteins.

Dehydrins are a group 2 member of late embryogenesis abundant (LEA) proteins, originally identified from land plants, and known to be associated with desiccation (water stress) tolerance ³³ . In fact LEA proteins were thought to be restricted to plants and other lower eukaryotes viz. cyanobacteria, algae, but now they have also been identified in other animals including insects ³⁴ . Dehydrins are evolutionarily conserved proteins acclimated to low-temperatures (LT) that allow efficient tolerance to drought and cold stress among photosynthetic as well as some non-photosynthetic organisms such as yeast ^{35– 38} . Dehydrins are characterized by lysine rich conserved K-segment comprising consensus amino acid sequence EKKGIMDKIKEKLPG, towards the C-terminus that may be repeated many times to encode 9 -200 kDa protein ^{39– 42} . This unique feature renders these proteins cationic, providing cryoprotective activity towards freezing sensitive enzymes ⁴³ . The biochemical characterization of a novel cryoprotective protein in freeze-tolerant Eurosta solidaginis larvae shows dehyrin like activities ⁴⁴ , but a true homolog of dehydrin is yet to be verified.

Mosquito dehydrins have not been reported so far, though a putative transcript AGAP000328 has been predicted from the mosquito A. gambiae genome, carrying (PF00257 domain) a signature of dehydrin like proteins ( Supplementary material S2). Finding PLT encoding proteins associated with dehydration stress e.g. dehydrin, aquaporin, expansin etc. encouraged us to further examine their possible functions in the mosquito A. culicifacies. A comprehensive molecular analysis of the identified transcript AcDehydrin showed 100% identity to the plant dehydrin, having two conserved lysine rich K-segments ( Figure 2A–D). In our relative gene expression analysis, we observed a constitutive expression of AcDehydrin, throughout the aquatic developmental as well as adult stages of the mosquito ( Figure 2E), indicating that AcDehydrin transcript may encode a putative functional protein.

For functional validation of AcDehydrin protein, we examined the developmental expression of the dehydrin protein through immuno-blotting assay using rabbit antiserum containing anti-dehydrin antibody, raised against conserved K-segment sequence TGEKKGIMDKIKEKLPGQH ( Figure 2D) of plant dehydrin ⁴⁰ (kind gift from Dr. Timothy Close). In our experiments we used wheat seedling protein sample as positive reference control. The anti-dehydrin antibody not only recognized the expected (28, 53 and 62 kDa) protein band in the wheat samples ⁴⁵ , but also identified at least two equivalent proteins (28 and 62kDa) abundantly expressing in different mosquito developmental stages viz. egg, adult male and female mosquitoes ( Figure 2F). Additionally, we were also able to observe multiple isoforms ranging from (~10 ->70 kDa) expressing at low level in different developmental stages, an expected unique feature of dehydrin to form macromolecular structures ^{39– 42} . Finally, immuno-florescence assay not only corroborated the abundant expression in the egg, but also suggested that mosquito encoded AcDehydrin protein may play a crucial role in the stress tolerance and survival of the embryo in the egg ( Figure 2G).

Like other LEAs, dehydrins accumulate to high amounts in plant embryos, but remain undetectable in other vegetative tissues until their exposure to dehydration stress. The stress exposure results in their rapid induction and binding to multiple proteins, probably through intramolecular hydrogen bonding to protect tissue damage from dehydration/cold stress ⁴⁶ . In fact, we also find another key transcript, encoding a putative protein named expansin, a member of plant cell wall-loosening proteins. These proteins are known to be involved in cell enlargement and developmental processes requiring cell-wall modification ⁴⁷ . Like dehydrin, a real-time PCR analysis of expansin also showed dominant expression in the egg, as compared to other developmental stages ( Supplementary material S2B). Taken together, we hypothesize that the mosquito A. culicifacies may have survival benefits of cold stress tolerance as well as developmental regulation, similar to plants. Future studies involving dsRNA mediated gene silencing may unravel molecular and functional relationship of the PLTs controlling feeding and adaptation phenotypes in the mosquito ^{48, 49} .

Next, to understand the possible evolutionary relationship we performed an extensive phylogenomic analysis of a few selected transcripts. To do this, first we retrieved and analyzed all plant-homolog putative transcripts (537 PLTs), and performed an extensive BLASTX analysis against either the NR database or an insect specific database at NCBI ( http://blast.ncbi.nlm.nih.gov/), and further characterized three categories of transcript(s) (i) one transcript: encoding highly conserved alpha-tubulin (cytoskeleton associated protein), showing highest identity to plant and insect (>95% and 85–90%, respectively); (ii) two transcripts: encoding aquaporin (water channel membrane protein) ⁵⁰ and active site of the cysteine protease (protein chewing enzyme) ⁵¹ showing highest identity to plant (>90%) and 40–52% identity to insect (iii) two transcripts: encoding dehydrin (cold stress response protein) ⁴⁶ and expansin (plant cell wall loosening protein) ⁵² only matched to plants, but remained unmatched to any insect database ( Supplementary material S3 A–C).

The above results prompted to follow up the associated evolutionary consensus, favoring plant-mosquito relationship: a parallelism setting where different species from unrelated taxa faces the common selective pressure ⁵³ . Initial multiple sequence alignment analysis revealed significant heterogeneity (substitution/deletion) of amino acid residues, but also indicated unique conservation of insect or plant specific residues within the mosquito A. culicifacies, result a clade formation with plant species ( Figure 3A,B; Supplementary material S4). Subsequently, we also tested whether the evolution of common traits from unrelated taxa owing to similar selection pressure favors adaptive significance.

A maximum likelihood (ML) estimation was applied to calculate and compare the site-wise likelihood (∆SSLS) values between the two hypotheses, i.e. mosquito-mosquito species evolution (H ₀) and mosquito-plant convergent adaptive evolution (H ₁), for the selected PLTs. The site-wise log likelihood plot indicator, i.e. divergence towards negative (∆SSLS) was compared with LRT (likelihood ratio test), using the parametric bootstrap at 1000 replicate analysis (cut off p-value 5%). Final data analysis and comparison statistics favored the convergent hypothesis ⁵⁴ , demonstrating that mosquito A. culicifacies PLTs followed a convergent model favoring (H ₁), an adaptive evolution for sugar feeding associated functional relationship with plants ( Figure 3C; Supplementary material S4). Our analysis also supports the previous observations noted for the evolution of echolocating gene clusters among bats and bottlenose dolphins ⁵⁵ . Additionally, the predicted 3D structural analysis revealed fine conservation of the active functional domains in the mosquito and plant proteins e.g. cysteine protease ( Figure 3D; Supplementary material S5). From these studies, we concluded that mosquito feeding associated genes are not only evolving actively, but also acquiring new genes (e.g. dehydrin, expansin), to adapt successfully over the plant host.

Insect-plant association represents one of the most dominant interactions over millions of years ^{56– 58} . These interactions are thought to play an important role in the co-evolution of molecular effector arms, enabling effective adaptation over each other ⁵⁹ . Uncovering of the molecular mechanisms of the herbivore insect-plant interaction has greatly facilitated the design of molecular strategies to save the valuable crops from insect pests ^{60– 63} . However, such studies have not given special attention to mosquitoes. From the unexpected findings of the mosquito PLTs, we interpreted that either studies in relation to the sugar feeding associated biology have largely been ignored ²⁸ or the mosquito A. culicifacies may have evolved with more complex genetic architecture favoring evolution of several environmentally-guided traits viz. carbon metabolism; light mediated photo conditions for mating, feeding, survival etc. Therefore, to predict sugar metabolism associated molecular and functional relationships of salivary PLTs, initially we analyzed all the putative PLTs against three databases (Reactome, KEGG, and Biocycles) annotated for Arabidopsis thaliana, using KOBAS online software, version 2.0 ( http://kobas.cbi.pku.edu.cn/home.do).

Notably, we observed that 18 transcripts encoding proteins related to at least five Biocyclic pathways linked to photosynthetic organelles viz. plastid in plants ( Supplementary material S6, T1). To verify the predicted ‘plastid’ related salivary transcripts, Fisher’s exact test was performed using BLAST2GO, revealing a pool of 11 transcripts differentially expressed in the sugar fed mosquitoes (Fisher test p<0.001; Supplementary material S6B) encoding important enzymes/proteins, associated with one of the key pathway “Carbon fixation in Photosynthetic Organisms” ( Figure 4A). Further, we also identified four unique salivary transcripts encoding different enzymes linked to three other secondary metabolite synthesis pathways, namely: ‘Trepenoid Backbone Biosynthesis’ (4-hydroxy-3-methylbut-2-enyl diphosphate reductase/E.C.1.17.1.2, LYTB); ‘Carotenoid Biosynthesis’ (Phytoene Synthase/E.C.2.5.1.32, PS); and ‘Flavonoid Biosynthesis’ (3-dioxigenase/E.C.1.14.11.9 & 3’ beta-hydroxylase/E.C.1.14.13.88) pathways restricted to plants ( Supplementary material S7). A comprehensive molecular and phylogenetic analysis of a few selected transcripts, encoding an enzyme 4-hydroxy-3-methylbut-2-enyl diphosphate reductase/E.C.1.17.1.2 (LYTB) and phytoene synthase/E.C.2.5.1.32 (PS) exclusively revealed a unique evolutionary relationship to the cyanobacteria, algae, plants and aphid Acyrithosiphon pisum ( Figure 4B,C; Supplementary material S7). In fact, during its early development mosquito larvae start to feed on diverse micronutrients e.g. bacteria, algae, fungi etc., and switch to feed on nectar sugars in adult mosquito stage. Thus, it could be possible that a long association and regular microbe-mosquito-plant interactions ^{64, 65} , might have favored insects (mosquitoes) to adapt, feed, and digest sugar as well as selective synthesis of secondary metabolites/pigments, essential for specific phenotype e.g. visual pigmentation/dark body coloration ⁶⁶ . A recent study on light-induced ATP synthesis from the chloroplastid-like carotene pigments in Acyrithosiphon pisum, a plant sap sucking aphid, provides the first molecular evidence that the aphid genome may carry plant like photosynthetic machinery components ⁶⁷ . A fungal mediated lateral HGT mechanism has been proposed for the evolution of the carotenoid biosynthesis gene in this aphid ²⁰ .

In nature, mosquitoes are regularly exposed to various environmental factors which have adverse effects on their reproductive success, longevity & vector competence ⁶⁸ . Gut bacterial endosymbionts are known to play a part in several functions including food digestion, metabolism, reproduction and immunity ⁶⁹ . Our recent metagenomic analysis of salivary microbiome identified several unique bacterial phyla, including Chlorobium, Cyanobacteria, Nitrospira and other phototrophic bacteria associated with salivary glands ( Figure 4D), but absent in the gut of laboratory reared 3–4 days old adult female mosquitoes A. culicifacies ⁷⁰ . Indirectly, the above findings further support the hypothesis that mosquitoes may have feeding associated distinct plant like molecular machinery components, partly shared by the residing symbiotic bacterial community for diverse carbon/nitrogen rich plant sugar source metabolism. For example, finding of prominent salivary associated Acidobacteria (2.4%), may facilitate the utilization of plant polymer viz. cellulose/xylan sugars of diverse origin ⁷¹ , as reported in the gut of the wood feeding larvae of Huhu Beetle ( Prionoplus reticulari) ⁷² .

The observation of a large pool of chloroplast and nuclear encoded plant genes in the mosquito transcriptome supported the previous finding of similar gene transfer of photosynthetic machinery components from algae to mollusk ¹⁹ . In addition to this, molecular analysis of PLTs also revealed a plant related class of secondary metabolites (see above) and immune genes i.e. Remorin (anti-fungal); osmotin/thaumatin (anti-fungal); and Vicilin (Antimicrobial) ( Supplementary material ST-1), a similar finding of active genes in the aphid genome ⁷³ . Although, in case of algal-mollusk or Aphid-plant interaction studies, the role of microbes is yet to be established, however, our metagenomic analysis provides initial evidence that tissue associated microbial flora may also share and facilitate optimal function. Thus, we believe that the accumulating data of genetic material transfer within metazoans are still at a premature stage, but emerging evidences strongly suggest that acquisition and retention of desired active functional genes for beneficial traits, may favor improved survival and adaptation values in changing ecologies ⁴ .

With the current available data, including the present investigation, we hypothesize that HGTs in metazoans may also play an important role in the evolution and acquisition of beneficial traits that facilitate feeding and survival adaptation over diverse ecologies. This hypothesis is further strengthened by our following new observations: (a) that PLTs expression seems to be restricted to the tissues, i.e. the feeding machinery components that facilitate digestion and metabolism, e.g. salivary glands, midgut olfactory tissues in case of the mosquito ( Figure 5); (b) absence of PLTs from other non-digestive tissues, e.g. hemocyte (mosquito blood cells); (c) the finding of dominantly associated unique bacterial species to the mosquito digestive tissues viz. salivary gland and midgut ⁷⁰ , e.g. Acidobacteria (sugar metabolism); Agromonas bacteria, a soil oligotroph (nitrogen fixing bacteria) that usually grow at extra low nutrient environments of the paddy field, complementing the high larval density of the mosquito A. culicifacies in paddy fields of the rural India ⁷⁴ . Indeed, Agromonas has been previously isolated from paddy fields ⁷¹ ; but largely remain unidentified from any insect species so far.

Furthermore, in mosquitoes it has long been accepted and proven that a significant variation exists in the chromosomal DNA as well as genome size within Anopheline and other mosquito species ^{75, 76} , but how these variations differentially affect the mosquito biology viz. behavior, physiology, immunity and vectorial capacity etc., are poorly understood at the molecular level.

Tissue collection: For the desired tissues viz. salivary glands, midgut, hemocytes, olfactory tissues collection, we followed essentially the same protocol as established for our recent metagenomic and transcriptomic studies ^{31, 70} . Briefly, prior to dissection, 3–4 day old sugar or blood fed adult female mosquitoes were surface sterilized using 70% ethanol for 1 min followed by dissection in a sterile water drop on a microscopic slide in a sterile working area under laminar flow. Sterile entomological needles/forceps were used to manually pick up and collect the tissues in the pre-sterilized 1.5 ml Eppendorff tubes containing 50 μl Trizol solution. For the hemocyte collection, a flushing method was opted for as described previously ⁸⁰ . Briefly 2–3 μl of Schneider’s (RPMI): FBS: citrate buffer (60:10:30) was injected into the lateral wall of the mesothorax of cold anesthetized mosquitoes, followed by flushing out the diluted hemolymph with an additional 3–5 μl of Schneider’s (RPMI), by clipping of the last abdominal segment. The diluted hemolymph was directly collected by pipette in Trizol. For the egg collection, a clean fine art paint brush was used to scrape the eggs from moistened filter paper, rinsed with sterile water and collected in the Trizol. The other aquatic developmental stages viz. larva, pupae were also manually picked up with a Pasteur pipette, washed with sterile water twice and collected in Trizol for RNA isolation.

RNA isolation, cDNA preparation and PCR analysis: The desired tissues viz. salivary glands, midgut and hemocyte ⁷⁸ or the whole body, were collected in Trizol. Total RNA was isolated using standard Trizol method, followed by first-strand cDNA synthesis using Oligo-dT or Random Hexamer primers (Verso kit). For differential expression analysis, routine RT-PCR and agarose gel electrophoresis protocols were used. Relative gene expression was assessed by QuantiMix SYBR green dye (Biotool Biolabs, Madrid, Spain) in Eco-Real-Time PCR Machine (Illumina). PCR cycle parameters involved an initial denaturation at 95°C for 5 min, 40 cycles of 10 s at 95°C, 15 s at 55°C, and 22 s at 72°C. Fluorescence readings were taken at 72°C after each cycle. A final extension at 72°C for 5 min was completed before deriving a melting curve, to confirm the identity of the PCR product. Actin gene was used as an internal control in all qPCR measurements, where minimum two technical replicates were used in each real-time experiment. To better evaluate the relative expression, each experiment was performed in three independent biological replicates. The relative quantification results were normalized with internal control Actin gene and analyzed by 2 ^–ΔΔCt method ⁷⁹ .

PLTs identification and phylogenomic analysis: In an attempt to clarify and improve the functional annotation of a cluster of unique sequences encoding plant like proteins, unexpectedly observed from our recently sequenced salivary transcriptomes ³¹ , we performed a comparative analysis for both the sugar fed as well as blood fed salivary transcriptomic databases. Initially, to do this we did a species distribution analysis and manually sorted and catalogued the sequences that best match to the plants from the FASTA file. The shortlisted transcripts were subjected to a similarity search against NCBI's NR database using the BLASTx algorithm ⁸⁰ , with a cut-off E-value of ≤10 ⁻³ using BLOSUM62 matrix as well as GO annotation/Interproscan analysis using BLAS2GO ⁸¹ . Biocyclic pathway analysis for PLTs KOBAS online ( http://kobas.cbi.pku.edu.cn/home.do) software, version 2.0 ⁸² . Following primary BLASTX analysis, the reference sequences from the selected top hits were retrieved and edited for subsequent analysis in the FASTA format. Multiple sequence alignment was performed using ClustaX2, version 2.0 ⁸³ . The CLC Sequence viewer ( http://www.clcbio.com) software (version 6.9.1) was used for better quality graphics. The phylogenetic relationship was inferred through MEGA5.1 ( http://www.megasoftware.net/) software. The evolutionary history was inferred using the Neighbor-Joining method, and percentage of replicate trees in which the associated taxa clustered together in the bootstrap test (1000 replicates). The evolutionary distances were computed using the p-distance method, presented in the units of the number of amino acid differences per site. A work flow for the Phylogenomic analysis has been presented in the Supplemental material S6. The following major steps were followed:

Modeling procedure & 3D structural prediction analysis: All structures of representative protein were retrieved from the Protein Data Bank ( www.rcsb.org) and aligned using the structure alignment program STAMP4.0 ⁸⁶ . Models using all four structures as templates were generated using Modeller9 version 10 ⁸⁷ . 3D representation of the model was prepared in VMD version 1.9 (Visual Molecular Dynamics tool) ⁸⁸ .

Genomic DNA isolation & PCR: For the genomic DNA extraction, a total of five adult female mosquitoes, decapitated with head and wing, were collected in extraction buffer and processed as described earlier ⁷⁰ . All the PCR amplification conditions and parameters were identical as described above for RT-PCR analysis.