Kin recognition: Neurogenomic response to mate choice and sib mating avoidance in a parasitic wasp.

Sib mating increases homozygosity, which therefore increases the risk of inbreeding depression. Selective pressures have favoured the evolution of kin recognition and avoidance of sib mating in numerous species, including the parasitoid wasp Venturia canescens. We studied the female neurogenomic response associated with sib mating avoidance after females were exposed to courtship displays by i) unrelated males or ii) related males or iii) no courtship (controls). First, by comparing the transcriptional responses of females exposed to courtship displays to those exposed to controls, we saw a rapid and extensive transcriptional shift consistent with social environment. Second, by comparing the transcriptional responses of females exposed to courtship by related to those exposed to unrelated males, we characterized distinct and repeatable transcriptomic patterns that correlated with the relatedness of the courting male. Network analysis revealed 3 modules of specific 'sib-responsive' genes that were distinct from other 'courtship-responsive' modules. Therefore, specific neurogenomic states with characteristic brain transcriptomes associated with different behavioural responses affect sib mating avoidance behaviour.

Introduction

Mate choice often depends on individual attractiveness. In many species, females discriminate between competing males [1]. Female mate choice sometimes occurs after elaborate courtship displays, and males are selected according to their highly heritable qualities (‘good genes’ hypothesis) and/or their genetic compatibility with the female (‘compatible genes’ hypothesis) [2]. The ‘good genes’ hypothesis predicts that females favour reproduction with males carrying traits that are honest indicators of good genes or as a result of sensory bias [3, 4], hence obtaining genetic benefits for their offspring [5]. The ‘compatible genes’ hypothesis predicts that offspring fitness is correlated with genetic dissimilarity between mating partners, where each female prefers males that possess genes compatible with their own genotypes. Sib mating avoidance fits into this second hypothesis [6, 7]; in this case, kin recognition is crucial for enabling females to discriminate between related and unrelated mates. Sib mating avoidance has been extensively documented in numerous animal species (including species that are social and solitary) and reduces the risk of inbreeding depression [8].

In colonial marine invertebrates, contact between 2 conspecific colonies may trigger fusion if the 2 colonies are identical or related, or aggressive rejection if the 2 colonies are not related. Molecular bases for kin recognition involve several highly polymorphic loci that have been well characterized in 2 species: the cnidaria Hydractinia symbiolongicarpus (reviewed in [9]) and the chordata Botryllus schosseri (reviewed in [10]). In terrestrial animals, most of the presently identified mechanisms usually rely on individual-specific olfactory cues and sensory organs able to discriminate these cues from one another. In vertebrates, major histocompatibility complex genes confer individual olfactory identities [11-14]. These olfactory cues bind to receptors often located on the vomeronasal organ neurons and are used by females to evaluate mate relatedness [15-19]. In insects, numerous species discriminate kin by using odour differences inherent to divergence in the cuticular hydrocarbon composition [20-23]. Odour molecules bind to specific sensory receptors located on sensilla basiconica, where dozens of olfactory sensory neurons project to the antennal lobes [24-26], where olfactory second-order neurons project to the lateral horn and mushroom bodies [27, 28]. To the best of our knowledge, no genome-wide study has addressed how kin recognition signals are processed in the context of mate choice and sib mating avoidance. Transcriptomic studies provided recent insight into female mating decisions [29]. Coordinated changes in the expression of many genes in female brains, i.e., a neurogenomic response, have been identified following courtship in Poeciliidae fishes [30-32]. This response depends on male attractiveness and is in accordance with female preferences.

Hymenopterans (including bees, ants, wasps and sawflies) are relevant models to explore the biology of mate choice, especially regarding sib mating avoidance and its underlying biochemical mechanisms. These insects share an ancestral and unusual sex determination system called single-locus complementary sex determination (sl-CSD), where sex is determined by heterozygosity at the sl-SCD locus. Heterozygous individuals develop as females, hemizygous individuals develop as males, but homozygous individuals develop into diploid unviable or sterile males [33-36]. When sib mating occurs, the risk of genetic incompatibility between partners sharing the same allele at the sl-CSD locus is one in two; in that case, half of the progeny will be homozygous diploid sterile males (with notable exceptions [37, 38]). Therefore, the sl-CSD has an additional immediate effect on inbreeding. Consequently, Hymenoptera are much more exposed to inbreeding depression than any other diploid species [39]. In this context, congruent with the compatible gene hypothesis, one expects that selection will favour females that prefer males with a distinct allele at the sl-CSD locus, thus excluding sibs. In the parasitoid wasp Venturia canescens, which has sl-CSD [40], females only mate once [41], making mate choice particularly decisive. Indeed, in this species, females are able to discriminate kin and non-kin during male courtship based on olfactory-mediated cues [42]. In laboratory conditions, the proportion of successful mate when a single female was in presence of 2 unrelated males was higher (79% of success) than in presence of related males (i.e. with 2 brothers) (41% of success) [42]. Furthermore, the latency to mate increased when the courting male was related [42]. Experiments involving choice between 2 males, one brother and one unrelated present in a same area have shown that female mates indifferently with the 2 males [42]. Hence, the two distinct social contexts (presence of a conspecific with and without relation) provoke two contrasting behavioural responses that correspond to the definition of behavioural states, i.e., the performance of a distinct and quantifiable behaviour for a measurable period of time [43].

We tested whether these contrasting behavioural responses to relatedness were correlated to transcriptomic changes in V. canescens females. We designed a behavioural experiment where females were exposed for 10 minutes to an unrelated male or a related male (a brother) or were without any social contact. Then, we used RNA sequencing to compare gene expression profiles from the entire heads of these females. This timing coincided with a period of active courtship, where female wasps perceived and evaluated males, experienced change in receptivity and decided whether or not to mate [44]. Previous studies demonstrated that this temporal frame is adequate to detect early transcriptional changes following mate exposure in different species [30, 45, 46]. We focused on the whole head, where sensory signals are processed and integrated to mediate complex decision making, including mating decisions.

We hypothesized that information processing during male courtship display and mate choice were related to gene expression changes. First, we examined if courtship perception was mediated by changes in gene expression. In this case, we expected the transcriptomic profiles of all courted females would be analogous, whatever the relatedness of the courting male. Then, we assessed whether the differences in female receptivity were associated with gene expression changes. Under this hypothesis, the expectation was that non-receptive females, i.e., solitary females and females courted by related males, would have similar transcriptomic profiles.

Materials and methods

Biological model

Venturia canescens G. (Hymenoptera: Ichneumonidae) is a solitary parasitic wasp; a maximum of one adult emerges from one host, regardless of the number of parasitoid eggs initially laid in the host. In the wild, V. canescens females parasitize pyralid moth larvae feeding on dried fruits, such as figs, carobs, almonds, dates or loquats [47]. Our knowledge of how mating partners encounter each other in the field is largely incomplete as a consequence of the small size of the species, which renders observations difficult. Virgin V. canescens females emit chemicals that, in combination with host kairomones, attract males [42]. In turn, males do not attract virgin females at a distance [41]. V. canescens females are monandrous and only mate once in their lifetime [41]. Conversely, males can mate more than once. These observations led us to consider that females choose their mate, and we thus focused our study on female transcriptomic responses to courtship. Moreover, females recognize sibs on the basis of a chemical signature emitted by males and can avoid mating with their brothers through kin recognition [42]. This kin recognition ability has also been robustly established in this species in the context of host choice, since females prefer to lay their eggs in hosts already parasitized by conspecifics rather than in hosts already parasitized by themselves [48]. In addition, males also have the ability to distinguish between non-sib and sib females by using chemical markers emitted by females [49].

Behavioural manipulations

Wild parasitic wasps were collected via a large field sampling during the summer of 2014 in orchards in southern France (Valence, N 44°58’21” E 4°55’39”). Wasps were maintained on the host Ephestia kuehniella (Lepidoptera: Pyralidae) in a climatic chamber (25°C, 60% humidity, DL 12:12) and fed with a 50/50 honey/water solution. Families were produced as described by Metzger et al. (2010). A total of 90 newly emerged virgin females were individually placed in a box and randomly submitted to one of the three following conditions for a maximum duration of 10 minutes: 1) a compatible male (i.e., unrelated male), 2) an incompatible male (i.e., a brother), or 3) isolated (controls). The experimental design has been conceived in order to minimize the influence of circadian rhythm and genetic background on the results. Newly emerged females were isolated every morning and numbered. A random draw was then made to establish the order of passage of the females and the treatment assigned. Behavioural experiments took place in the afternoon between 1:00 and 4:00 p.m. for all females that were captured in the morning in the order established by the random draw. A maximum of 2 females from each family have been kept daily, and randomized in one condition. For each given condition, the ten females belonged to different families to avoid an effect of genetic homogeneity on transcriptomic results. New families have been produced for every biological replicate. The 10 minutes period coincides with an active male courtship behaviour [44]. For each female under conditions 1 and 2, one sequence of active courtship parade was observed within the 10 minutes following the introduction of the male to the box for all pairs. Whenever this behaviour sequence was not observed, the pair was eliminated from the study. After this time, or within 3 seconds of mounting, the female and male were separated. For each condition, females were caught in a small tube (3 x 7 cm) and were immediately anaesthetized with CO2 and decapitated with a scalpel. Antennae were removed from the head to focus on the cerebral transcriptome of the decision-making centre rather than the sensory centre. The head is a heterogeneous tissue divided into many functionally distinct regions that could affect inferences about the functional significance of gene expression patterns [50]. However, the study of more specific brain regions remains challenging given the small size of the wasps. Collected heads were immediately flash-frozen and stored at -80°C until RNA extraction. A total of 10 individual heads were collected per experimental condition and pooled to collect higher RNA amounts and to average the expression state across individuals, thus mitigating individual variability during transcriptome comparison [51]. The experiment was sequentially repeated 3 times to obtain biological replicates.

RNA extraction, library preparation and sequencing

Total RNA extraction was performed in one batch on the 9 biological samples collected (3 experimental conditions x 3 biological replicates). Mechanical cell lysis was performed with metallic beads added to frozen microtubes and shaken with a TissueLyser (Qiagen, 25 hertz 45 seconds). RNA was extracted using Qiazol and an RNeasy mini kit according to the manufacturer’s protocol (Qiagen), with optional DNase treatment (Thermo Fisher Scientific). The quality and quantity of total RNAs were assessed on a Bioanalyzer (Agilent) and a Qubit fluorometer (Life Technologies). Polyadenylated RNAs were enriched from 1 μg of high-quality total RNAs with oligo-dT magnetic beads, fragmented and converted to cDNAs (Illumina TruSeq Stranded mRNA Library Prep kit). After adding an A to the 3’ end of each cDNA, adapters were ligated, and fragments were amplified by PCR to generate DNA colonies. Each library was labelled, multiplexed and pooled for sequencing on a HiSeq 2500 Illumina sequencer (Fasteris, Switzerland), with a paired-end protocol (2x125 bp). For each of the 9 biological samples, a second library was prepared independently, sequenced and considered a technical replicate. Overall, a total of 18 libraries were sequenced (3 experimental conditions x 3 biological replicates x 2 technical replicates). A minimum of 18 million paired-reads were obtained per library, representing a final dataset of more than 535 million paired-reads.

RNA-seq quality control, mapping, transcriptome assembly and annotation

