Molecular signatures of sexual communication in the phlebotomine sand flies

Abstract

Phlebotomine sand flies employ an elaborate system of pheromone communication wherein males produce pheromones that attract other males to leks (thus acting as an aggregation pheromone) and females to the lekking males (sex pheromone). In addition, the type of pheromone produced varies among populations. Despite the numerous studies on sand fly chemical communication, little is known of their chemosensory genome. Chemoreceptors interact with chemicals in an organism's environment to elicit essential behaviors such as the identification of suitable mates and food sources. Thus, they play important roles during adaptation and speciation. Major chemoreceptor gene families, odorant receptors (ORs), gustatory receptors (GRs) and ionotropic receptors (IRs) together detect and discriminate the chemical landscape. Here, we annotated the chemoreceptor repertoire in the genomes of Lutzomyia longipalpis and Phlebotomus papatasi, major phlebotomine vectors in the New World and Old World, respectively. Comparison with other sequenced Diptera revealed a large and unique expansion where over 80% of the ~140 ORs belong to a single, taxonomically restricted clade. We next conducted a comprehensive analysis of the chemoreceptors in 63 L. longipalpis individuals from four different locations in Brazil representing allopatric and sympatric populations and three sex-aggregation pheromone types (chemotypes). Population structure based on single nucleotide polymorphisms (SNPs) and gene copy number in the chemoreceptors corresponded with their putative chemotypes, and corroborate previous studies that identified multiple populations. Our work provides genomic insights into the underlying behavioral evolution of sexual communication in the L. longipalpis species complex in Brazil, and highlights the importance of accounting for the ongoing speciation in central and South American Lutzomyia that could have important implications for vectorial capacity.

Fig 1. Chemoreceptors in the phlebotomine sand flies.

(A) Chemoreceptor repertoire size in the sand flies L. longipalpis (Llon) and P. papatasi (Ppap) compared with An. gambiae (Agam), C. quinquefasciatus (Cqui), Ae. aegypti (Aaeg), Ma. destructor (Mdes), Mu. domestica (Mdom), G. morsitans (Gmor) and D. melanogaster (Dmel). Species tree estimated using OrthoFinder with the multiple sequence alignment option and 252 single-copy orthologs. (B) Phylogenetic analysis of chemoreceptors in five Diptera revealed a large taxonomically-restricted clade comprising over 80% of the L. longipalpis and P. papatasi ORs (highlighted in green). (C) Several smaller lineage expansions are evident in the GRs, while (D) only one IR lineage (Ir7c) was expanded in the sand flies. Phylogenies were estimated using L. longipalpis, P. papatasi, A. gambiae, M. destructor and D. melanogaster protein sequences aligned with ClustalX. The JTT model of protein substitution and Maximum Likelihood method in RAxML v.8.2.4 were used for tree estimation. The trees are rooted at the branch leading to Orco, the CO₂ receptors, and the Ir25a and Ir8a clades for ORs, GRs and IRs, respectively. Branch support based on 500 bootstrap replications. Scale bars indicate the number of amino acid substitutions per site.

Fig 2. Population structure of 63 L. longipalpis from four sites in Brazil based on SNPs in 245 chemoreceptor genes.

(A) Geographic distribution, sex-aggregation pheromones (9MGB, 3MαH, sobralene) and copulatory songs based on previous studies of L. longipalpis in Brazil (16). (B) Individuals from Sobral (with 2 spots: S2S), Marajo and Lapinha formed discrete clades, while individuals from Sobral (with one spot, S1S) and Jacobina split into two clades each. Unrooted tree estimated using SNPs in the exons of all 100 single-copy orthologs and neighbor-joining method in Tassel v5.2.57. (C) Principal component analysis (PCA) was used to identify loci associated with population structure and conducted using pcadapt (explained variance EV). The first two principal components accounted for 73.6% of the total variation and grouped individuals into four clusters. Putative chemotypes were assigned based on previous studies and PCA clustering patterns. (D) Principal components 1–5 and genes with the highest number of SNPs based on component-wise outlier analysis in pcadapt. (E) Ancestry proportions within individual sand flies for ADMIXTURE models from K = 1 to K = 7 ancestral populations. Each vertical bar represents the proportion of ancestry within a single individual, with colors corresponding to ancestral populations. Data are the average of the major q-matrix clusters derived by CLUMPAK analysis. (F) Violin plot of ADMIXTURE cross-validation error for each of 30 replicates for each K value from 1 to 7.

Fig 3. Visualization of reads aligned to chemoreceptor loci revealed variation in coverage that indicated potential copy number variation.

(A) For example, JAC01 Or109 had much deeper coverage than Or98, while Or94 had reads mapped only at the end of the second exon. The Tablet software program was used for visualization of the mapped reads. (B) To quantify these differences for comprehensive analysis of all 245 chemoreceptor loci, we calculated background-normalized sequencing depth of each gene using the modal depth across the exons in all protein coding genes. (C) A central tendency of ~2 is expected for single copy genes with two intact alleles. Normalized depth was rounded to the nearest whole number as a proxy for copy number (CN). (D) The number of intact chemoreceptors (CN≥2) in all individuals of a chemotype ranged from 141 (Sobralene) to 170 (3MαH). (E) The mean number of absent (CN = 0) genes differed among all chemotypes (P <0.001), with Sobralene individuals having the most and 3MαH individuals having the fewest. (F) Accordingly, the number of single-copy (CN = 2) genes differed among all three chemotypes (P <0.001), with 3MαH individuals having the largest number and Sobralene individuals having the fewest. (G) The number of duplicated genes (CN>2) differed only between 3MαH and sobralene individuals.

Fig 4. Relationships among 63 L. longipalpis from four sites in Brazil (Lapinha Cave, Sobral, Jacobina and Marajó) based on copy number (CN) of 245 chemoreceptor genes.

(A) Hierarchical analysis of gene CN clustered individuals according to their putative chemotype as determined previously using SNPs. A heatmap of CN illustrates the large number of genes with CNV among the odorant receptors. (B) PCA of gene CN showed a similar pattern, wherein Sobral 2S (2 spots), Marajó and six Jacobina clustered together (Sobralene); seven Sobral 1S (1 spot) and eight Jacobina clustered together (3MαH); and Lapinha Cave and six Sobral 1S (one spot) were in separate clusters (9MGB), which is consistent with the genetic differentiation observed by Hamilton et al. (2005) thus leading them to classify these as different chemotypes (9MGB and 9MGB+) due to the larger quantity of 9MGB produced by males from Lapinha [10].

Fig 5. V_ST was used to identify the most differentiated genes based on CNV.

Heatmaps of copy number (CN) and Manhattan plots of V_ST between (A) sobralene and 9MGB, (B) sobralene and 3MαH, and (C) 9MGB and 3MαH (V_ST >0.5 highlighted in red). The dashed lines indicate the threshold for significance (0.99) based on 1,000 permutations. Heatmaps illustrate CN of genes with V_ST>0.5. Of the 60 genes with V_ST>0.5 in all three pairwise analyses, 43 were ORs, 16 were GRs and only one was an IR.