Adaptation of Fig Wasps (Agaodinae) to Their Host Revealed by Large-Scale Transcriptomic Data

Simple Summary Research on fig wasps has made a considerable contribution to the understanding of insect–plant interactions. However, the molecular mechanisms underlying fig wasp host specificity are poorly understood. This study reports on a relatively large-scale transcriptomic dataset of 25 fig wasp species. We outline potential genetic mechanisms underlying the specific host adaptation by investigating changes in a gene family, in evolutionary rates, and in genes under positive selection. The transcriptome datasets reported here (1) provide new insights into the evolutionary diversification and host specificity of fig wasps and (2) contribute to a growing dataset on fig wasp genomics. Abstract Figs and fig wasps are highly species-specific and comprise a model system for studying co-evolution and co-speciation. The evolutionary relationships and molecular adaptations of fig wasps to their fig hosts are poorly understood, and this is in part due to limited sequence data. Here, we present large-scale transcriptomic datasets of 25 fig wasp species with the aim of uncovering the genetic basis for host specificity. Our phylogenetic results support the monophyly of all genera associated with dioecious figs, and two genera associated with monoecious figs, Eupristina and Platyscapa, were revealed to be close relatives. We identified gene loss and gain, potentially rapidly evolving genes, and genes under positive selection. Potentially functional changes were documented and we hypothesize as to how these may determine host specificity. Overall, our study provides new insights into the evolutionary diversification of fig wasps and contributes to our understanding of adaptation in this group.