The V. canescens genome has been sequenced and is currently available (http://bipaa.genouest.org/sp/venturia_canescens/, v.1.0) but thus far, RNA-seq has not been performed on heads [52, 53]. To include all loci detected in the head transcriptome that might be missing from the current annotation, we constructed a de novo transcriptome assembly using the genome reference with the TopHat-Cufflinks pipeline [54] (v2.2.1). The global quality of sequences was assessed with FastQC-0.10.1. Low-quality bases at the ends of reads were trimmed using Trimmomatic-0.36. Adaptor-containing reads, low quality reads (scores phred <30) and N-containing reads were filtered out. Curated reads were then processed through the TopHat-Cufflinks pipeline. First, paired-reads of each sample were aligned to the V. canescens genome. TopHat used the Bowtie algorithm to align reads to the genome. Unmapped reads were cut into segments that aligned far apart from one another (between 100 bp to several hundred kb) to predict potential intron/exon structures. Next, an index of hypothetical splice sites was built without any prior information. Each site had to be confirmed by several read segments consistently showing the same alignment pattern. TopHat was used with each of the 18 libraries as input and with the following parameters: stranded libraries (first strand), 5 mismatches/indels allowed, and report the best alignment possible for each read. After mapping, resulting alignment files were provided to Cufflinks, which generated a transcriptome assembly for each sample annotated into genes, transcripts and isoforms. Then, the 18 transcriptomes were merged into one master transcriptome with Cuffmerge. This final transcriptome assembly contained a total of 16,752 genes and provided a uniform basis for calculating gene expression in each condition. The transcriptome was then annotated to gain insight into the functions of the transcripts and the proteins they encode. Sequence similarity was searched for each of the 16,752 predicted genes by comparing the six-frame translation putative products of the nucleotide sequences using BLASTX (v2.2.29+) against the NCBI non-redundant protein database. Transcriptome analysis was completed using gene ontology (GO) annotation, which associated genes with functions in 3 categories: molecular function (molecular activities of gene products), cellular component (where gene products are active) and biological process (pathways and larger processes made up of the activities of multiple gene products). This classification enabled functional interpretation of a large group of genes via enrichment analysis.

Sample-based clustering and PCA

To obtain an initial overview of gene expression patterns across samples, multivariate analyses were performed on the 18 transcriptomic libraries, including 3 biological replicates for each of the 3 conditions representing a total of 9 samples, each one represented by 2 technical replicates. Those 18 transcriptomes were represented by isolated females (6), females courted by unrelated males (6) or females courted by brothers (6). A raw count table was obtained by using HTSeq [55] (v0.5.4p1). The gene model produced with Cuffmerge was combined with the 18 mapping files previously obtained with TopHat for each gene in each library. This final dataset was exported to R (v3.4.3) [56] for downstream statistical analysis with DESeq2 [57] (v1.30.0). Counts were normalized with the variance stabilizing transformation method (VST), which produced transformed data on the log2 scale and normalized data with respect to library size. The overall variation of expression levels among samples was evaluated with a two-dimensional PCA and with hierarchical clustering, both based on the expression of the 500 genes with the highest variance across libraries. A hierarchical clustering dendrogram was produced based on the sample-to-sample Euclidian distance matrix to obtain an overview of the similarities and dissimilarities between samples. Uncertainty in hierarchical clustering was assessed with Pvclust [58] (v2.0–0) using multiscale bootstrap resampling with an approximately unbiased P-value to measure statistical support for each cluster (1,000,000 bootstraps; average agglomerative method; correlation method distance). Two technical replicates from one biological sample of females courted by brothers were obvious outliers in the PCA and hierarchical clustering, most likely due to a problem that occurred during sample conditioning and were thus excluded from the analysis (S1 Fig).

Differential expression, co-expression network and functional enrichment analysis

All pairs of technical replicates were merged before proceeding to differential expression analysis and network analysis as recommended by Love et al. (2014), keeping a total of 8 samples (3 isolated females, 3 females courted by unrelated males and 2 females courted by brothers). To identify genes with different expression patterns across conditions, we performed pairwise comparisons between 1) females courted by related and unrelated males to controls (i.e., isolated females); and 2) females courted by unrelated males to females courted by related males. Differential expression was tested by using negative binomial generalized linear models implemented in the program DESeq2. We tested for differential expression of all transcripts with an average level of expression superior to 10 reads per library (n = 14,034). After normalization, within-group variability (i.e., the variability between biological replicates) was modelled for each gene by the dispersion parameter, which describes the variance of counts by sharing information across genes, assuming that genes of similar average expression strength have a similar dispersion. Among the 14,034 genes tested, only 7 presented an outlier status (i.e., an inconsistent expression pattern across the dataset) and were thus excluded from the test, keeping 14,027 transcript in the final reference transcriptome. For the subset of genes that passed the filtering test, a Wald test P-value was calculated and finally adjusted for multiple testing [59]. A gene was considered differentially expressed when the false discovery rate (FDR)-adjusted P-value was less than 0.01. We did not apply any log fold change threshold.

We performed a weighted gene co-expression network analysis (WGCNA [60], v1.63) to identify subgroups of genes that shared common expression patterns across the experimental conditions and potentially drove the differences in mate choice. WGCNA is a data reduction technique that regroups genes with similar expression patterns into modules of co-expressed genes and tests the correlations between modules and traits. First, log2-transformed and VST-normalized counts were used to construct a gene co-expression network with the blockwiseModules function. A correlation matrix was computed for the genes, and the correlations were weighted using a power function ß. Then, genes sharing similar patterns of variation across conditions were regrouped using hierarchical clustering and a dynamic tree-cutting algorithm to define modules of co-expressed genes. For our analysis, the parameters used were as follows: maximum block size = 15,000 genes; power (ß) = 10, minimum module size = 30 genes. The remaining parameters were kept at the default settings. A colour name was assigned to each module, and biologically interesting modules were identified by correlating a summary profile for each module to external experimental conditions; P-values <0.05 were considered significant and numbered from 1 to 11. Finally, the potential ‘hub genes’ in every significant module were identified. So-called hub genes may influence the expression of other genes in their module and may be causal factors for a trait of interest. Such hub genes are potentially biologically relevant by driving phenotypic variations [61, 62]. The identification of hub genes relies on both connectivity with other genes from the module and the correlation to the trait. Accordingly, genes were ranked according to their module membership values in each module. The top 5 hub genes of every module were annotated, and their expression pattern, which was representative of the module they belonged to, was detailed.

Functional characterization of gene sets (i.e., genes with differential expression, or modules of co-expressed genes) was analysed using enrichment analysis and gene ontology annotations. An enrichment test was performed on test sets compared to the full transcriptome with Blast2GO [63] (v5.0). The proportion of genes associated with GO terms was compared between the test set and the transcriptome (14,027 transcripts) with a unilateral Fisher’s exact test (one-sided), P-values < 0.01 were considered significantly enriched.

Results

The transcriptomic response to courtship is consistent with the mate preference of females

We first determined whether the male courtship display provoked a transcriptomic response in females. To do so, we built the first V. canescens head transcriptome by sequencing more than 500 million reads using a genome guide assembly (18 to 55 million paired-reads from 18 libraries, S1 Table) and used it as a reference to compare the different head transcriptomes of females. After quality filtering a mean of 98.4% of paired-reads were kept, on which 70.8% were successfully mapped to the genome (representing a mean of 20 millions per sample, for a total of 363 million paired-reads, S1 Table). Such values corresponded to the high-quality standards observed in other Hymenopteran species with an annotated genome [64]. The transcriptome constructed with these sequences encompassed 16,752 genes. Overall, 76% of predicted genes get a blast hit (12,740) while 4,012 sequences get no hit. Among genes matching with blast, 89.4% presented their best hit with an insect sequence, of which 84.4% match more specifically to a hymenopteran insect sequences. Finally, most of the genes were uncharacterized, since only 33.4% of the predicted genes (5,589) were associated with at least one Gene Ontology (GO) functional annotation.

Based on this unique reference, head transcriptomes from females 1) courted by an unrelated male, 2) courted by a related male, and 3) isolated (controls) were first analysed without a priori knowledge. Principal component analysis (PCA, Fig 1A) revealed that transcriptomes shared higher similarity within a given experimental condition than between different experimental conditions. PCA defined three consistent clusters according to the social environment proposed: the control group where females were kept in isolation, the group of females courted by unrelated males, and the group of females courted by related males. Principal component 1 separated females courted by related males from the two other conditions and explained the largest fraction of variance in gene expression (46%). Principal component 2 separated females courted by unrelated males from the two other conditions and accounted for 17% of the variance in gene expression (Fig 1A). These results were supported by hierarchical clustering analysis of sample-to-sample distance, which established that samples were clustered by the experimental conditions that they were exposed to (Fig 1B). Bootstrap resampling provided strong statistical support for this result: the control group (isolated females) formed a cluster with 100% support (Fig 1B), while the group of females courted by unrelated males constituted a cluster with 100% support. Finally, females courted by related males constituted a group with 95% support (Fig 1B). Together, these results showed that i) male courtship had an effect on the female transcriptome, and ii) being courted by a related or an unrelated male provoked two distinct transcriptomic responses.

Fig 1

Multivariate analysis of the 16 RNA-seq libraries based on the gene expression profiles of the 500 genes with the highest variance across samples.

The samples are clustered by the experimental conditions that females were exposed to. Experimental conditions are represented by the three different colours (grey: control; orange: females courted by related males; blue: courted by unrelated males), shapes indicate the three biological replicates, and the filling of the shape shows the 2 technical replicates. (A) Two-dimensional principal component analysis (PCA) plot. (B) Hierarchical clustering dendrogram based on sample-to-sample distances. Statistical support is indicated by an approximately unbiased P-value with one million multiscale bootstrap resamplings (all bootstrap values were > 50%, those < 80% were not shown for clarity).

Functional characterization of differentially expressed genes following courtship

The transcriptomes of females courted by related and unrelated males were compared to those of isolated females to identify the female neurogenomic response to courtship. We identified 1,001 differentially expressed genes (DEGs), representing 7.1% of the 14,027 gene sets tested (Fig 2A, S2 Table). Among the 1,001 DEGs, 463 had higher expression in isolated females (3.3% of total transcriptome), gene ontology enrichment analysis reveals that this set of gene was enriched in DNA-binding Transcription Factor Activity (full list in S2 Table). In contrast, 538 DEGs were overexpressed in courted females (3.8% of total transcriptome), gene ontology enrichment analysis reveals that this set of gene was enriched in Reproductive Behaviour (full list in S2 Table).

Fig 2

Transcriptomic response to courtship: 7.1% of the total transcriptome (grey) was regulated in response to courtship (1,001 DEGs).

(A) Comparison of females courted by related and unrelated males to controls revealed 538 genes were overexpressed in courted females (negative fold change values) and 463 genes were overexpressed in isolated females (positive fold change values). (B) Boxplot showing significant expression changes following courtship in highlighted genes mentioned in the text. The X axes indicate gene names, the Y axes show normalized counts after log transformation, and the boxplot whiskers show the range of read counts between biological replicates. PDF, pigment-dispersing factor; DAT, Dopamine transporter; 5-HT1A, Serotonin transporter 1A.

Among the DEGs, we noticed three neuropeptides associated with female receptivity, namely, Capability, PDF and SIFamide, were all downregulated in courted females (Fig 2). The orthologue to muscleblind, required for normal regulation of female sexual receptivity [65], was also detected among the genes overexpressed in courted females. Notably, we also reported two genes associated with dopamine, ebony (downregulated in courted females), DAT (dopamine transporter, overexpressed in courted females), and one gene associated with serotonin (5-HT1A, serotonin receptor 1A, overexpressed in courted females) (Fig 2).

Relatedness of the courting male influences the female head transcriptome

We compared females courted by unrelated males to those courted by related males to determine whether the relatedness to the courting male had an impact on the female transcriptomic response. By comparing these two groups, we identified 831 DEGs representing 5.9% of the tested genes (Fig 3A, S3 Table). Among these DEGs, 481 genes presented higher expression in females courted by related males (3.4%), gene ontology enrichment analysis reveals that this set of gene was enriched including ATP metabolism and Ribosome (full list in S3 Table). Moreover, 350 DEGs were overexpressed in females courted by unrelated males (2.5%). Enrichment analysis revealed 22 associated GO terms including Reproductive Behaviour and Male Mating Behaviour (full list in S3 Table). These 2 GO terms both refer to the same 2 genes, which belong to the yellow-major royal jelly protein family (Fig 3B, S2 Fig), an insect gene family notably associated with behaviour and caste specification [66]. Interestingly, we also noticed the regulation of the transcription factor PAX6, which is required for normal brain structure and function, notably locomotory behaviour [67] (Fig 3B).

Fig 3

Influence of relatedness of courting males on females’ transcriptome expression: 5.9% of the total transcriptome (grey) was regulated according to the relatedness of the courting male.

(A) Comparison of females courted by related males to females courted by unrelated males showed 481 genes overexpressed in females courted by related males (positive fold change values) and 350 genes overexpressed in females courted by unrelated males (negative fold change values). (B) Boxplot showing significant expression changes according to the relatedness of the courting male for highlighted genes mentioned in the text. The X axes indicate gene names, Y axes show normalized counts after log transformation, and boxplot whiskers show the range of read counts between biological replicates. Yellow—mrjp: yellow major royal jelly proteins; PAX6: paired box 6.

Co-expression network and characterization of ‘courtship-responsive’ modules and ‘sib-responsive’ modules

To better characterize the variations in gene expression according to the experimental conditions experienced by females, we applied a co-expression network analysis using WGCNA on the 14,027 genes that passed through the expression filter. The gene co-expression network grouped genes that shared a similar expression pattern across different experimental conditions into modules.

Overall, the 14,027 genes were organized into 50 modules of highly correlated genes symbolized by a colour, with sizes varying from 32 to 2,675 genes (Fig 4, Table 1). Among the 50 modules defined by the cluster analysis, 11 modules (numbered from 1 to 11) had significant correlations with at least one of the experimental conditions experienced by the females (Table 1, Fig 4B, S3 Fig); these were considered biologically relevant and were further analysed.

Fig 4

Gene co-expression network of the full transcriptome (14,027 genes) in response to courtship and according to relatedness with the courting male.

(A) Clustering dendrogram of 14,027 genes (top side), with assignation to 50 modules of co-expressed genes represented by colors (bottom side). B) Summary plot of modules dendrogram and relationship with social environment experienced. The right panel shows the dendrogram of modules from gene network. The heatmap (left panel) illustrates pairwise correlation between modules of co-expressed genes and social environment: red denotes high positive correlation, white indicates no correlation, while blue indicates high negative correlation. The 11 modules significantly correlated with at least one of the social environment experimented by females were numbered from 1 to 11. Social environments were labelled with symbols: §; control (isolated females), $; females courted by related males, and &; females courted by unrelated males.

Table 1

Significant correlations between 11 modules and courtship experienced by females.

Module name	Module size	representative enriched GO terms	Controls (isolated females)	Females courted by related males	Females courted by unrelated males
1 (purple)	227	molybdenum ion binding	NS	NS	-0.79 (0.02)
2 (green)	935	nucleotide metabolic process	NS	NS	-0.73 (0.04)
3 (pale turquoise)	77	mitochondrion, respiration	NS	NS	-0.8 (0.02)
4 (blue)	2108	ribosome, sensory perception	NS	0.94 (5e-04)	NS
5 (red)	643	Methylation	NS	-0.76 (0.03)	NS
6 (dark magenta)	73	GDP-mannose metabolic process	NS	-0.74 (0.04)	NS
7 (turquoise)	2675	mating behaviour	0.91 (0.002)	NS	NS
8 (dark green)	105	cell surface receptor signaling pathway	-0.77 (0.03)	NS	NS
9 (floral white)	38	oxidoreductase activity	NS	-0.8 (0.02)	0.82 (0.01)
10 (brown)	1989	response to stimulus, protein kinase activity	-0.73 (0.04)	NS	0.81 (0.01)
11 (yellow)	1022	carbohydrate catabolic process, metal ion binding	-0.8 (0.02)	0.79 (0.02)	NS

NS: not significant

Correlation coefficients and P-values in parentheses are indicated when significant (P-value <0.05). The results of enrichment analyses are indicated in the GO terms column (full lists in S4 Table).

First, three modules were associated with courtship display by an unrelated male (1, 2 and 3), which grouped 1,239 genes. Then, three modules were associated with sib-responsive genes (4, 5 and 6; S3 Fig; S4 Table), which grouped 2,824 genes responding only in the presence of related males. Next, two modules were associated with response to courtship, whatever the degree of relatedness of the courting males (modules 7 and 8; S3 Fig; S4 Table), which grouped 2,780 genes. Finally, three modules were associated with a gene expression pattern peculiar to each of the 3 social environments (9, 10 and 11; S3 Fig; S4 Table).

Discussion

In the current study, we characterized the female neurogenomic response associated with sib mating avoidance by identifying remarkable differences in the head transcriptome triggered by courtship display that differed according to the relatedness of the courting male. In V. canescens females, mate relatedness influences female sexual receptivity and is estimated during male courtship displays through chemical cues [42]. Unrelated males induce sexual receptivity in females, whereas related males induce weak sexual receptivity. Hence, sib mating avoidance can be considered a behavioural state, similar to other transitory behaviours such as aggressiveness [68] or singing [69]. We showed that sib mating avoidance is associated with distinct and repeatable cerebral transcriptomic patterns involving a significant part of the transcriptome (>5%). Such results fit the definition of neurogenomic state, i.e., a distinct and repeatable pattern of gene expression in the brain revealed by contrasting brain transcriptomes of individuals across different behavioural states. Despite the quite low number of biological replicates, the highly contrasted transcriptomes observed in the different social contexts suggest that sib mating avoidance behaviour could be considered a neurogenomic state. This research paves the way for further study on neurogenomic effects of sib mating avoidance in many species where such behaviours have been described and, thus, may contribute to the understanding of the molecular mechanisms underlying the evolution of avoiding consanguinity.

We measured major transcriptomic modifications occurring in the female head following a courtship display. Females exhibited transcriptomic pattern changes following an encounter with a partner, regardless of the relatedness of the courting male. This suggests that transcriptomic shifts immediately arise following the male display, even before the eventual copulation. These transcriptomic changes can be triggered by the presence of a male, mate evaluation, or from a social encounter. With the current experimental protocol, it is not possible to distinguish the origin of these changes. The addition of two other controls such as females in the presence of i) an unrelated female or ii) a related female would allow further specification of the impact of social environment on neurogenomic responses. Overall, 7.1% of the total transcriptome was differentially expressed at most within ten minutes after the courtship started (1,001 DEGs). Such a neurogenomic response on the time scale of minutes following environmental change is mediated by immediate early genes [70]. For such genes, near-instantaneous transcription is allowed by the presence of RNA polymerase II, which stalls in the promoter regions of these genes [71].

A recent transcriptomic study conducted on female mate preference in the guppy (Poecilia reticulata) showed that the presence of a potential mate induced changes in the brain transcriptome after only 10 minutes of exposure [30]. In insects, few studies have attempted to identify transcriptional changes associated with courtship displays; to the best of our knowledge, all have focused on Drosophila [45, 46]. Immonen and Ritchie submitted D. melanogaster females to a courtship song diffused by a speaker for 15 minutes and found only 41 DEGs (0.3% of the transcriptome) between courted and control females (in the presence of a male that was unable to perform a courtship display). Veltsos et al. (2017) submitted D. pseudoobscura females to male courtship and examined head transcriptomes immediately after mounting. They identified only 16 DEGs between courted and control females (virgin isolated females) (0.1% of the transcriptome). Together, these studies suggest that the female neurogenomic response to courtship only affects a very small set (>1%) of genes in Drosophila. Interestingly, in contrast with those studies, we found that the neurogenomic response to courtship was much greater in the wasp V. canescens. The differences may rely on the contrasted mating systems in these species. Despite the cost of reproduction, females mate multiple times in the majority of animal species in the wild including D. melanogaster [72] and D. pseudoobscura [73], most often with different males [74-76]. The benefits of polyandry for females include an adequate sperm supply [77], an increase in sperm competition [74] and a reduction in the cost of inbreeding [7]. In species with sl-CSD sex determination, mate choice is particularly determinant for female fitness, given the risk of genetic incompatibility. Furthermore, in monandrous species, such as V. canescens and a majority of parasitoid wasps (80%) [78], all progeny will have the same genitor. Thus, it is likely that selective pressures regarding mate choice should be stronger compared to those of polyandrous females. We suggest that the mating system might be an important determinant influencing the extent of neurogenomic response to courtship.

From a functional point of view, transcription factor activity and reproductive behaviour were some of the functions regulated following courtship. Indeed, we identified numerous transcription factors regulated following courtship, consistent with the large number of DEGs observed. In particular, we emphasized the transcription factors orthologous to Thd1 and muscleblind and the kinase Pink1. All of these genes were also differentially expressed following courtship in D. melanogaster [45]. These candidates may have a conserved regulation pattern following courtship and might be associated with the response to courtship in insects. Our results also indicated the regulation of genes related to neurotransmitters, such as dopamine. The dopamine transporter DAT is overexpressed in courted females and the ebony gene involved in dopamine catabolism is downregulated in courted females, which is compatible with an increase in dopamine concentration following courtship. Dopamine is notably implicated in the control of motivation, movement and memory in the fruit fly D. melanogaster [79], and increased dopamine levels result in increased responsiveness to courtship cues [80]. Furthermore, we detected three neuropeptides downregulated in females courted by an unrelated male: Capability, PDF and SIFamide. These neuropeptides are involved in female sexual behaviour in the fruit fly, since reduction or absence of SIFamide makes females extremely receptive [81], while PDF-mutant females show an increased frequency of re-mating compared to wild-type females [82]. These three neuropeptides are candidates for involvement in the modulation of female receptivity. Together, these results demonstrate that the neurogenomic response to courtship in V. canescens involves neurotransmitters and neuropeptides. These genes are prime targets for further functional analyses.