Introduction
Intimate inter-specific interactions are pervasive in nature. Species embedded within these complex networks have consumed each other, provided provisions for each, and competed over ecological and evolutionary time [1].
Evidence for co-evolution in the strict sense is rare [2], but insects and plants clearly form part of each other's selective landscapes. In some cases, reciprocal selection appears to trigger increased rates of diversification [3]. It is likely that "diffuse" co-evolution acting among groups of individuals is more widespread [4]. In any case, a mechanistic understanding of how selection shapes genomic architecture is hard to achieve in a multispecies setting. Simpler systems that offer a degree of phenotypic matching therefore represent attractive opportunities.
Obligate mutualisms between plants and their pollinators are useful here because the fitness of each partner is closely linked, and trait mismatch or association with an incompatible partner can be expected to result in greatly reduced fitness [5][6][7]. Reciprocal selection likely results in a strong signal, with much of the background "noise" associated with complex lifecycles involving multiple partners removed.
Here we focus on the interaction between figs and their pollinating wasps. The genus Ficus (Moraceae) is pantropical and consists of over 750 species [8]. Each species of fig is pollinated by one to two wasps from the chalcid family Agaonidae. However, as many as nine pollinators can occur across a single host [9]. The inflorescence of a fig tree is known as the "syconium"; this receptacle has only one entrance (the bract-lined ostiole). Therefore, only specialized fig wasps can enter the syconia to pollinate the female flowers inside (although specialized parasitic fig wasps also exploit this entrance). Once inside, the fig wasp extends its ovipositor into the style of an individual female flower to lay eggs in the ovary within [10,11]. The mutualism between fig and fig wasp has endured for tens of millions of years. Strict mutual adaptation in morphology, behavior, physiology, and development has been documented in this system [12,13].
Olfaction plays an essential role in forming and maintaining the highly specific mutualism between a fig and its corresponding pollinating fig wasps [14,15]. All fig species use chemical cues to attract their specific pollinators, which may include a mixture of compounds or even a single compound, a "private channel" [15,16]. In turn, fig wasps must detect and filter these cues from the surrounding chemical landscape. Besides olfaction, fig wasps can also use short-range tactile and visual cues to determine whether the host is suitable [11,[17][18][19]. Detoxification and the immune response of fig wasps also play an essential role in determining host specificity at the larval stage. Fig wasps are also exposed to pathogens including bacteria, fungi, and nematodes and viruses within syconia (often vectored by the wasps themselves). Therefore, fig wasps should be able to defend themselves against pathogens [20][21][22][23].
Adaptive trait matching has been observed between figs and fig wasps [24]. For example, there is a strong correlation between ovipositor length in wasps and style length in figs. It is reported that figs and fig wasps have a consistent relationship of co-cladogenesis and co-evolution with the same subgenus and same section/subsection of figs.  [13,[25][26][27][28][29]. In total, 19 subgroups of Ficus have been described and can be distinguished based on distinct morphological and reproductive characteristics [28,30]. Consequently, it has been predicted that fig wasps are under selection to adapt to changes in their hosts. For example, the thorax, abdomen, and forefoot of fig wasps in the genus Ceratosolen all have enlarged spiracles to compensate for the low oxygen environment within the fluid of their hosts' syconia [30]. In general, wasp heads are flattened and elongated to fit within the narrow bract lining the entrance to otherwise enclosed figs. The arrangement of bracts is also subgenus-or group-specific, are corresponding adaptations are seen in wasps: when the bracts are linear, the head and mandible appendages are longer and thinner, while the pollinators of figs with bracts that are interlocked into a spiral have heads that are flatter, with a soft area for folding, and the mandible appendages are shorter and firmer [18,31]. We suggest that genomic footprints of selection vary among wasps associated with different lineages of figs. Sexual system (monoecy vs. dioecy) is sometimes correlated with other traits in figs and dioecious species tend to be understory specialists. In contrast, pollinators of monoecious figs disperse using above-canopy winds.
Adult female fig wasps are short-lived and non-feeding, and selection should act to favor individuals capable of quickly locating their host fig using species-specific chemical cues from a distance or other visual and tactile signals from a short range. In general, we predicted that these particular organisms would have a reduced genomic architecture to avoid detection of non-target scents. Genes encoding proteins related to feeding, environmental perception, and the immune response would be expressed and/or show signs of positive selection. Variation in the evolutionary rates of genes and gene families were also predicted to be consistent in closely related species and genera of fig wasps when compared to, for example, allo-generics.
Recent studies building on the first fig wasp genome [6] have used an omics approach to greatly enhance our understanding of how selection leaves footprints in expressed genes. For example, reciprocal selection has shaped signal (volatile organic carbon) and receptor (olfactory and gustatory genes) in fig wasps [32,33], while wasps exposed to their host cues actively alter gene regulation of receptors [34]. Here we took a phylogenetically structured approach and compared baseline gene expression in newly emerged adults among (i) a species complex of five pollinating wasps associated with one host (Valisia); (ii) one species associated with five hosts (Blastophaga sp); (iii) a selection of fig wasps from a single genus spread across several host figs (Ceratosolen); (iv) three additional genera sampled for between one to three species; and (iv) the family Agaonidae. Identifying genes capable of species differentiation and evidence for adaptive evolution at the genomic level will assist with understanding the mechanisms shaping reciprocal adaptation, and phylogenetic estimates should be improved through the consideration of many more markers.
Specifically, we used transcriptomic data from newly emerged adult female wasps and   Table 1). One species, Ficus hirta, is pollinated by nine fig wasp species that occupy distinct geographical regions [9]. Eight of these nine fig wasp species share a recent common ancestor. One species, V. esquirolianae, enters a close relative of F. hirta, F. triloba, in certain parts of its range. In this study, we selected four of the eight pollinators Valisia sp. 1, sp. 2, sp. 7, and sp. 8, and V. esquirolianae as a related species group. In addition, five of the taxa that pollinate F. pyriformis, F. variolosa and F. erecta var. beecheyana, F. formosa, and F. abeli have been identified as a single species by morphology and gene sequencing [35][36][37][38][39]. We considered these to be a monophyletic group.  [9]. It has been reported that the pollinating wasps in F. pyriformis, F. erecta var. beecheyana, F. formosa, F. variolosa, and F. abeli are the same species according to the morphology and DNA sequences [35][36][37][38][39]. Our results also confirmed that they form a single clade. We considered them as one species with five hosts: B.

mRNA-Seq Library Construction, Illumina Sequencing, Assembly, and Annotation
Total RNA of each species was extracted from 50-100 individuals using a modified CTAB method [40]. The quality of RNA was assessed through gel electrophoresis and using an Eppendorf AG 2231 Bio Photometer Plus (Hamburg, Germany). For each sample, a messenger RNA-Seq library was constructed using an Illumina TruSeq™ RNA Sample Preparation Kit (Illumina, San Diego, CA, USA) following the manufacturer's recommendations. The isolation of messenger RNA (mRNA), fragment interruption, complementary DNA (cDNA) synthesis, adaptor ligation, PCR amplification, and RNA-Seq were performed by Novogene Bioinformatics Technology Co., Ltd. (Beijing, China). The RNA-Seq library was sequenced on an Illumina Hiseq 4000 platform.
Raw reads were filtered using Trimmomatic v0.38 [41] and adaptors trimmed. Then, the high-quality clean reads were de novo assembled to get contigs using Trinity v2.8.5 [42] with default parameters. The expression of the contigs was calculated using RSEM v1.3.1 [43]. TransDecoder v5.5.0 (https://github.com/TransDecoder/; accessed on 9 October 2020) was used to predict the coding sequence (CDS) for each isoform of a gene and the isoform sequence with the highest expression was selected as a unigene. Finally, the protein sequences of all the sampled species were compared with the 5991 Benchmarking Universal Single-Copy Orthologs (BUSCO) in the Hymenoptera_odb10 database to evaluate the integrity of transcriptome using BUSCO v4.1.2 [44] with default settings. The raw sequence data have been deposited in the Genome Sequence Archive (GSA) in the National Genomics Data Center, Chinese Academy of Sciences (https://bigd.big.ac.cn/gsa; accessed on 23 March 2021), under accession number PRJCA004756.

Ortholog Identification and Alignment
To find orthologous genes, CDS and protein sequences of 6 Hymenoptera species and 1 Diptera species, Drosophila melanogaster, were collected from NCBI. These species have relatively complete BUSCO and their gene functions have been fully studied. We made a pairwise comparison of the genome or transcriptome protein sequences among these 32 species using the blastp command in diamond v2.0.2.140 [51], and then filtered the blast results using an in-house perl script. Orthologous genes in these filtered data were analyzed using OrthoMCL [52] and clustered with the MCL algorithm [53]. This selection of species includes other Hymenoptera with more complex life histories (including parasitoids and social insects) that also occupy more complex habitats, and they therefore provide a useful baseline for comparison. Based on protein sequence similarity and the mutual best-hit algorithm of all 32 species, orthologous and paralogous gene pairs were identified and clustered into 38,762 orthologous cluster groups (OCGs). Of these, 18,008 OCGs were contained in at least 4 species, and 11,809 OCGs were contained in at least 7 species. These OCGs were selected for codon alignment and downstream analysis.

Phylogenetic Tree and Divergence Time
In total, 661 single/low-copy orthologous genes were found across the 25 species; in 60% of species they were single-copy. The genes with conserved codon sequences of less than 60 bp were filtered by Gblocks v0.91b [54]. The remaining 625 genes were concatenated using an in-house perl script. We estimated the phylogenetic relationships among the species using the maximum likelihood (ML) criterion as implemented in RAxMLv8.2.12 [55].
Two clades in the ML tree, Valisia associated with F. hirta and F. triloba and the five Blastophaga taxa, were recovered with relatively low support (90% and 70% respectively). To help resolve these clades, we generated a dataset including only these taxa and recovered 3189 and 5528 single/low-copy orthologous genes with which to generate an ML phylogeny using the above methods. This phylogeny was better supported by and rooted to Ceratosolen gravelyi or Ceratosolen constrictus. For these two clades, the overlapping 625 single-copy homologous genes were analyzed using the baseml command in PAML v4.9i [56] to generate individual gene phylogenies. These 625 individual phylogenies were integrated into a final phylogenetic tree using an in-house perl script with the branches standardized by the median method.
The divergence time between species was calculated using mcmtree in PAML based on the final topology [56]. The following nodes were calibrated: 127-192 MYA between Polistes canadensis and Apis mellifera, 172-241 MYA between Nasonia vitripennis and Acromyrmex echinatior, and 308-366 MYA between Apis mellifera and Drosophila melanogaster. This was undertaken using constraints derived from "Time Tree" (http://timetree. org/; accessed on 20 October 2020) and the nodes were used to convert relative divergence times to absolute divergence times. Mcmctree was run twice to test convergence (estimates between runs varied by <1%).

Gene Family Expansion and Contraction
In this study, the genes obtained through transcriptome sequencing were all functional genes expressed across study species. We identified gene families using CAFE [57], which employs a random birth and death model to study gene gains and losses in gene families across a user-specified phylogeny. The global parameter λ, which describes both the gene birth (λ) and death (µ = −λ) rate across all branches in the tree for all gene families, was estimated using maximum likelihood. A conditional p-value was calculated for each gene family. Only families with conditional p-values less than a fixed threshold (0.01) were considered as having an accelerated rate of gain or loss. The expanded and contracted gene families for each branch were tested for enrichment in the GO and KEGG databases. We compared the proportions of enrichment from GO and KEGG between the foreground and background gene families. Significance was tested using Fisher exact tests.
Five species included in our phylogenetic tree had whole genome data available as well as information relating to gene family expansion and contraction [6]. We summarized the number of expanded and contracted gene families expressed for each species in comparison to our phylogenetic tree and assessed how transcriptome data compared to genome level sequences [6]. When considering whole genomes, any change in the quantity of a gene family was regarded as either an expansion or contraction. Here, we used the same standard to summarize the number of such families in the transcriptome data. The correlation coefficient of the number of expansions and contractions and ratio of contractions/expansions between genomic and transcriptomic data were analyzed for each pair using R (version 3.5.2; R Development Core Team, www.R-project.org; accessed on 22 October 2020).

Evolutionary Rate and Positive Selection Analyses
The number of rapid-evolution/positively selected genes detected would have been greatly reduced if only the homologous genes existing in all species or only the conservative blocks in multiple sequence alignments were selected. Therefore, we performed analysis across the genes orthologous in at least four species.
The codeml program in the PAML package [56] with the free ratio model (model = 1) was used to estimate the evolutionary rate following Yang et al. 2014 [58] and Wang et al. 2015 [59]. The lineage-specific mean ratios of non-synonymous (dN) to synonymous (dS) substitutions rates (ω = dN/dS) were estimated for each homologous gene, using the final phylogenetic tree as the guide tree. The values of dN, dS, and dN/dS were obtained for each branch. Assessments of the statistical significance of the differences in the dN/dS ratios along different lineages were conducted using the Wilcoxon rank sum test. To find genes that potentially experienced positive selection, the branch-site model (model = 2 and NSsites = 2) of the PAML package was used, with each branch specified as the foreground branch according to the following rigorous criteria: dN/dS ratio (ω) of the foreground branch greater than the background; p-value ≤ 0.05 in the likelihood ratio test [60]; or positively selected sites with a posterior probability greater than 0.95 [61]. The functions of genes with rapidly evolving rates and positive selection were estimated from GO and KEGG. We again compared the proportion of enriched genes from GO and KEGG between the foreground and background gene families. Significance was tested using Fisher's exact test.

Annotation of Chemosensory Genes
We explored the genetic basis for chemosensory variation among wasps [62]. The number of protein sequences for the following gene families was compared in 25 fig wasps and 7 non-fig wasp insect species: odorant binding proteins (OBPs), olfactory receptors (Ors), chemosensory proteins (CSPs), ionotropic receptors (Irs), and gustatory receptors (Grs). We searched for these families in the Pfam A database using the hmmscan command inHMM v3.3.2, and the results were filtered using a GA bitScore threshold with an e-value of 1e-5 and 25% HMM coverage.

Comparison of Transcriptome Sequencing and Assembly among 25 Fig Wasp Species
We sequenced transcriptomes of 25 fig wasp species comprising six representative genera of the family Agaonidae (Table 1)

Unigenes Annotated against Public Databases
Blasting of BLASTP or HMM against protein databases (Nr, Swiss-Prot, KOG, eggNOG, Pfam, KEEG, and GO) revealed that 7859 to 52,340 unigenes (median: 10,437) were successfully annotated for 25 fig wasps (Figure 1; Supplementary Materials, Table  S2), and the proportions of unigenes annotated ranged from 75.99% to 89.53% (median: 84. 70%). The proportions annotated against Nr and eggNOG were relatively high, with   Table S2), and the proportions of unigenes annotated ranged from 75.99% to 89.53% (median: 84. 70%). The proportions annotated against Nr and eggNOG were relatively high, with Nr from 71.45% to 85.89% (median: 78.97%) and eggNOG from 61.03% to 82.72% (median: 76.16%). The proportion annotated against the PFAM was the third highest (from 61.99% to 74.83%; median: 69.14%). The proportions annotated against other databases are given in Table S2. In total, 17,952 and 21,285 OCGs could be annotated by GO and KEGG, respectively, for downstream enrichment analysis.

Phylogenetic Relationships
The phylogenetic tree generated from 625 single-copy homologous genes was further integrated into one final tree with all branches supported by 100% bootstrap support (Figure 2b

Gene Family Expansion and Contraction at the Level of Expression
The patterns of expansion and contraction and the ratios of contraction/expansion of the gene family for five insect species, A. mellifera, C. floridanus, N. vitripennis, C. solmsi, and D. melanogaster, showed no consistent correlation (correlation coefficients (r) of 0, −0.34, and −0.55, respectively; Supplementary Materials, Excel S1). For the genome, the numbers of expanded gene families ranged from 226 to 1426, while for transcriptomes they ranged from 44 to 420; for genome-level family contractions, the range was from 1186 to 2808, while for transcriptomes this number ranged from 365 to 7067. For transcriptomes the number of contracted gene families for all five species was larger than

Gene Family Expansion and Contraction at the Level of Expression
The patterns of expansion and contraction and the ratios of contraction/expansion of the gene family for five insect species, A. mellifera, C. floridanus, N. vitripennis, C. solmsi, and D. melanogaster, showed no consistent correlation (correlation coefficients (r) of 0, −0.34, and −0.55, respectively; Supplementary Materials, Excel S1). For the genome, the numbers of expanded gene families ranged from 226 to 1426, while for transcriptomes they ranged from 44 to 420; for genome-level family contractions, the range was from 1186 to 2808, while for transcriptomes this number ranged from 365 to 7067. For transcriptomes the number of contracted gene families for all five species was larger than that of the expanded gene family, with the ratio of contraction/expansion of the gene family being between 6.84 and 27.5; this was also the case for genome-level data, which ranged from 1.55 to 12.42, except for one species, N. vitripennis, which has a ratio of 0.83.
Compared to the most closely related outgroup taxon, Nasonia vitripennis, there are 10 gene families that have undergone contraction in the family Agaonidae. These contracted gene families have functions related mainly to immune response and signal transduction ( Table 2   The numbers of expanded and contracted gene families varied greatly among fig wasp species. This was in contrast to their closest clade with 5-770 (mean ± SE: 184.2 ± 178.8) and 0-2532 (440.0 ± 841.0), respectively (Table 2; Figure 2a,b). While the gene families enriched in GO and KEGG were rarely shared among species, their functions were mainly related to signal transduction, immune system, drug resistance, endocrine, energy metabolism, digestive system, protein production, cytoplasmic translation, and regulation (Supplementary Materials, Excels S3 and S4). The numbers of contracted gene families in genus clades in the phylogenetic tree were greater than that the numbers of expanded genes. The numbers of contracted gene families in most species were lower, except for two related species of V. javana sp. 7 and sp. 2 and four taxa of Blastophaga (Figure 2a,b). For contracted gene families, the GO-and KEGG-enriched gene families in V. javana sp. 7 and sp. 2 were related to amino acid metabolism, signal transduction, energy metabolism or carbohydrate metabolism, and the nervous system, but the specific gene families and metabolic pathways (and, indeed, the proteins that they produce) were different (Supplementary Materials, Excels S3 and S4). Four enzymes or gene families related to protein synthesis (e.g., serine and threonine) enriched in KEGG were shared between B. sp.-F. abeli and B. sp.-F. pyriformis (Table 2; Supplementary Materials, Excel S4). Ribosome, valine, leucine, and isoleucine biosynthesis-related genes as well as neurodegenerative disease-related genes were shared between B. sp.-F. formosa and B. sp.-F. erecta var. beecheyana. No KEGG pathway was shared among the four taxa.

Contraction of Genes Involved in Chemosensory
It has been reported that some chemosensory gene families in fig wasps have experienced dramatic contractions in relation to other insects [6]. Therefore, we compared the numbers of genes in OBP, Or, CSP, Ir, and Gr families among fig wasps and other insect species (Table 3; Figure 3a

Rapidly Evolving Genes (REGs)
In the family Agaonidae, we detected 857 rapidly evolving genes (REGs) ( Table 2; Figure S1b). Of these REGs, 484 and 45 were enriched in GO and KEGG, respectively; their gene functions were related mainly to signal transduction, immune response, and antibacterial systems (Table 2; Supplementary Materials, Excels S5 and S6). For the Eu/Pl, Ceratosolen, Valisia, and Blastophaga genus clade, we detected higher numbers REGs, from 2073 to 2572 (2262.8 ± 222.5), than the numbers of positively selected genes, which ranged from 10 to 85 (62.8 ± 35.4). Low numbers of REGs were enriched in GO in the clades of Valisia (17), Eu/PL (14), Blastophaga (1), and Ceratosolen (4) (Supplementary Materials, Excel S5) while no REGs were enriched in KEGG. There were 13 REGs enriched in GO that were shared between Valisia and Eu/Pl, and these mainly related to mitochondrial function. The REGs enriched in GO in Blastophaga were related to the vesicle-tethering complex, while four of the Ceratosolen REGs encoded antibacterial peptides (Supplementary Materials, Excel S5).

Rapidly Evolving Genes (REGs)
In the family Agaonidae, we detected 857 rapidly evolving genes (REGs) ( Table 2; Figure S1b). Of these REGs, 484 and 45 were enriched in GO and KEGG, respectively; their gene functions were related mainly to signal transduction, immune response, and antibacterial systems (Table 2; Supplementary Materials, Excels S5 and S6). For the Eu/Pl, Ceratosolen, Valisia, and Blastophaga genus clade, we detected higher numbers REGs, from 2073 to 2572 (2262.8 ± 222.5), than the numbers of positively selected genes, which ranged from 10 to 85 (62.8 ± 35.4). Low numbers of REGs were enriched in GO in the clades of Valisia (17), Eu/PL (14), Blastophaga (1), and Ceratosolen (4) (Supplementary Materials, Excel S5) while no REGs were enriched in KEGG. There were 13 REGs enriched in GO that were shared between Valisia and Eu/Pl, and these mainly related to mitochondrial function. The REGs enriched in GO in Blastophaga were related to the vesicle-tethering complex, while four of the Ceratosolen REGs encoded antibacterial peptides (Supplementary Materials, Excel S5).

Positively Selected Genes (PSGs)
In the family Agaonidae, we detected 68 PSGs and none of these were GO-or KEGGenriched (Table 2; Figure S1b and Excels S7 and S8). These PSGs mainly coded for proteins related to signal transduction (CAMK1, FOXG, GRIN, IRAK4, PLK2, STAT5B, and KCNH8), immune response, antibacterial systems, genetic information processing (transport, translation, transcription, membrane trafficking, replication, and repair), development and regeneration, amino acid metabolism, and energy metabolism (Supplementary Materials, Excel S9). When making comparisons at the genus level, we found 75-110 (87.8 ± 15.4) genes under positive selection; of these genes, 0-13 were enriched in GO and 0-2 were enriched in KEGG ( Figure S1b and Excels S7 and S8). These genes mainly coded for proteins involved in energy metabolism, genetic information processing, and environmental adaptation. There was a single common KEGG, genetic information processing, shared among the three genera associated with dioecious hosts.
In total, the numbers of PSGs were lower than those of REGs for Agaonidae, genera, and species (Table 2; Figure S1a,b). In contrast to REGs, the numbers of PSGs were much higher among the five related Valisia species (128-207 (t = 8.773, p < 0.001)) and the five taxa of Blastophaga (227-288 (t = 15.227, p < 0.001)) compared to those of other species (22-96) ( Table 2; Figure S1). The number of GO-and KEGG-enriched PSGs in most species was 0 (Supplementary Materials, Excels S7 and S8). Among the five related Valisia species, only the PSGs of V. esquirolianae were enriched in GO and KEGG which were related to genetic information processing, amide synthesis, and peptide metabolism. There was one common PSG among the five related

Comparisons among Agaonidae
Our results demonstrate a dramatic contraction across a wide range of gene families in pollinating fig wasps compared to other Hymenoptera, offering support for results from previous studies [6,33,63,64]. It has been reported that in the genomes of Ceratosolen solmi, Eupristina verticillata, and Wiebesia pumila, many gene families related to chemosensory, detoxification, and innate immune response are reduced [6,33,63,64]. This finding is reflected in our transcriptome study: for the whole family of Agaonidae, the contracted genes families at the level of expression in newly emerged adult fig wasps were mainly enriched in relation to signal transduction and immune response in comparison to other parasitic wasps (Nasonia vitripennis and Copidosoma floridanum), ants, bees, and other Hymenoptera. Genes encoding for proteins used in environmental information processing, antibacterial defenses, amino acid and energy metabolism, cell growth and death, and the nervous system evolve quickly and are under positive selection. In the genome of Ceratosolen solmi, 13 genes were identified as rapidly evolving genes and they often act as transporters of signals or substances [6].
Figs provide a relatively safe nursery for developing larvae, shielding them from the external environment and many antagonists [6,10] . Fig wasps display morphological, behavioral and physiological adaptations related to life inside the syconium, where nutrients (starch and amino acid) are produced within endosperm tissue. It is important, however, to note the context of this study. We generated transcriptomes from adult females and did not explicitly account for age or external stimuli; as such, our main findings are relevant to newly emerged wasps only. The inconsistent pattern of gene family expansion and contraction between genomic data and our transcriptomic data also implies that gene expression is easily influenced by the environment, although the number of gene family contractions was typically larger at both levels.
Chemosensory genes play an important role for insects in foraging, ovipositioning, mating, avoiding natural enemies, and searching for hosts [65]. Our results showed that Ors, Grs, and Irs in fig wasp species are reduced to a greater extent in comparison to other insect species, with similar numbers as Ceratosolen solmi, Eupristina verticillate, and Wiebes pumila-Ors: 43-46, 48, and 67; Grs: 5, 4, and 4; Irs: 11-31, 30, and 16, respectively [6,33,63]. Genes in the Or family are related to the detection of volatile signals, while those in the Gr family generally contribute to the detection of soluble chemical substances. Groups of receptors are essential for survival and reproduction, e.g., for host location and avoidance of natural enemies and while searching for mates [66]. In addition to detecting volatile odor molecules, such as ammonia, phenylacetaldehyde, and hexanal, the genes in the Ir family are also involved in taste, temperature, and humidity sensing when used alongside the genes in the Or family [67].

Comparisons at Genus Level
At the genus level, the numbers of contracted gene families were far greater than those of expansions. Moreover, the more phylogenetically distant the genera under comparison were, the greater the number of contracted gene families detected. In contrast to Valisia (the pollinator of a dioecious species and a close relative in the current sampling regime), the monoecious Eu/Pl clade had 1182 contracted gene families mainly enriched in "environmental information processing". We speculate that this was due to the adaptation of pollinators to the monoecious reproduction of their hosts. There were more REGs and PSGs, but their functions were rarely enriched in GO and KEGG, which indicates that the adaptation of these functional genes at the level of the genus was random and scattered. In terms of PSGs, only one enriched pathway of genetic information processing (ko03010) was shared among the three dioecious genera; this is a candidate pathway for adaptation to the dioecious breeding system, but more dioecious genera are needed to test this hypothesis.

Comparisons among Species
Our phylogenetically structured sampling strategy also allowed us to assess differences among species within Agaonidae. The GO-and KEGG-enriched genes across gene families for species may have related to the enclosed nature of the syconia, and functions included signal transduction, immune response, drug resistance, endocrine regulation, energy metabolism, digestion, protein production, cytoplasmic translation, and regulation. However, the exact genes expressed were rarely the same across species, which was consistent with the highly species-specific nature of the interaction with their host. REG and PSG genes among species were mainly related to energy metabolism, drug resistance, environmental information processing, genetic information processing, P450 function, and carbohydrate metabolism. However, few genes could be enriched with GO and KEGG at the species level, a similar result to the genus-level comparisons.

Comparisons among Closely Related Species
Among closely related Vallisia species or taxa of Blastophaga, the numbers of contracted gene families varied greatly, which may have been due to the low expression of these genes, implying that fig wasps can quickly respond to host changes through gene expression. Future studies focusing on genomic sequencing would provide the context needed to confirm this. Moreover, the KEGG-enriched genes in contracted gene families were rarely shared across close relatives; this indicates differential expression and potentially subtle differences related to host identity. Fig wasps can likely respond quickly to environmental changes through gene expression. Meanwhile, the PSG numbers increased significantly in these comparisons, suggesting that these species/taxa were under selection.
The ANK was the only gene under positive selection that was shared among the five related Valisia species, and it was also one of the only rapidly evolving genes enriched in V. javana sp 1. The function of ANK relates to the IMD pathway and it can be activated by the infection of Gram-negative bacteria. Activation of ANK subsequently leads to the activation of NF-κB signal transduction pathways, which induce the synthesis of antibacterial peptide genes and other genes involved in the immune response (for example, see [69]). The hosts of the V. javana clade are two related species, F. hirta and F. triloba, in the same Ficus subsection as F. hirta. These figs have a broad geographical distribution, which implies that their syconia have similar bacterial infections countered by the same adaptation in five pollinator species. Further, the ANK gene is the only gene family found in all currently sequenced bracoviruses (BVs) and ichnoviruses (IVs; [70]) of parasitoid wasps [71], where they are likely used to manipulate host physiology and defense [72].

Conclusions
Our large-scale transcriptomic dataset was derived from 25 fig wasps, 6 other Hymenoptera species, and 1 Diptera species. First, we reported the genomic mechanisms underlying this highly species-specific interaction at the levels of family, genus, and species.