Although courtship accounted for the majority of the detected DEGs, our experimental design nonetheless highlighted the major influence of relatedness between partners on female response to courtship. In total, 9.1% of the transcriptome was differentially expressed when comparing females courted by related males to those courted by unrelated males. In addition to the 2,780 courtship-responding genes that exhibited the same expression pattern whatever the relatedness of the courting male (modules 7 and 8), the network analysis highlighted 3 modules of sib-responsive genes (modules 4, 5 and 6; 2,824 genes) and 3 modules of genes regulated only in females courted by unrelated males (modules 1, 2 and 3; 1,239 genes) that could be related to female receptivity. Our results clearly showed that the relatedness of the courting male had a major influence on the female response to courtship and that sib mating avoidance behaviour observed in this species is correlated to complex and massive changes in gene expression.

Concerning the functions associated with genes that vary according to the relatedness of the courting male, the GO terms Reproductive Behaviour and Male Mating Behaviour are notable. All genes underlying both terms are homologous to the yellow—major royal jelly proteins (yellow-mrjp). The yellow gene is unique to insects and some bacteria, while the mrjp gene family is restricted to Hymenopteran species and evolved from recent duplications of the yellow gene [83]. In the honey bee, the mrjp gene family is involved in both group social behaviour (royal jelly is constituted with 90% of MRJP proteins), and in individual sexual behaviour, with sex- and caste-specific gene expression patterns [81]. Yellow-mrjp functions in parasitic wasps are unknown, even though the largest expansion of the family described so far came from the Nasonia parasitic wasp genome, where genes are expressed broadly in different tissues and life stages [84]. In this study, we identified 6 members of the yellow-mrjp family that were DEGs following courtship, of which one presented differential expression according to the relatedness of the courting male. Further functional characterization will be required to test whether the yellow-mrjp gene family is involved in sib mating avoidance and female receptivity. Among the regulated genes, we also highlighted transcription factors such as PAX6 that could drive the transcriptomic changes accompanying female mate choice. Very few studies have explored the molecular pathways underlying kin recognition. In the amphibian Xenopus laevis, tadpoles exhibited plasticity in social preferences according to exposure to kinship odourants [85]. Sustained kin odourants exposure during development drives changes in neurotransmitter expression from GABA to dopamine neurons, which are stimulated from an increase in the expression of the transcription factor PAX6 and accompanied by a behavioural preference for kin odourants [85]. Here, we observed that the V. canescens PAX6 orthologue is downregulated in females courted by brothers. Further studies to characterize PAX6 function in V. canescens, particularly in the context of kin recognition, could test whether common molecular pathways could be elicited for kin recognition in distant taxa such as amphibians and insects.

We had formulated two non-mutually exclusive hypotheses. First, the perception of courtship was mediated by a change in gene expression, that would result in similar expression patterns in all females regardless of their relationship to the courting male. We identified such patterns for 2,780 genes (modules 7 and 8). Second, changes in female receptivity could result in changes in transcriptomic profiles. In this case, similar expression patterns would be expected for isolated females and females courted by their brothers. We have identified such expression patterns for 1,239 genes (modules 1, 2 and 3). Thus our results suggest that both courtship perception and changes in female receptivity induce a different neurogenomic response. In addition to sib mating avoidance, V. canescens females express kin recognition in the context of host choice, since females prefer to lay their eggs in a host parasitized by others than by a relative [48]. Neurogenomic analysis of responses to the presence of a relative in distinct ecological contexts would determine whether there are molecular bases associated with kin recognition.

Exclusion of one biological replicate from the analysis due to outlier status.

Multivariate analysis based on the expression profiles of the 500 genes with the highest variance across all samples showed that one biological replicate corresponding to females courted by related males (2 empty circles representing 2 technical replicates) is far from the other points (full circles) A) in the plan defined by the two first axes of the principal component analysis and B) in the sample hierarchical clustering dendrogram. Significant statistical support of outlier status is indicated by an approximately unbiased p-value with one million multiscale bootstrap replicates.

(TIF)

Click here for additional data file.

Boxplot showing expression patterns of the 6 members of the yellow–mrjp family were regulated following courting displays and/or according the relatedness between the female and the courting male.

Boxplot colours indicate biological condition: grey, isolated females; blue, females courted by unrelated males; orange, females courted by related males. The Y axes show the normalized counts after log transformation and VST normalization, and the boxplot whiskers show the range of reads between biological replicates. *, P-adj (FDR) < 0.01; NS, not significant.

(TIF)

Click here for additional data file.

Boxplots showing the expression patterns of the top 5 hub genes for each of the 11 modules varied according to experimental conditions.

The Y axes show the normalized counts after log transformation and VST normalization, and the boxplot whiskers show the range of reads between biological replicates.

(TIF)

Click here for additional data file.

Summary statistics of the sequencing runs and transcriptome assembly.

(XLSX)

Click here for additional data file.

Control females versus females courted by both related and unrelated males: A) list of 1,001 DEGs; B) 7 GO terms enriched in control females; C) 36 GO terms enriched in courted females.

(XLSX)

Click here for additional data file.

Females courted by related males versus females courted by unrelated males: A) list of 831 DEGs; B) 22 GO terms enriched in females courted by unrelated males; C) 66 GO terms enriched in females courted by related males.

(XLSX)

Click here for additional data file.

Enriched GO terms associated with the 11 significant modules of co-expressed genes.

(XLSX)

Click here for additional data file.

(CSV)

Click here for additional data file.

(CSV)

Click here for additional data file.

(TXT)

Click here for additional data file.

(CSV)

Click here for additional data file.

(R)

Click here for additional data file.

21 Feb 2020

PONE-D-20-01602

Kin recognition: neurogenomic response to mate choice and sib mating avoidance in a parasitic wasp

PLOS ONE

Dear Dr. Gallot,

Thank you for submitting your manuscript to PLOS ONE. After careful consideration, we feel that it has merit but does not fully meet PLOS ONE’s publication criteria as it currently stands. Therefore, we invite you to submit a revised version of the manuscript that addresses the points raised during the review process.

Two external referees and me have reviewed your manuscript and found some interesting results, but we agree that they paper needs major revisions and much clarification before possible publication. The reviewers have provided copious criticisms and comments. I will emphasize their major concerns of a rather small experimental design and sample sizes, and in some places very difficult to follow descriptions of the results. The main conclusions regarding gene modules and the actual transcriptional differences are not always well described and need to be further explained. One way to simplify the most important results would to be use a more conservative FDR P=0.01 instead of 0.05. Figure 4 is almost unreadable, and I agree the many of the figures need to be redrawn. The color schemes are not that helpful, and as one reviewer pointed out, are not likely to be viewer friendly for readers with color vision issues.

We would appreciate receiving your revised manuscript by Apr 06 2020 11:59PM. When you are ready to submit your revision, log on to https://www.editorialmanager.com/pone/ and select the 'Submissions Needing Revision' folder to locate your manuscript file.

If you would like to make changes to your financial disclosure, please include your updated statement in your cover letter.

To enhance the reproducibility of your results, we recommend that if applicable you deposit your laboratory protocols in protocols.io, where a protocol can be assigned its own identifier (DOI) such that it can be cited independently in the future. For instructions see: http://journals.plos.org/plosone/s/submission-guidelines#loc-laboratory-protocols

Please include the following items when submitting your revised manuscript:

A rebuttal letter that responds to each point raised by the academic editor and reviewer(s). This letter should be uploaded as separate file and labeled 'Response to Reviewers'.

A marked-up copy of your manuscript that highlights changes made to the original version. This file should be uploaded as separate file and labeled 'Revised Manuscript with Track Changes'.

An unmarked version of your revised paper without tracked changes. This file should be uploaded as separate file and labeled 'Manuscript'.

Please note while forming your response, if your article is accepted, you may have the opportunity to make the peer review history publicly available. The record will include editor decision letters (with reviews) and your responses to reviewer comments. If eligible, we will contact you to opt in or out.

We look forward to receiving your revised manuscript.

Kind regards,

William J. Etges

Academic Editor

PLOS ONE

Journal Requirements:

When submitting your revision, we need you to address these additional requirements.

1. Please ensure that your manuscript meets PLOS ONE's style requirements, including those for file naming. The PLOS ONE style templates can be found at

http://www.journals.plos.org/plosone/s/file?id=wjVg/PLOSOne_formatting_sample_main_body.pdf and http://www.journals.plos.org/plosone/s/file?id=ba62/PLOSOne_formatting_sample_title_authors_affiliations.pdf

[Note: HTML markup is below. Please do not edit.]

Reviewers' comments:

Reviewer's Responses to Questions

Comments to the Author

1. Is the manuscript technically sound, and do the data support the conclusions?

The manuscript must describe a technically sound piece of scientific research with data that supports the conclusions. Experiments must have been conducted rigorously, with appropriate controls, replication, and sample sizes. The conclusions must be drawn appropriately based on the data presented.

Reviewer #1: Partly

Reviewer #2: Partly

**********

2. Has the statistical analysis been performed appropriately and rigorously?

Reviewer #1: Yes

Reviewer #2: No

**********

3. Have the authors made all data underlying the findings in their manuscript fully available?

The PLOS Data policy requires authors to make all data underlying the findings described in their manuscript fully available without restriction, with rare exception (please refer to the Data Availability Statement in the manuscript PDF file). The data should be provided as part of the manuscript or its supporting information, or deposited to a public repository. For example, in addition to summary statistics, the data points behind means, medians and variance measures should be available. If there are restrictions on publicly sharing data—e.g. participant privacy or use of data from a third party—those must be specified.

Reviewer #1: No

Reviewer #2: Yes

**********

4. Is the manuscript presented in an intelligible fashion and written in standard English?

PLOS ONE does not copyedit accepted manuscripts, so the language in submitted articles must be clear, correct, and unambiguous. Any typographical or grammatical errors should be corrected at revision, so please note any specific errors here.

Reviewer #1: Yes

Reviewer #2: Yes

**********

5. Review Comments to the Author

Reviewer #1: The manuscript by Gallot and colleagues provides gene expression data from female heads that change in expression in response to courting and also that differ in their change in expression when courting is between related and unrelated parasitoid wasps. Females only mate once in this species, making the distinction between sib and non-sib very important, and the manuscript adds to phenotypic data showing females are able to distinguish kin and avoid mating with their brothers. Some good candidate genes whose expression changes accordingly are also reported.

One major issue with the paper is that the courtship assay is unclear. In particular, did all pairs produce courtship behaviour in each sample? It seems not from line 136, and line 103 mentions some females were not receptive. The experimental design thus seems to conflate receptivity with mating with sibs, since females have less receptivity to sibs. Was the proportion of successful courtship the same for sibs and non-sibs in each sample pool? Could the results be interpreted as a difference in harassment females experience between the treatments? A more detailed description of the assay is required in order to be able to comment on the discussion of the results.

In addition, the removal of one of the three biological replicates is problematic. The PCA analysis shows the samples have a more extreme profile towards the direction that their biological condition tends to. More discussion on the reasons the sample is considered unreliable is needed. Did less reads map to the reference? Was the RNA for this treatment collected at a different time of day than the other samples? More generally, the design of sample collection is unclear. Is it possible that circadian effects may have affected the design (for example if the courtship assays were performed at different times of the day or different days for the contrasted conditions?) Were the females from the same families in each pool, or is it possible that some conditions over-represent some families?

The GO enrichment tests were done against the full transcriptome (line 266). They should be done against the transcriptome with sufficient counts (i.e. only the 14,034 genes tested - 7 outliers that were excluded). This would be a more fair test.

I don't think it is a problem to avoid using a logFC threshold, but the authors should mention it explicitly in the text (line 240), and also provide some visualisation of how it would affect the results if one was employed. An easy way to do that is to provide a volcano plot (logCPM vs logFC) indicating the significant transcripts, so the reader can visualise the effect of any logFC threshold. The same figure would also be useful to illustrate the genes that are discussed in detail, such as the three neuropeptides in line 326.

Finally it would be good to provide the script and counts used for RNAseq analysis.

Overall I think this would be an interesting study to many readers. The main issue is to help interested readers compare the results to the rest of the literature. Both the provision of the script used for differential expression analysis, and a more detailed description of the conditions under which RNA was collected would help in comparing the results of this study to the literature.

Minor comments

Lines 71, 72, 76: sl-CSD -> sl-CSD locus (ie should indicate whether "sl-CSD" refers to the sex determination system or the genetic locus).

Line 79: please mention whether these experiments involved choice between two males, or not.

Please mention some statistics on the number of reads used in the RNAseq analysis, i.e. after filtering.

Line 321: should this number be 14,034 - 7?

Line 323: delete "levels", and clarify that 6.2% refers to tested gene set (I expected it to refer to the proportion 865/1,884)

Line 433: important -> significant

Line 436: I think the sentence "Thus ... state" is unnecessary as it repeats the message from previous sentences.

Line 447: delete "more precicely"

Line 451: within ten minutes -> "at most within 10 minutes". I think it is better to give a name to the condition captured, which is a better description of the biology e.g. "upon female acceptance of a mate" or "upon mounting by the male"

The red-green colour scheme in the figures is not colour-blind friendly.

I don't think panel A is needed for Fig 2 and Fig 3, the description in the text is sufficient. The legend is missing.

Fig 4: It is unclear what the "related control unrleated" text refers to in panel B. Also the text on the left part of panel B reads from top to bottom, while in other figures it reads from bottom to top. It would be best to rotate it by 180 degrees. More generally, please consider referring to the modules by a number rather than colours, which are too many to follow, especially for colour-blind readers.

This could be a major comment but I leave it to the authors to decide on whether it is worth following up. I think there is potential for an analysis of the interaction between genes that responded to courtship, and those that responded to courtship only with sibs or non-sibs.

For example one would make two virgin/courted contrasts: courted by sibs with virgins, and courted by non-sibs with virgins, and compare these results. There are many ways to illustrate this, for example

1. a scatterplot logFC virgin-sib courted vs logFC virgin -non sib courted or

2. 4-way venn: up courted with sib, up courted win non-sib, down courted with sib, down courted with non-sib.

This could indicate overall patterns of gene expression. For example are genes that are higher expressed when courted by sibs also expressed higher in courted by non-sibs, compared to virgins? Or do genes that are upregulated when mating with non-sibs resemble the expression pattern in virgins? This seems like the question in lines 102-104.

I appreciate Figure 4 is an alternative attempt to do this but if so it has not been discussed sufficiently to answer this question. One issue is that it is more difficult to appreciate the direction of change in gene expression in gene networks, but the information is there. At a minimum the authors could link the discussion of gene modules more directly with the original question in lines 102-104.

Reviewer #2: In general, this is an interesting study examining differential gene expression in females courted by related or unrelated males in a parasitoid wasp. This type of study is timely, and the authors do a good job presenting their research in the context of current work on both kin recognition and brain transcriptomics.

I am a little concerned about the small sample size (the authors have an N of 3 per treatment, plus technical replicates), and the authors apparent use of the technical replicates as independent samples. They need to be sure to pool the technical replicates for their statistical analyses, and to better articulate their actual sample size for their statistical analyses.

While their sample size is low, the results remain potentially interesting. I say potentially because I found the authors’ description of the different modules of genes associated with the different social conditions tested particularly difficult to read. The figures do not do a good job of illustrating what the clusters are, or how different gene expression is between the three treatments. The written description of these clusters was also difficult to read. This may be partially because the authors use the term “define” to describe these modules- but aren’t these modules of genes defined by the cluster analysis, not the authors themselves? Thus a better way to talk about these modules may be “X modules were associated with courtship….” That, coupled with a figure that better illustrates the different clusters, would make the results section seem less subjective. I suggest the authors revise their figures to show the different modules of genes that are expressed in these three different treatments. The authors should also describe which GO terms are associated with the different modules, not just state that there are GO terms associated with the modules.

A last major concern was that the authors were light on citations throughout their manuscript. I have mentioned a few of the many sentences missing citations in my minor comments below, but this is by no means exhaustive.

Minor/Specific points:

Line 37: Mates are also selected as a result of sensory bias (Fuller et al, 2005; Ryan & Cummings 2013). This hypothesis should be included in the introduction as well.

Line 59: missing a citation here.

Line 77: missing a citation here.

Line 79: missing a citation here.

Line 166: The number of biological replicates is quite low, particularly considering the number of treatments. That doesn’t mean the study isn’t useful, but does suggest a high likelihood of type II errors.

Line 204 and throughout: I suggest the authors clarify that they have 9 samples of 10 pooled heads each, plus 9 technical replicates. This is really a sample size of 9, with technical replicates, and should be analyzed accordingly. The technical replicates are not independent, thus the authors do not have a sample size of 18. Since these are technical replicates instead of independent individuals, they may skew the results, as they will amplify the differences between the three individuals in each treatment, and give an illusion of low within-treatment variance.

Lines 220-221: It is particularly disconcerting that the authors excluded one of their three samples from the sib-courting treatment. This brings the sample size for this treatment down to two, which is really too small for statistical analyses. Given that each of these samples contains the RNA from 10 heads, this sample dissimilarity suggests to me there was an error in the RNA extraction. However, an N of 2 is really too small for statistical analyses.

Line 283: Only 37.5% of the genes were most similar to Hymenopteran proteins? What were the rest of the genes most similar to? Drosophila or other insect proteins? A more detailed description of the gene set would be appreciated here.

Line 284: Only 33.4% of the genes were associated with at least one GO term? Is that because most of the genes were uncharacterized? Please clarify.

Line 288: What happens when the authors run these analyses with those outliers? Do most of the outlier reads map back to Hymenoptera? Or are they something else?

Line 291: Please define the groups by treatment and not color in the text (color should be in the figure legend, and remember, some readers will print this paper in black and white).

Line 384-388: The explanation of the visualization of these 50 different modules should be in the figure legend, not the main text.

Line 392: Naming the modules by colors is not particularly descriptive. Are there more biologically relevant names the authors could use for these modules?

Line 428: citation needed here.

Line 522: The genes constitute 90% of the royal jelly? Or are expression patterns of these genes responsible for 90% of the production and content of the royal jelly?

Line 538: Citation needed here.

Figure 4: The trees in panel A for figure 4 are so dense they are not readable. I suggest either removing them or changing the thickness of the lines.

**********

6. PLOS authors have the option to publish the peer review history of their article (what does this mean?). If published, this will include your full peer review and any attached files.

If you choose “no”, your identity will remain anonymous but your review may still be made public.

Do you want your identity to be public for this peer review? For information about this choice, including consent withdrawal, please see our Privacy Policy.

Reviewer #1: No

Reviewer #2: No

[NOTE: If reviewer comments were submitted as an attachment file, they will be attached to this email and accessible via the submission site. Please log into your account, locate the manuscript record, and check for the action link "View Attachments". If this link does not appear, there are no attachment files to be viewed.]

While revising your submission, please upload your figure files to the Preflight Analysis and Conversion Engine (PACE) digital diagnostic tool, https://pacev2.apexcovantage.com/. PACE helps ensure that figures meet PLOS requirements. To use PACE, you must first register as a user. Registration is free. Then, login and navigate to the UPLOAD tab, where you will find detailed instructions on how to use the tool. If you encounter any issues or have any questions when using PACE, please email us at figures@plos.org. Please note that Supporting Information files do not need this step.

7 Aug 2020

5. Review Comments to the Author

Reviewer #1: The manuscript by Gallot and colleagues provides gene expression data from female heads that change in expression in response to courting and also that differ in their change in expression when courting is between related and unrelated parasitoid wasps. Females only mate once in this species, making the distinction between sib and non-sib very important, and the manuscript adds to phenotypic data showing females are able to distinguish kin and avoid mating with their brothers. Some good candidate genes whose expression changes accordingly are also reported.

One major issue with the paper is that the courtship assay is unclear. In particular, did all pairs produce courtship behaviour in each sample? It seems not from line 136, and line 103 mentions some females were not receptive. The experimental design thus seems to conflate receptivity with mating with sibs, since females have less receptivity to sibs. Was the proportion of successful courtship the same for sibs and non-sibs in each sample pool? Could the results be interpreted as a difference in harassment females experience between the treatments? A more detailed description of the assay is required in order to be able to comment on the discussion of the results.

-The courtship behavior, i.e. a courtship parade produced by the male nearby the female, has been observed in all the pairs included in the study. Females did not experienced harassment by males, since males produce only a maximum of one courtship parade sequence. Whenever this behaviour sequence was not observed, the couple was eliminated from the study.

This information has been added to the manuscript (l. 147) “The 10 minutes period coincides with an active male courtship behaviour(42). For each female under conditions 1 and 2, one sequence of active courtship parade was observed within the 10 minutes following the introduction of the male to the box for all pairs. Whenever this behaviour sequence was not observed, the pair was eliminated from the study.”

Female receptivity has not been assessed the in this study. Therefore, the results on female receptivity (i.e. the increase in female receptivity in presence of unrelated males compared to in the presence of a related males) came from a study previously published by Metzger et al. (2010).

In addition, the removal of one of the three biological replicates is problematic. The PCA analysis shows the samples have a more extreme profile towards the direction that their biological condition tends to. More discussion on the reasons the sample is considered unreliable is needed. Did less reads map to the reference? Was the RNA for this treatment collected at a different time of day than the other samples? More generally, the design of sample collection is unclear. Is it possible that circadian effects may have affected the design (for example if the courtship assays were performed at different times of the day or different days for the contrasted conditions?) Were the females from the same families in each pool, or is it possible that some conditions over-represent some families?

-Concerning the removal of one sample (brother, replicate 1), we were unable to establish any difference explaining the extreme profile of this peculiar sample neither during sample collection neither in the data mapping process. The hypothesis that seems most plausible to us is that a problem occurred during the extraction of RNA for this peculiar sample. However, we were not able to identify differences in the reads mapped to the reference (one-way anova, p=0.895).

This sample, just as all the other is constituted of a pool of individual heads, and every individual in every samples have been randomized in the same way: “The experimental design has been conceived in order to minimize the influence of circadian rhythm and genetic background on the results. Newly emerged females were isolated every morning and numbered. A random draw was then made to establish the order of passage of the females and the treatment assigned. Behavioural experiments took place in the afternoon between 1:00 and 4:00 p.m. for all females that were captured in the morning in the order established by the random draw. A maximum of 2 females from each family have been kept daily, and randomized in one condition. For each given condition, the ten females belonged to different families to avoid an effect of genetic homogeneity on transcriptomic results. New families have been produced for every biological replicate.” We have now added these precisions in this second version of the manuscript (l. 138).

The GO enrichment tests were done against the full transcriptome (line 266). They should be done against the transcriptome with sufficient counts (i.e. only the 14,034 genes tested - 7 outliers that were excluded). This would be a more fair test.

-Thanks for this relevant suggestion. We have performed a new set of GO enrichment tests using the 14,027 transcripts (i.e. 14,034 minus 7 outliers) as a reference. All the results are included in supplementary tables, and the manuscript has been modified accordingly to these new results. These changes did not result in any major qualitative changes in our key results, which gives us confidence in the robustness of these results.

I don't think it is a problem to avoid using a logFC threshold, but the authors should mention it explicitly in the text (line 240), and also provide some visualisation of how it would affect the results if one was employed. An easy way to do that is to provide a volcano plot (logCPM vs logFC) indicating the significant transcripts, so the reader can visualise the effect of any logFC threshold. The same figure would also be useful to illustrate the genes that are discussed in detail, such as the three neuropeptides in line 326.

-We have added an explicit mention that we didn’t used a FC threshold (l. 259). We have also provided volcano plots for the 2 contrasts (figures 2A and 3A) where we have highlighted the genes discussed in detail.

Finally it would be good to provide the script and counts used for RNAseq analysis.

Overall I think this would be an interesting study to many readers. The main issue is to help interested readers compare the results to the rest of the literature. Both the provision of the script used for differential expression analysis, and a more detailed description of the conditions under which RNA was collected would help in comparing the results of this study to the literature.

-Raw data count table and R script have been added in supplementary material.

Minor comments

Lines 71, 72, 76: sl-CSD -> sl-CSD locus (ie should indicate whether "sl-CSD" refers to the sex determination system or the genetic locus).

-In the paragraph from l. 66 to l. 79 we have rephrased and indicated “sl-CSD” when it refers to the sex determination system and “sl-CSD locus” when it refers to the genetic locus.

Line 79: please mention whether these experiments involved choice between two males, or not.

-These experiments included the proportion of successful mate when a single female is in presence of sibs (i.e. presence of 2 brothers), or in presence of non sibs (i.e. presence of 2 unrelated males). Experiments involving choice between 2 males, one brother and one unrelated present in a same area have shown that female mates indifferently with the 2 males (Metzger et al., 2010). These points has now been explained lines 80-86.

Please mention some statistics on the number of reads used in the RNAseq analysis, i.e. after filtering.

-An overview of the RNAseq analysis statistics are summarized in the table S1. We have added some of those statistics in the results (l. 300): “After quality filtering a mean of 98.4% of paired-reads were kept, on which 70.8% were successfully mapped to the genome (representing a mean of 20 millions per sample, for a total of 363 million paired-reads, S1 table). Such values corresponded to the high-quality standards observed in other Hymenopteran species with an annotated genome(62).”

Line 321: should this number be 14,034 - 7?

-Indeed, we have corrected with the value 14,027 (l. 343)

Line 323: delete "levels", and clarify that 6.2% refers to tested gene set (I expected it to refer to the proportion 865/1,884)

-We have done the following modification l. 344-349 : “Among the 1,001 DEGs, 463 had higher expression in isolated females (3.3% of total transcriptome), gene ontology enrichment analysis reveals that this set of gene was enriched in DNA-binding Transcription Factor Activity (full list in S2 table). In contrast, 538 DEGs were overexpressed in courted females (3.8% of total transcriptome), gene ontology enrichment analysis reveals that this set of gene was enriched in Reproductive Behaviour (full list in S2 table).”

Line 433: important -> significant

-We have done the modification (l. 451)

Line 436: I think the sentence "Thus ... state" is unnecessary as it repeats the message from previous sentences.

-The sentence has been deleted (l. 454)

Line 447: delete "more precicely"

-We have done the modification (l. 467)

Line 451: within ten minutes -> "at most within 10 minutes". I think it is better to give a name to the condition captured, which is a better description of the biology e.g. "upon female acceptance of a mate" or "upon mounting by the male"

-The correction has been done (l. 470)

The red-green colour scheme in the figures is not colour-blind friendly.

-We have modified the color palette of all figures, using a colour-blind friendly palette (with blue, orange and grey colors).

I don't think panel A is needed for Fig 2 and Fig 3, the description in the text is sufficient. The legend is missing.

-Figures 2 and 3 have been modified by removing the A panels, that have been replaced by the volcano plots as suggested above. The legend is now included within the results section, in accordance with Plos one requirements.

Fig 4: It is unclear what the "related control unrleated" text refers to in panel B. Also the text on the left part of panel B reads from top to bottom, while in other figures it reads from bottom to top. It would be best to rotate it by 180 degrees. More generally, please consider referring to the modules by a number rather than colours, which are too many to follow, especially for colour-blind readers.

-Figure 4 has been modified “related” was replaced by “females courted by related males”; “unrelated” by “females courted by unrelated males”; and “control” by “control (isolated females)”. We also replaced this text previously located on the top of the heatmap, by the symbols “&”, “§” and “$” to increase clarity. Colours associated to modules have been conserved since the network illustration is based on colours attribution, but we considered the use of number rather than colour name, by numbering the all significant modules from 1 to 11. Then we refer to the gene module all along the manuscript only by using the numbered nomenclature. Thus, the manuscript as well as the figure 4 are now easier to follow, including for colour-blind readers.

This could be a major comment but I leave it to the authors to decide on whether it is worth following up. I think there is potential for an analysis of the interaction between genes that responded to courtship, and those that responded to courtship only with sibs or non-sibs.

For example one would make two virgin/courted contrasts: courted by sibs with virgins, and courted by non-sibs with virgins, and compare these results. There are many ways to illustrate this, for example

1. a scatterplot logFC virgin-sib courted vs logFC virgin -non sib courted or

2. 4-way venn: up courted with sib, up courted win non-sib, down courted with sib, down courted with non-sib.

This could indicate overall patterns of gene expression. For example are genes that are higher expressed when courted by sibs also expressed higher in courted by non-sibs, compared to virgins? Or do genes that are upregulated when mating with non-sibs resemble the expression pattern in virgins? This seems like the question in lines 102-104.

I appreciate Figure 4 is an alternative attempt to do this but if so it has not been discussed sufficiently to answer this question. One issue is that it is more difficult to appreciate the direction of change in gene expression in gene networks, but the information is there. At a minimum the authors could link the discussion of gene modules more directly with the original question in lines 102-104.

-Indeed, the purpose of the network analysis and the figure 4 (illustrating this analysis) is precisely to bring insight to this point. We are convinced that in the context of this study, such network approach is much more powerful than a set of bilateral comparisons. We considered this comment, together with the ones from the other reviewer and editor, and we understand that the results of this analysis need to be clarified in this manuscript. Thus we have proposed in this revised version of the manuscript a range of modifications accordingly:

o a new version of the figure 4 to illustrate the network analysis.

o a new description of the result paragraph

o modifications of the discussion directly oriented to the original question in the introduction (l. 102-104)

Concerning the last point, a part of the discussion has been added (l. 563-571):“We had formulated two non-mutually exclusive hypotheses. First, the perception of courtship was mediated by a change in gene expression, that would result in similar expression patterns in all females regardless of their relationship to the courting male. We identified such patterns for 2,780 genes (modules 7 and 8). Second, changes in female receptivity could result in changes in transcriptomic profiles. In this case, similar expression patterns would be expected for isolated females and females courted by their brothers. We have identified such expression patterns for 1,239 genes (modules 1, 2 and 3). Thus our results suggest that both courtship perception and changes in female receptivity induce a different neurogenomic response.”

Reviewer #2: In general, this is an interesting study examining differential gene expression in females courted by related or unrelated males in a parasitoid wasp. This type of study is timely, and the authors do a good job presenting their research in the context of current work on both kin recognition and brain transcriptomics.

I am a little concerned about the small sample size (the authors have an N of 3 per treatment, plus technical replicates), and the authors apparent use of the technical replicates as independent samples. They need to be sure to pool the technical replicates for their statistical analyses, and to better articulate their actual sample size for their statistical analyses.

-It seems we were not clear enough in the first version, since we did not used the technical replicates as independent samples. The technical replicates have been used independently only for the PCA and hierarchical clustering of libraries (figure 1). Then, each pair of technical replicates has been merged before proceeding to differential expression analysis and network analysis, as recommended. We have now added precisions on sample size in the material and methods (l. 217-222 and l. 241-244).

While their sample size is low, the results remain potentially interesting. I say potentially because I found the authors’ description of the different modules of genes associated with the different social conditions tested particularly difficult to read. The figures do not do a good job of illustrating what the clusters are, or how different gene expression is between the three treatments. The written description of these clusters was also difficult to read. This may be partially because the authors use the term “define” to describe these modules- but aren’t these modules of genes defined by the cluster analysis, not the authors themselves? Thus a better way to talk about these modules may be “X modules were associated with courtship….” That, coupled with a figure that better illustrates the different clusters, would make the results section seem less subjective. I suggest the authors revise their figures to show the different modules of genes that are expressed in these three different treatments. The authors should also describe which GO terms are associated with the different modules, not just state that there are GO terms associated with the modules.

-We understand that the results of this analysis need to be clarified in the manuscript.

o The description of the different modules of genes associated with the different social conditions has been modified following the reviewer suggestions (l. 400-419).

o The figure 4 illustrating the gene network analysis has been modified: i) the gene tree and gene network in the panel A have been changed in order to improve the clarity of how the clusters have been defined; ii) in the panel B significant modules have been numbered from 1 to 11; iii) the name of the social conditions have been modified.

o Some enriched GO terms associated with the significant modules have been highlighted in the table 1 within the manuscript, but exhaustive lists are still contained in supplementary table 4.

A last major concern was that the authors were light on citations throughout their manuscript. I have mentioned a few of the many sentences missing citations in my minor comments below, but this is by no means exhaustive.

-The revised version of the manuscript has been enriched with numerous citations throughout the manuscript.

Minor/Specific points:

Line 37: Mates are also selected as a result of sensory bias (Fuller et al, 2005; Ryan & Cummings 2013). This hypothesis should be included in the introduction as well.

-We have added this hypothesis in the introduction (l. 37): “The ‘good genes’ hypothesis predicts that females favour reproduction with males carrying traits that are honest indicators of good genes or as a result of sensory bias(3,4), hence obtaining genetic benefits for their offspring(5).”

Line 59: missing a citation here.

-A citation has been added here (l. 59-63) “Transcriptomic studies provided recent insight into female mating decisions(27). Coordinated changes in the expression of many genes in female brains, i.e., a neurogenomic response, have been identified following courtship in Poeciliidae fishes(28–30). This response depends on male attractiveness and is in accordance with female preferences.”

Line 77: missing a citation here.

-A citation has been added here (l. 77-79): “In the parasitoid wasp Venturia canescens, which has sl-CSD(38), females only mate once(39), making mate choice particularly decisive.”

Line 79: missing a citation here.

-A citation has been added here (l. 79-80): “Indeed, in this species, females are able to discriminate kin and non-kin during male courtship based on olfactory-mediated cues(40).”

Line 166: The number of biological replicates is quite low, particularly considering the number of treatments. That doesn’t mean the study isn’t useful, but does suggest a high likelihood of type II errors.

Line 204 and throughout: I suggest the authors clarify that they have 9 samples of 10 pooled heads each, plus 9 technical replicates. This is really a sample size of 9, with technical replicates, and should be analyzed accordingly. The technical replicates are not independent, thus the authors do not have a sample size of 18. Since these are technical replicates instead of independent individuals, they may skew the results, as they will amplify the differences between the three individuals in each treatment, and give an illusion of low within-treatment variance.

-We have clarified the number of samples (9 biological samples), each represented by 2 technical replicates (18 transcriptomic libraries) in material and methods. The 18 libraries were used only for PCA and hierarchical clustering. Then technical replicates have been merged for differential expression test and network analysis. We have now precisions on these points in the material and methods (l. 217-222 and l. 241-244).

Lines 220-221: It is particularly disconcerting that the authors excluded one of their three samples from the sib-courting treatment. This brings the sample size for this treatment down to two, which is really too small for statistical analyses. Given that each of these samples contains the RNA from 10 heads, this sample dissimilarity suggests to me there was an error in the RNA extraction. However, an N of 2 is really too small for statistical analyses.

-Concerning the removal of one sample (brother, replicate 1), we also speculate that the more probable explanation is that a problem occurred during RNA extraction. However, we were unable to establish any difference explaining the extreme profile of this peculiar sample neither during sample collection neither in the data mapping process (see our reply to a comment of referee 1). There were no differences in the proportion of reads mapped to the reference (one-way anova, p=0.895). This sample, just as all the other is constituted of a pool of 10 individual heads, and every individual in every samples have been randomized in the same way.

We are aware of the limitations related to the small sample size. We now clearly mentioned this limit in the discussion (l. 454-460) : “Despite the quite low number of biological replicates, the highly contrasted transcriptomes observed in the different social contexts suggest that sib mating avoidance behaviour could be considered a neurogenomic state. This research paves the way for further study on neurogenomic effects of sib mating avoidance in many species where such behaviours have been described and, thus, may contribute to the understanding of the molecular mechanisms underlying the evolution of avoiding consanguinity.”

Line 283: Only 37.5% of the genes were most similar to Hymenopteran proteins? What were the rest of the genes most similar to? Drosophila or other insect proteins? A more detailed description of the gene set would be appreciated here.

-The reviewer pointed a mistake in the data that has been now corrected in this version of the manuscript (l.337): “Overall, 76% of predicted genes get a blast hit (12,740) while 4,012 sequences get no hit. Among genes matching with blast, 89.4% presented their best hit with an insect sequence, of which 84.4% match more specifically to a hymenopteran insect sequences.”

Line 284: Only 33.4% of the genes were associated with at least one GO term? Is that because most of the genes were uncharacterized? Please clarify.

-Indeed only 33.4% of the predicted genes in V. canescens genome get at least one GO annotation. This means that while a majority (76%) of predicted genes shown homology with protein sequences that have been previously deposited in databases, i.e., get a blast hit; only a low proportion of these sequences are associated to functional annotation (molecular function, cellular component, biological process). Such gene biological identity are scarce and essentially come from experiments performed on model organisms. Most of the genes were uncharacterized, notably in non-model metazoan species. Nonetheless, this proportion is rapidly increasing. In comparison, in human, only 32% of genes got at least one GO annotation in 2004; but in 2015, 65% of genes got at least one GO annotation (Tomcsak et al., 2018).

We have added precision in the manuscript (l. 308): “Finally, most of the genes were uncharacterized, since only 33.4% of the predicted genes (5,589) were associated with at least one Gene Ontology (GO) functional annotation.”

Line 288: What happens when the authors run these analyses with those outliers? Do most of the outlier reads map back to Hymenoptera? Or are they something else?

-We run again these analyses without those 7 outliers, since the 7 corresponding sequences do not match with any known sequences. All the results relative to differential expression and network analysis in this new version of the manuscript have been done on the dataset containing 14,027 genes.

Line 291: Please define the groups by treatment and not color in the text (color should be in the figure legend, and remember, some readers will print this paper in black and white).

-The correction has been done (l. 315)

Line 384-388: The explanation of the visualization of these 50 different modules should be in the figure legend, not the main text.

-The correction has been done (l. 408)

Line 392: Naming the modules by colors is not particularly descriptive. Are there more biologically relevant names the authors could use for these modules?

-We have now proposed to use number to characterize modules rather than color. Particularly, we have numbered from 1 to 11 the modules significantly correlated with at least one experimental condition, while other modules which are not mentioned in the manuscript have not been numbered. However, color names have been conserved in parenthesis, because the colors are used in the figure 4 to illustrate the gene network.

Line 428: citation needed here.

-A citation has been added here (l. 447) “I In V. canescens females, mate relatedness influences female sexual receptivity and is estimated during male courtship displays through chemical cues(40).”

Line 522: The genes constitute 90% of the royal jelly? Or are expression patterns of these genes responsible for 90% of the production and content of the royal jelly?

-This sentence was incorrect and has been rectified (l. 544): “royal jelly is constituted with 90% of MRJP proteins”.

Line 538: Citation needed here.

-A citation has been added here (l. 557): “Sustained kin odourants exposure during development drives changes in neurotransmitter expression from GABA to dopamine neurons, which are stimulated from an increase in the expression of the transcription factor PAX6 and accompanied by a behavioural preference for kin odourants(83).”

Figure 4: The trees in panel A for figure 4 are so dense they are not readable. I suggest either removing them or changing the thickness of the lines.

-We proposed a new illustration of the gene network in the figure 4. The tree panel resolution has been increased.

7 Sep 2020

PONE-D-20-01602R1

Kin recognition: neurogenomic response to mate choice and sib mating avoidance in a parasitic wasp

PLOS ONE

Dear Dr. gallot,

Thank you for submitting your manuscript to PLOS ONE.  We invite you to submit a revised version of the manuscript that addresses the points raised during the review process.

I agree with both external referees that the authors have done a solid job of responding to the reviewers' comments and criticisms. Both reviewers have suggested a few more issues that need addressing before a final decision can be made.  I think the MS is much clearer and the results are better described for your readers.

Please submit your revised manuscript by Oct 22 2020 11:59PM. If you will need more time than this to complete your revisions, please reply to this message or contact the journal office at plosone@plos.org. When you're ready to submit your revision, log on to https://www.editorialmanager.com/pone/ and select the 'Submissions Needing Revision' folder to locate your manuscript file.

Please include the following items when submitting your revised manuscript:

A rebuttal letter that responds to each point raised by the academic editor and reviewer(s). You should upload this letter as a separate file labeled 'Response to Reviewers'.

A marked-up copy of your manuscript that highlights changes made to the original version. You should upload this as a separate file labeled 'Revised Manuscript with Track Changes'.

An unmarked version of your revised paper without tracked changes. You should upload this as a separate file labeled 'Manuscript'.

If you would like to make changes to your financial disclosure, please include your updated statement in your cover letter. Guidelines for resubmitting your figure files are available below the reviewer comments at the end of this letter.

If applicable, we recommend that you deposit your laboratory protocols in protocols.io to enhance the reproducibility of your results. Protocols.io assigns your protocol its own identifier (DOI) so that it can be cited independently in the future. For instructions see: http://journals.plos.org/plosone/s/submission-guidelines#loc-laboratory-protocols

We look forward to receiving your revised manuscript.

Kind regards,

William J. Etges

Academic Editor

PLOS ONE

[Note: HTML markup is below. Please do not edit.]

Reviewers' comments:

Reviewer's Responses to Questions

Comments to the Author

1. If the authors have adequately addressed your comments raised in a previous round of review and you feel that this manuscript is now acceptable for publication, you may indicate that here to bypass the “Comments to the Author” section, enter your conflict of interest statement in the “Confidential to Editor” section, and submit your "Accept" recommendation.

Reviewer #1: All comments have been addressed

Reviewer #2: (No Response)

**********

2. Is the manuscript technically sound, and do the data support the conclusions?

The manuscript must describe a technically sound piece of scientific research with data that supports the conclusions. Experiments must have been conducted rigorously, with appropriate controls, replication, and sample sizes. The conclusions must be drawn appropriately based on the data presented.

Reviewer #1: Yes

Reviewer #2: Yes

**********

3. Has the statistical analysis been performed appropriately and rigorously?

Reviewer #1: Yes

Reviewer #2: Yes

**********

4. Have the authors made all data underlying the findings in their manuscript fully available?

The PLOS Data policy requires authors to make all data underlying the findings described in their manuscript fully available without restriction, with rare exception (please refer to the Data Availability Statement in the manuscript PDF file). The data should be provided as part of the manuscript or its supporting information, or deposited to a public repository. For example, in addition to summary statistics, the data points behind means, medians and variance measures should be available. If there are restrictions on publicly sharing data—e.g. participant privacy or use of data from a third party—those must be specified.

Reviewer #1: Yes

Reviewer #2: Yes

**********

5. Is the manuscript presented in an intelligible fashion and written in standard English?

PLOS ONE does not copyedit accepted manuscripts, so the language in submitted articles must be clear, correct, and unambiguous. Any typographical or grammatical errors should be corrected at revision, so please note any specific errors here.

Reviewer #1: Yes

Reviewer #2: Yes

**********

6. Review Comments to the Author

Please use the space provided to explain your answers to the questions above. You may also include additional comments for the author, including concerns about dual publication, research ethics, or publication ethics. (Please upload your review as an attachment if it exceeds 20,000 characters)

Reviewer #1: The authors have responded to all my comments in a satisfactory manner.

There is one exception with regards to the sample that was removed. The text should mention why the sample represents a problem during preparation, and not due to biological variation, given that both technical replicates are very similar, and both are outliers.

Minor comments

Line 259:did not applied -> did not apply

Line 289: this should not be the full transcriptome, but the part of the transcriptome that was included in the analysis. The authors mention that they corrected for this, so I suspect this is text that has not been updated.

Lines 356-358: Please mention in the text the direction of high expression for the genes that are discussed.

Line 425: Showed -> shows, illustrated -> illustrates, denoted -> denotes, indicated -> indicates

Line 480: Veltos -> Veltsos

Fig 3A: Please check the top left outlier gene. Is it reliable? If not consider removing it to shorten the y axis.

Fig 4: the symbols and numbers overap, please use smaller text.

Reviewer #2: This is a much improved version of a previously reviewed manuscript: “Kin recognition: neurogenomic response to mate choice and sib mating avoidance in a parasitic wasp”. I have a few remaining minor comments that the authors should address.

Minor points:

Line 241: “All pairs of technical replicates were merged before….” Not “has been merged”.

Line 249: missing the words “of all” between expression and transcripts. It should read: “tested for differential expression of all transcripts with an average level of expression above 10 reads per library”

Figure 1 legend: Please state which color represents which treatment in the figure legend instead of stating that treatments are represented by different colors.

Line 381: This section is a little unclear. Are all 22 of these GO terms only associated with 2 genes? Or are they associated with some of the other 479 differentially expressed genes as well?

Lines 549-550: The authors need to mention the decades of research on the molecular pathways involved in kin-recognition in Hydractinia symbiolongicarpus and Botryllus schosseri. This topic has been well studied in colonial marine invertebrates since at least 1994 (Mokady & Buss, 1994; Fagan & Weissman, 1997; Cadavid et al., 2004; Litman, 2006; Nicotra et al., 2009; Nydam et al., 2017, to name a few of the many papers on this topic).

It looks like the authors included both the original 4 figures and the revised figures, but failed to mention which figures were the old figures and which figures were the new figures. It looks like the second set of figures are the new figures. Assuming that is true, the color scheme is definitely more color-blind friendly, though Figure 4b is still a little difficult to read. I suggest stretching the y-axis so the cluster numbers are easier to match to bars.

The new S1 figure is useful, as is S3. I appreciate their addition.

**********

7. PLOS authors have the option to publish the peer review history of their article (what does this mean?). If published, this will include your full peer review and any attached files.

If you choose “no”, your identity will remain anonymous but your review may still be made public.

Do you want your identity to be public for this peer review? For information about this choice, including consent withdrawal, please see our Privacy Policy.

Reviewer #1: No

Reviewer #2: No

[NOTE: If reviewer comments were submitted as an attachment file, they will be attached to this email and accessible via the submission site. Please log into your account, locate the manuscript record, and check for the action link "View Attachments". If this link does not appear, there are no attachment files.]

While revising your submission, please upload your figure files to the Preflight Analysis and Conversion Engine (PACE) digital diagnostic tool, https://pacev2.apexcovantage.com/. PACE helps ensure that figures meet PLOS requirements. To use PACE, you must first register as a user. Registration is free. Then, login and navigate to the UPLOAD tab, where you will find detailed instructions on how to use the tool. If you encounter any issues or have any questions when using PACE, please email PLOS at figures@plos.org. Please note that Supporting Information files do not need this step.

5 Oct 2020

Dear William J. Etges,

We want to thank the reviewers as well as the editor for their corrections. We hope that this revised version includes all the changes requested.

A comprehensive response to the specific items provided by the 2 reviewers has been addressed. We have responded to each item in details in the following section (part 6).

On behalf of the authors,

Aurore Gallot, PhD.

PONE-D-20-01602R1

Kin recognition: neurogenomic response to mate choice and sib mating avoidance in a parasitic wasp

PLOS ONE

Dear Dr. gallot,

Thank you for submitting your manuscript to PLOS ONE. We invite you to submit a revised version of the manuscript that addresses the points raised during the review process.

I agree with both external referees that the authors have done a solid job of responding to the reviewers' comments and criticisms. Both reviewers have suggested a few more issues that need addressing before a final decision can be made. I think the MS is much clearer and the results are better described for your readers.

Please submit your revised manuscript by Oct 22 2020 11:59PM. If you will need more time than this to complete your revisions, please reply to this message or contact the journal office at plosone@plos.org. When you're ready to submit your revision, log on to https://www.editorialmanager.com/pone/ and select the 'Submissions Needing Revision' folder to locate your manuscript file.

Please include the following items when submitting your revised manuscript:

* A rebuttal letter that responds to each point raised by the academic editor and reviewer(s). You should upload this letter as a separate file labeled 'Response to Reviewers'.

* A marked-up copy of your manuscript that highlights changes made to the original version. You should upload this as a separate file labeled 'Revised Manuscript with Track Changes'.

* An unmarked version of your revised paper without tracked changes. You should upload this as a separate file labeled 'Manuscript'.

If you would like to make changes to your financial disclosure, please include your updated statement in your cover letter. Guidelines for resubmitting your figure files are available below the reviewer comments at the end of this letter.

If applicable, we recommend that you deposit your laboratory protocols in protocols.io to enhance the reproducibility of your results. Protocols.io assigns your protocol its own identifier (DOI) so that it can be cited independently in the future. For instructions see: http://journals.plos.org/plosone/s/submission-guidelines#loc-laboratory-protocols

We look forward to receiving your revised manuscript.

Kind regards,

William J. Etges

Academic Editor

PLOS ONE

Reviewers' comments:

Reviewer's Responses to Questions

Comments to the Author

1. If the authors have adequately addressed your comments raised in a previous round of review and you feel that this manuscript is now acceptable for publication, you may indicate that here to bypass the “Comments to the Author” section, enter your conflict of interest statement in the “Confidential to Editor” section, and submit your "Accept" recommendation.

Reviewer #1: All comments have been addressed

Reviewer #2: (No Response)

________________________________

2. Is the manuscript technically sound, and do the data support the conclusions?

The manuscript must describe a technically sound piece of scientific research with data that supports the conclusions. Experiments must have been conducted rigorously, with appropriate controls, replication, and sample sizes. The conclusions must be drawn appropriately based on the data presented.

Reviewer #1: Yes

Reviewer #2: Yes

________________________________

3. Has the statistical analysis been performed appropriately and rigorously?

Reviewer #1: Yes

Reviewer #2: Yes

________________________________

4. Have the authors made all data underlying the findings in their manuscript fully available?

The PLOS Data policy requires authors to make all data underlying the findings described in their manuscript fully available without restriction, with rare exception (please refer to the Data Availability Statement in the manuscript PDF file). The data should be provided as part of the manuscript or its supporting information, or deposited to a public repository. For example, in addition to summary statistics, the data points behind means, medians and variance measures should be available. If there are restrictions on publicly sharing data—e.g. participant privacy or use of data from a third party—those must be specified.

Reviewer #1: Yes

Reviewer #2: Yes

________________________________

5. Is the manuscript presented in an intelligible fashion and written in standard English?

PLOS ONE does not copyedit accepted manuscripts, so the language in submitted articles must be clear, correct, and unambiguous. Any typographical or grammatical errors should be corrected at revision, so please note any specific errors here.

Reviewer #1: Yes

Reviewer #2: Yes

________________________________

6. Review Comments to the Author

Please use the space provided to explain your answers to the questions above. You may also include additional comments for the author, including concerns about dual publication, research ethics, or publication ethics. (Please upload your review as an attachment if it exceeds 20,000 characters)

Reviewer #1: The authors have responded to all my comments in a satisfactory manner.

There is one exception with regards to the sample that was removed. The text should mention why the sample represents a problem during preparation, and not due to biological variation, given that both technical replicates are very similar, and both are outliers.

We explicitly added in the manuscript why the outlier sample is problematic and has been removed (l. 248-249).

Minor comments

All the minor comments has been taken in account in the new version of the manuscript

Line 259:did not applied -> did not apply

We have corrected (l. 271).

Line 289: this should not be the full transcriptome, but the part of the transcriptome that was included in the analysis. The authors mention that they corrected for this, so I suspect this is text that has not been updated.

We have corrected (l. 302).

Lines 356-358: Please mention in the text the direction of high expression for the genes that are discussed.

The direction of high expression is now described for all the discussed genes. We added this information for the genes ebony, DAT, 5-HT1A and serotonin receptor 1A (l. 372-375).

Line 425: Showed -> shows, illustrated -> illustrates, denoted -> denotes, indicated -> indicates

We have corrected (l. 444-447).

Line 480: Veltos -> Veltsos

We have corrected (l. 505).

Fig 3A: Please check the top left outlier gene. Is it reliable? If not consider removing it to shorten the y axis.

We controlled the outlier gene but we do not considered that it is aberrant given that the corresponding sequence encode and open reading frame that perfectly match to the mitochondrial gene CO1. Considering that the gene annotation is reliable, we did not remove the outlier point, and did not rescale the figure either.

Fig 4: the symbols and numbers overap, please use smaller text.

We stretched the y axis and used a smaller text in the figure 4b.

Reviewer #2: This is a much improved version of a previously reviewed manuscript: “Kin recognition: neurogenomic response to mate choice and sib mating avoidance in a parasitic wasp”. I have a few remaining minor comments that the authors should address.

Minor points:

Line 241: “All pairs of technical replicates were merged before….” Not “has been merged”.

We have corrected (l. 252).

Line 249: missing the words “of all” between expression and transcripts. It should read: “tested for differential expression of all transcripts with an average level of expression above 10 reads per library”

We have corrected (l. 261).

Figure 1 legend: Please state which color represents which treatment in the figure legend instead of stating that treatments are represented by different colors.

We added which color represent which treatment in the figure legend (l. 350-351).

Line 381: This section is a little unclear. Are all 22 of these GO terms only associated with 2 genes? Or are they associated with some of the other 479 differentially expressed genes as well?

Only the 2 GO ‘Reproductive Behaviour’ and ‘Male Mating Behaviour’ are associated with 2 genes. We clarified this section (l. 401).

Lines 549-550: The authors need to mention the decades of research on the molecular pathways involved in kin-recognition in Hydractinia symbiolongicarpus and Botryllus schosseri. This topic has been well studied in colonial marine invertebrates since at least 1994 (Mokady & Buss, 1994; Fagan & Weissman, 1997; Cadavid et al., 2004; Litman, 2006; Nicotra et al., 2009; Nydam et al., 2017, to name a few of the many papers on this topic).

We mentioned research on kin recognition in colonial invertebrates in the introduction (l. 47-51).

It looks like the authors included both the original 4 figures and the revised figures, but failed to mention which figures were the old figures and which figures were the new figures. It looks like the second set of figures are the new figures. Assuming that is true, the color scheme is definitely more color-blind friendly, though Figure 4b is still a little difficult to read. I suggest stretching the y-axis so the cluster numbers are easier to match to bars.

We stretched the y axis and used a smaller text in the figure 4b.

The new S1 figure is useful, as is S3. I appreciate their addition.

Submitted filename: Response_to_reviewers.docx

Click here for additional data file.

9 Oct 2020

Kin recognition: neurogenomic response to mate choice and sib mating avoidance in a parasitic wasp

PONE-D-20-01602R2

Dear Dr. gallot,

We’re pleased to inform you that your manuscript has been judged scientifically suitable for publication and will be formally accepted for publication once it meets all outstanding technical requirements.

Within one week, you’ll receive an e-mail detailing the required amendments. When these have been addressed, you’ll receive a formal acceptance letter and your manuscript will be scheduled for publication.

An invoice for payment will follow shortly after the formal acceptance. To ensure an efficient process, please log into Editorial Manager at http://www.editorialmanager.com/pone/, click the 'Update My Information' link at the top of the page, and double check that your user information is up-to-date. If you have any billing related questions, please contact our Author Billing department directly at authorbilling@plos.org.

If your institution or institutions have a press office, please notify them about your upcoming paper to help maximize its impact. If they’ll be preparing press materials, please inform our press team as soon as possible -- no later than 48 hours after receiving the formal acceptance. Your manuscript will remain under strict press embargo until 2 pm Eastern Time on the date of publication. For more information, please contact onepress@plos.org.

Kind regards,

William J. Etges

Academic Editor

PLOS ONE

Additional Editor Comments (optional):

Reviewers' comments:

16 Oct 2020

PONE-D-20-01602R2

Kin recognition: neurogenomic response to mate choice and sib mating avoidance in a parasitic wasp

Dear Dr. gallot:

I'm pleased to inform you that your manuscript has been deemed suitable for publication in PLOS ONE. Congratulations! Your manuscript is now with our production department.

If your institution or institutions have a press office, please let them know about your upcoming paper now to help maximize its impact. If they'll be preparing press materials, please inform our press team within the next 48 hours. Your manuscript will remain under strict press embargo until 2 pm Eastern Time on the date of publication. For more information please contact onepress@plos.org.

If we can help with anything else, please email us at plosone@plos.org.

Thank you for submitting your work to PLOS ONE and supporting open access.

Kind regards,

PLOS ONE Editorial Office Staff

on behalf of

Dr. William J. Etges

Academic Editor

PLOS ONE