Introduction
============

Understanding the evolutionary mechanisms underlying sexual dimorphism has been a very challenging task.  In this regard, an important question is how two individuals of different sexes in a species may have conspicuous sexual variation, even when both sexes share practically the same genome.  Transcriptome studies indicate that most morphological sex differences are caused by divergent patterns of gene expression between sexes ([@evy009-B24]).  These are referred to as sex-biased genes, which have consistently shown rapid sequence evolutionary rates across taxa ([@evy009-B61]; [@evy009-B67]; [@evy009-B45]; [@evy009-B111]; [@evy009-B73]).  In *Drosophila*, male-biased expressed genes evolve particularly fast ([@evy009-B24]), which is mainly caused by adaptive evolution ([@evy009-B80]). 

Potential explanations for such phenomenon involve sperm competition, sexual selection, and/or sexual conflict ([@evy009-B101]).  If this hypothesis is at least partially true, some of the products of these genes might elicit pre- or postmating barriers which may, ultimately, play important roles reinforcing species boundaries ([@evy009-B97]; [@evy009-B31]).  Indeed, accessory gland proteins secreted in males' seminal fluid in *Drosophila* influence the females' physiology and behavior ([@evy009-B22]; [@evy009-B82]) and tend to evolve under positive selection ([@evy009-B100]), which may reflect a role on reproductive isolation in the first stages of speciation.  Furthermore, several female reproductive proteins from *Drosophila* have also been reported to evolve under positive selection ([@evy009-B102]; [@evy009-B72]).  Female proteins that interact with rapidly evolving male proteins may evolve faster because of coevolution ([@evy009-B36]).  In addition, proteins in the external layer of the eggshell (chorion) have been reported to evolve adaptively due to possible role in the sperm--egg and/or egg--environment interactions ([@evy009-B48]). 

Despite the evidence of the contribution of positive selection on sex-biased genes in *Drosophila*, there are alternative evolutionary explanations for rapid evolution on such genes.  Studies have demonstrated that selection is weakened when trait (or gene) expression is limited to a fraction of individuals such as sex-biased genes, resulting in an increased segregation of slightly deleterious variation, which can reach fixation by genetic drift ([@evy009-B107]; [@evy009-B81]).  As a consequence, not only polymorphism levels are increased on such genes, but they also evolve faster.  In fact, it has been demonstrated that relaxed constraints, genetic drift, or/and an increased segregation of slightly deleterious variation have an important impact on the evolution of male-specific genes ([@evy009-B32]; [@evy009-B37]).  Furthermore, sex-biased genes tend to have narrower expression pattern than unbiased genes, that may imply less pleiotropy and functional constraints, showing faster evolution because of relaxed purifying selection ([@evy009-B62]). 

Such as in reproductive tissues, other tissues may also express reproductive-related proteins.  Sex pheromones are important in controlling reproductive behavior ([@evy009-B42]).  Pheromones and other environmental olfactory cues are perceived as taste and olfactory stimuli and then processed by the chemosensory system in organs located mainly in the head (such as antennae) ([@evy009-B54]).  Among the genes involved in this process, there are sets of gene families that encode for proteins involved in ligand-binding (odorant binding proteins and chemosensory proteins) and receptor functions (odorant receptors, gustatory receptors, ionotropic receptors, and sensory neuron membrane proteins) ([@evy009-B49]; [@evy009-B90]).  The molecular evolution of these protein families has been widely studied in insects and has revealed that several of these genes evolve under positive selection under a birth-and-death process that leads to a rapid gene turnover ([@evy009-B91]; [@evy009-B7]; [@evy009-B10]). 

Here we investigate genes expressed in reproductive and head tissues of two closely related species, South American fruit flies (*Anastrepha fraterculus*) and West Indies fruit flies (*Anastrepha obliqua*), which belong to the *fraterculus* group ([@evy009-B70]).  Taxonomic identification of some species within this group based only on morphology is difficult due to overlapping variation even in the aculeus, which is one of the key traits in the systematics of this group ([@evy009-B123]; [@evy009-B77]).  Molecular phylogeny of the *fraterculus* group based on the mitochondrial gene COI showed polyphyly for these two species ([@evy009-B96]).  However, phylogenetic analyses using nuclear *loci* revealed that *A. obliqua* is a monophyletic lineage ([@evy009-B92]) and not as closely related to *A. fraterculus* as other species in the group, though there is evidence of historical introgression between these lineages ([@evy009-B92]; Díaz F, Lima ALA, Nakamura AM, Fernandes F, Sobrinho I, de Brito RA, unpublished data).  Furthermore, these species may produce viable hybrids in laboratory with descendants of some combinations obeying Haldane's rule ([@evy009-B21]), and as they are found in sympatry in several regions, it is possible that current introgression may still occur in nature.  Therefore, it is possible that *A. fraterculus* and *A. obliqua* have diverged recently while retaining some gene flow, emphasizing the importance of identifying genomic regions that responded to selection and may have had a leading role on their differentiation as has been proposed for other organisms with similar speciation patterns ([@evy009-B28]). 

We generated transcriptomes of reproductive tissues from *A. fraterculus* and *A. obliqua* and compared with RNA-seq data produced from head tissues of the same populations and species ([@evy009-B84]).  In this study, we estimated differentiation indexes and evolutionary rates from pairwise comparisons between both species and among seven species of Tephritidae.  In addition, we tested for signals of natural selection and relaxed constraints.  These results enabled us to identify which tissues, reproductive or cephalic, and sex, would show genes with higher evolutionary rates and whether this is due to positive selection or nonadaptive evolution.  Answers to these questions not only contribute to the understanding of the evolutionary mechanisms affecting sex-biased genes, but also may offer clues to the differentiation process influencing these fruit flies even in the presence of gene flow. 

Materials and Methods
=====================

Sampling and Laboratory Procedures
----------------------------------

Individuals of *A. fraterculus* were collected from the field from guava (Myrtaceae) fruits (22°01′03″S, 47°53′27″W) and *A. obliqua* from jocote (Anacardiaceae) fruits (16°41′58″S, 49°16′35″W).  These populations were maintained in laboratory under the following controlled conditions: 26 ± 1°C of temperature, 60--90% of humidity, and natural photoperiod.  Reproductive tissues of virgin adult (8--12 days) male (testis, accessory glands, and phallus) and female (ovaries, accessory glands, spermatheca, uterus, and ovipositor) flies of *A. fraterculus* and *A. obliqua* were collected.  Total RNA was extracted from a pool of reproductive tissues of five individuals following the protocol proposed by [@evy009-B14].  After extraction, each sample was formed by an equimolar mix of two pools, totaling samples from ten individuals in every mix.  Biological replicates were prepared for each profile (species, sex, and tissue), making for a total of eight samples.  RNA-seq libraries were constructed using the TruSeq RNA Sample Preparation Kit (Illumina) protocol according to the manufacturer\'s instructions.  Libraries of 2 × 100 bp paired-end reads were sequenced on Illumina HiSeq2000 and HiScan at the Laboratory of Functional Genomics Applied to Agriculture and Agri-energy, ESALQ-USP, Brazil. 

Cleaning and Assembly
---------------------

Reads obtained from sequencing of reproductive tissues as well as published transcriptomes from head tissues of the pair of species studied here ([@evy009-B84]) were trimmed using the program Trimmomatic v.0.33 ([@evy009-B6]), setting the parameters LEADING: 5 TRAILING: 5 SLIDINGWINDOW: 5: 20 MINLEN: 50.  This program also searches for and removes any remaining TrueSeq Illumina adapters in the reads.  Unpaired reads were also discarded.  After this censoring, reads from the same species were joined to produce two assemblies.  Each group of reads was normalized by coverage and assembled using default parameters of Trinity v.2.4.0 ([@evy009-B35]). 

Unigene Prediction and Assessment of the Quality of Assemblies
--------------------------------------------------------------

We searched for potential coding sequences (CDSs) in all six open reading frames of each transcript using the software TransDecoder v.3.0.1 (<http://transdecoder.github.io>) following three steps.  First, TransDecoder. LongOrfs was used to retain all potential CDSs coding peptides longer than 100 amino acids.  In the second step, these peptides were submitted to the hmmscan tool included in the HMMER v.3.1b2 package ([@evy009-B23]) to search for protein signatures in the Pfam-A database and BLASTP v.2.6 ([@evy009-B9]) to search for similar sequences in the nonredundant database of the GenBank (nonredundant \[nr\]) including only proteins of arthropods.  In the third step, the program TransDecoder. Predict uses the information produced by the other steps to predict the CDSs.  Redundancy of the obtained CDSs was reduced using Cd-hit-est ([@evy009-B29]) with a similarity threshold of 0.99.  To obtain the final set of putative unigenes, transcripts with these CDSs were filtered using the Trinity assembly information and only the isoform with the highest expression per trinity component was retained.  For that, the reads from each species were mapped to the respective assembly using Bowtie2 ([@evy009-B57]) and the abundance of each transcript was estimated by eXpress v.1.5.1 ([@evy009-B85]).  These steps were performed by the script align_and_estimate_abundance.pl included in the Trinity package, adding no bias correction option for the eXpress program.  The completeness and redundancy level of the raw and filtered assemblies of each species was evaluated by BUSCO (Benchmarking Universal Single-Copy Orthologs) ([@evy009-B94]) using the Arthropoda database as reference. 

Functional Annotation
---------------------

Predicted unigene CDSs were compared against the GenBank nr protein database including only arthropod proteins, the *Drosophila melanogaster* protein database (r6.14), and the Eukaryotic Orthologous Groups of proteins database (KOG) ([@evy009-B55]) using NCBI BLASTP v.2.6 ([@evy009-B9]).  To all these analyses, we set an *e*-value threshold of 10^−6^. We also searched for conserved protein domains using InterProScan 5.24-63.0 ([@evy009-B50]).  Annotations against nr and conserved protein domains databases were submitted to Blast2GO ([@evy009-B15]) to obtain a list of gene ontology (GO) terms associated with the annotated genes.  Frequencies of GO terms at the level 2 were obtained using the program WEGO ([@evy009-B116]) and their distributions were plotted using GO terms with frequencies greater than 1%. 

Identifying Sex- and Tissue-Biased Unigenes
-------------------------------------------

The sex-biased genes were identified by comparing the expressions of male and female profiles in each tissue (reproductive and head) and species (*A. fraterculus* and *A. obliqua*).  In addition, tissue-biased genes were identified by comparing the expressions of head and reproductive tissues in each sex and species.  Expression analysis was performed by using the scripts align_and_estimate_abundance.pl, abundance_estimates_to_matrix.pl, PtR, and analyze_diff_expr.pl provided by the Trinity package ([@evy009-B35]).  In the align_and_estimate_abundance.pl, we used Bowtie2 ([@evy009-B57]) and eXpress v1.5.1 ([@evy009-B85]) to map the reads back to each species\' assemblies (set of unigenes) and to estimate abundances of each unigene, respectively.  This script was run adding no bias correction option for the eXpress program and very-sensitive option to Bowtie2.  The abundance_estimates_to_matrix.pl script put the abundances values estimated to each RNA-seq library in a matrix.  The PtR program was used to verify the quality of the biological replicates using Pearson correlation and principal component analysis of the unigenes expression across samples measured as log~2~ transformed of counts per million.  Differential gene expression analysis among sexes and tissues was performed in edgeR ([@evy009-B86]) using the TMM (trimmed mean of M-values) normalized abundances.  Expression values are shown in transcripts per million (TPM).  Unigenes with fold-changes greater than 4 and a significance of FDR corrected *P*-values smaller than 0.05 were considered as differentially expressed. 

Single Nucleotide Polymorphism Calling, Differentiation Indexes, and the McDonald--Kreitman Test
------------------------------------------------------------------------------------------------

*Anastrepha fraterculus* unigenes were used as reference for single nucleotide polymorphism (SNP) calling.  Filtered reads from each library were mapped to each assembly according to tissue using Bowtie2 ([@evy009-B57]).  Mapped reads were converted to mpileup format and filtered based on minimum mapping quality of 20 and minimum PHRED quality of 30 using mpileup tool provided by Samtools v.1.3.1 package ([@evy009-B58]), minimum coverage of 20, minimum reads of 1 to call the variant and strand filter (removed variants with more than 90% supported by only reads of one strand) using the tool mpileup2snp included at VarScan v2.4.2 ([@evy009-B53]).  We considered only SNPs found in at least two libraries, regardless of sex, for further analysis. 

Allelic frequency for each SNP was calculated as the average of the frequencies estimated by VarScan in each library.  Hence, the frequency of each SNP was estimated based on 20--40 individuals depending on the number of samples that detected a particular SNP.  We determined whether SNPs promoted synonymous or nonsynonymous amino acid changes using the prediction of complete CDSs and a custom python script.  Allele frequencies were used to calculate the index of interspecific differentiation (*D*) defined as the absolute value of the difference in allele frequencies of an SNP in *A. fraterculus* and *A. obliqua* ([@evy009-B83]; [@evy009-B2]).  The statistical comparison of *D* distributions inferred for each type of SNP (synonymous, nonsynonymous, and noncoding) was performed by applying Kolmogorov--Smirnov tests.  We also estimated the average *D* using all SNPs in each CDS ($\overset{-}{D}$~CDS~), using only synonymous SNPs ($\overset{-}{D}$~S~), and using only nonsynonymous SNPs ($\overset{-}{D}$~NS~). 

McDonald--Kreitman test ([@evy009-B65]) (MKT) was performed for each CDS by comparing the number of synonymous and nonsynonymous substitutions of polymorphic and almost fixed SNPs (*D* \> 0.95) using a custom python script.  We removed variants with a frequency smaller than 0.05 in both species to avoid biases produced by segregation of slightly deleterious mutations ([@evy009-B76]).  Only genes that had a value of at least one in all four classes of SNPs were included in the analysis, and significant departures from neutrality were estimated by Fisher's exact test (two-tail *P*-value \< 0.05).  CDSs with significant Fisher's exact test and neutrality index (NI) lower than 1 were considered to evolve under positive selection.  The script also calculates the direction of selection (DoS), where a signature of positive selection is observed when DoS \> 0 ([@evy009-B99]). 

Calculating Evolutionary Rates
------------------------------

Complete CDSs were submitted to the reciprocal best hit strategy in BLASTn with an *e*-value threshold of 10^−6^ to obtain the potential pairs of ortholog CDSs between *A. fraterculus* and *A. obliqua*.  This strategy seeks to obtain the pairs of genes that produce best hit scores in a bidirectional BLAST comparison (interchanging the CDSs of the species as query and subject).  Pairs of sequences that showed a length difference greater than 5% were removed because there was a higher chance of being different isoforms or different genes with only similar domains.  We aligned the DNA sequences of putatively orthologs from the two species by their amino acid translations using the MAFFT algorithm ([@evy009-B52]) and back converted to DNA implemented in the program TranslatorX ([@evy009-B1]).  Resulted alignments were submitted to KaKs_Calculator ([@evy009-B121]) to calculate the pairwise nonsynonymous (Ka) to synonymous substitution rate (Ks) ratio of the *fraterculus* group lineage using the Model Selection framework ([@evy009-B79]).  To decrease the chance of poorly alignments or saturation, we removed pairs with outlier Ks values, defined as values greater than the average plus three times the standard deviation, which was 0.62.  Moreover, all the alignments with Ka/Ks \> 1 were visually checked. 

We also calculated the evolutionary rates of ortholog genes in Tephritidae and tested for selection using a phylogenetic approach.  For that, the CDSs of *Ceratitis capitata* (GCF_000347755.2) ([@evy009-B74]), *Rhagoletis zephyria* (GCF_001687245.1), *Zeugodacus cucurbitae* (GCF_000806345.1) ([@evy009-B93]), *Bactrocera dorsalis* (GCF_000789215.1), and *Bactrocera oleae* (GCF_001188975.1) were downloaded from GenBank.  To avoid using miss-annotated and miss-assembled sequences, we removed CDSs with more than one stop codon and reduced the redundancy using Cd-hit-est ([@evy009-B29]) with a similarity threshold of 0.99.  The putative cluster of orthologs was predicted using reciprocal best hit strategy in BLASTn with an *e*-value threshold of 10^−6^ and the CDSs of *A. fraterculus* as reference.  The complete clusters (seven sequences) were submitted to the filtering and alignments steps of POTION program ([@evy009-B41]).  This pipeline excludes the sequences based on relative sequence length and identity, then aligns the clusters, trims the alignments using trimAl v.1.2 ([@evy009-B11]), and detects recombination using three methodologies (Phi, NSS, and MaxChi2) implemented in PhiPack ([@evy009-B8]).  All parameters used in POTION are available in [supplementary file 1](#sup1){ref-type="supplementary-material"}, [Supplementary Material](#sup1){ref-type="supplementary-material"} online.  The maximum-likelihood phylogenies were inferred for each remained complete cluster of orthologs using GTRCAT model and 200 bootstrap replicates in the program RAxML v.8.2.9 ([@evy009-B98]). 

We used trimmed alignments and the phylogenies to estimate the global nonsynonymous/synonymous rate (d*N*/d*S*) ratio (ω) and performed the strict branch-site test implemented by CODEML included in the PAML v. 4.9 package ([@evy009-B115]).  The ω for the Tephritidae lineage was estimated using the M0 model (model = 0).  We removed clusters of orthologs with d*S* higher than the average d*S* plus three times the standard deviation (d*S* \> 7).  The ancestral branch of *A. fraterculus* and *A. obliqua* in each phylogeny was set as foreground for the branch-site test.  In order to statistically test whether a gene is evolving under positive selection, we compared the likelihoods of MA (model = 2, NSsites = 2) and MA1 (model = 2, NSsites= 2, fix_omega = 1) models using likelihood-ratio tests (LRTs) ([@evy009-B117]).  After the LRTs, we used the χ2 distribution to obtain *P*-values.  We also detected variation in selection strength across the cluster of orthologs using RELAX ([@evy009-B108]).  This program compares LRT of alternative and null (*k = 1*) models, where *k* is selection intensity defined as ω~foreground~ = ω~background~^*k*^. Significant comparisons with *k *\> 1 and *k *\< 1 indicate selection intensification and relaxation, respectively.  To perform the phylogenies and selection tests in parallel, we used a custom python script which uses the script raxml_bs_wrapper.py ([@evy009-B114]) and functions of the ete3 module ([@evy009-B44]). 

Comparing Sexes and Tissues
---------------------------

As *A. fraterculus* and *A. obliqua* are phylogenetically closely related and displayed similar patterns of gene expression, we compared the patterns of sequence evolution of sex-biased genes found in one and both species using ortholog information.  This approach allows the evaluation of genes with ancestral expression control and generalizes the results for both species and perhaps to other related species in the *fraterculus* group as well.  Besides, the genes with species-specific expression enabled the analysis of recent evolutionary patterns after the change in expression pattern.  Statistical comparisons of differentiation indexes and evolutionary rates among biased (sex and tissue) and unbiased categories were performed using Wilcoxon rank-sum test corrected by Holm approach ([@evy009-B40]).  Furthermore, differences between proportions of genes evolving under different selective regimes in biased and unbiased categories were estimated using Fisher's exact test. 

Results
=======

Sequencing, Cleaning, and Assessment of the Quality of the Assemblies
---------------------------------------------------------------------

We produced 28,808,966 × 2 reads for *A. fraterculus* and 28,020,776 × 2 reads for *A. obliqua*, from males and females, two replicates each totalizing eight RNA libraries of reproductive tissue (approximately 7M × 2 reads per library).  These libraries along with the previously sequenced samples of head tissue totalized 58,551,775 × 2 reads and 55,626,431 × 2 for *A. fraterculus* and *A. obliqua*, respectively.  The cleaning step removed an average of 13.94% of reads for *A. fraterculus* and 14.71% for *A. obliqua*.  Summary statistics for the two assemblies produced similar N50, mean, median, and length distributions ([supplementary table S1](#sup1){ref-type="supplementary-material"}, [Supplementary Material](#sup1){ref-type="supplementary-material"} online).  Furthermore, from 1,066 conserved Arthropoda ortholog groups, BUSCO identified 95% as complete orthologs and around 3% as fragmented orthologs in the transcriptome assemblies of *A. fraterculus* and *A. obliqua* ([supplementary table S2](#sup1){ref-type="supplementary-material"}, [Supplementary Material](#sup1){ref-type="supplementary-material"} online).  Moreover, roughly 50% of the complete orthologs were duplicated in the raw assembly; however, the redundancy in filtered unigenes was almost zero (approximately 0.5%). 

Functional Annotation
---------------------

Around 70% of the CDSs were annotated using *D. melanogaster* protein database and over 90% were matched with a protein from the GenBank nr protein including only Arthropoda entries ([supplementary table S3](#sup1){ref-type="supplementary-material"}, [Supplementary Material](#sup1){ref-type="supplementary-material"} online).  The comparison against the nr database showed that most frequent top hits were to Tephritidae species ([supplementary fig.  S1](#sup1){ref-type="supplementary-material"}, [Supplementary Material](#sup1){ref-type="supplementary-material"} online).  Blast2GO successfully mapped approximately 67% of the CDSs.  The distributions of the level 2 GO terms of both species and tissues were similar ([supplementary fig.  S2](#sup1){ref-type="supplementary-material"}, [Supplementary Material](#sup1){ref-type="supplementary-material"} online).  KOG functional classification also showed a similar representation of the categories in reproductive tissues of the two species ([supplementary fig.  S3](#sup1){ref-type="supplementary-material"}, [Supplementary Material](#sup1){ref-type="supplementary-material"} online). 

SNP Calling, Differentiation Indexes, and MKT
---------------------------------------------

A total of 226,827 and 140,504 intra- and interspecific SNPs were identified in reproductive and head transcriptomes, respectively.  We found 109,828 SNPs (79,947 of them associated with synonymous and 29,881 with nonsynonymous changes) in 3,662 coding regions expressed in reproductive tissues and 63,489 SNPs in 2,602 CDSs expressed in head tissues, of which 48,288 were synonymous and 15,201 nonsynonymous.  SNP frequency distribution showed that more than 50% of the SNPs have rare alleles ([supplementary fig.  S4](#sup1){ref-type="supplementary-material"}, [Supplementary Material](#sup1){ref-type="supplementary-material"} online).  We rejected the hypothesis that frequency distributions of *D* for synonymous, nonsynonymous, and noncoding SNPs were drawn from the same distribution (Kolmogorov--Smirnov test, *P* \< 0.01 for each of the three pairwise comparisons).  These distributions also revealed that over 6% of the total of SNPs were fixed or almost fixed between species (*D* \> 0.95) and the proportion of this type of SNPs is greater in nonsynonymous variants in both tissues ([fig.  1*A* and *B*](#evy009-F1){ref-type="fig"}).  The histograms of $\overset{-}{D}$~CDS~, $\overset{-}{D}$~S~, and $\overset{-}{D}$~NS~ showed that there is a greater proportion of highly differentiated unigenes using nonsynonymous than synonymous SNPs ([fig.  1*C* and *D*](#evy009-F1){ref-type="fig"}).  After excluding rare alleles, we retained 906 CDSs which met the minimum requisites to perform the MKT, that is, at least one synonymous and one nonsynonymous fixed SNP and one synonymous and one nonsynonymous polymorphic SNP.  Fifty-one CDSs showed significant statistical departure from neutrality and NI \< 1, thus were considered evolving under positive selection ([supplementary file 2](#sup1){ref-type="supplementary-material"}, [Supplementary Material](#sup1){ref-type="supplementary-material"} online). 

![ ---Frequency distributions of differentiation index. (* A*) and (*B*) Distributions of *D* (absolute allele frequency differences between *A. fraterculus* and *A. obliqua*) of SNPs found in reproductive and head transcriptomes, respectively.  Light blue, yellow, and blue bars represent the distribution of noncoding, synonymous, and nonsynonymous SNPs, respectively. (* C*) and (*D*) Distributions of average *D* per CDS using all SNPs ($\overset{-}{D}$~CDS~), using only synonymous ($\overset{-}{D}$~S~) and only nonsynonymous ($\overset{-}{D}$~NS~) found in reproductive and head transcriptomes, respectively.  Gray, yellow, and blue bars represent $\overset{-}{D}$~CDS~, $\overset{-}{D}$~S,~ and $\overset{-}{D}$~NS,~ respectively.](evy009f1){#evy009-F1}

Patterns of Gene Expression across Sexes and Tissues
----------------------------------------------------

There were 12,887 and 13,605 unigenes expressed (TPM \> 1) in the *A. fraterculus* transcriptome in head and reproductive tissues, respectively, with similar values also found in *A. obliqua*: 12,073 (head tissue) and 13,455 (reproductive tissue).  Biological replicates are strongly correlated, with coefficients ranging from 0.96 to 0.98 and 0.95 to 0.98 for *A. fraterculus* and *A. obliqua*, respectively ([fig.  2*A* and *B*](#evy009-F2){ref-type="fig"}).  Moreover, female and male samples of both species showed Pearson correlation coefficients higher than 0.96, establishing well-defined clusters in the principal component analysis ([fig.  2*C* and *D*](#evy009-F2){ref-type="fig"}).  Patterns of differential gene expression across sexes and tissues were similar in both studied species ([fig.  3](#evy009-F3){ref-type="fig"}).  A total of 21.3% and 28.7% of the unigenes showed biased expression between sexes in *A. fraterculus* and *A. obliqua*, respectively ([supplementary table S4](#sup1){ref-type="supplementary-material"}, [Supplementary Material](#sup1){ref-type="supplementary-material"} online).  This difference is mainly affected by 6% more genes which are upregulated in *A. obliqua* males.  Interestingly, less than 1% of the genes in head transcriptomes are sex-biased in both species ([supplementary table S4](#sup1){ref-type="supplementary-material"} and [fig.  S5](#sup1){ref-type="supplementary-material"}, [Supplementary Material](#sup1){ref-type="supplementary-material"} online).  In the comparison across tissues, approximately 27% of the genes were tissue-biased in male transcriptomes of both species ([supplementary table S5](#sup1){ref-type="supplementary-material"}, [Supplementary Material](#sup1){ref-type="supplementary-material"} online).  Female transcriptomes showed 22% of tissue-biased genes in *A. fraterculus* and 33% in *A. obliqua*.  Besides, we also found variation in the magnitude of differential expression of biased expressed genes.  Male-biased genes displayed greater fold-change average (measure by log~2~) than female-biased genes in both species (Wilcoxon rank-sum test *P*-value \< 0.01; [supplementary fig.  S6](#sup1){ref-type="supplementary-material"}, [Supplementary Material](#sup1){ref-type="supplementary-material"} online).  Tissue-biased genes expressed in males and females showed opposite patterns of magnitude of gene expression, whereas in males, the genes with greater differences in gene expression are reproductive-biased, in females, they are head-biased. 

![ ---Analysis of expression of biological replicates.  Heatmap of Pearson correlations and hierarchical cluster of samples from *A. fraterculus* (*A*) and *A. obliqua* (*B*).  Principal component analysis of all samples from *A. fraterculus* (*C*) and *A. obliqua* (*D*).  RM and RF: Samples from male and female reproductive transcriptomes.  HM and HF: Samples from male and female head transcriptomes.](evy009f2){#evy009-F2}

![ ---Heatmap and hierarchical clustering of differentially expressed genes on head and reproductive transcriptomes from *A. fraterculus* and *A. obliqua*.  Differentially expressed genes between male and female reproductive transcriptomes from *A. fraterculus* (*A*) and *A. obliqua* (*D*).  Differentially expressed genes between male reproductive and head transcriptomes from *A. fraterculus* (*B*) and *A. obliqua* (*E*).  Differentially expressed genes between female reproductive and head transcriptomes from *A. fraterculus* (*C*) and *A. obliqua* (*F*).  RM and RF: Samples from male and female reproductive transcriptomes.  HM and HF: Samples from male and female head transcriptomes.](evy009f3){#evy009-F3}

Evolutionary Patterns of Sex- and Tissue-Biased Genes
-----------------------------------------------------

Most comparisons between population differentiation index averages ($\overset{-}{D}$~CDS~, $\overset{-}{D}$~NS~, and $\overset{-}{D}$~S~) across sex and tissues failed to show significant differences between biased and unbiased genes, but the few that did, involved contrasts to nonsynonymous mutations ([fig.  4](#evy009-F4){ref-type="fig"}).  Male-biased unigenes in both species displayed greater levels of differentiation than unbiased using the parameter $\overset{-}{D}$~NS~ (Wilcoxon rank-sum test *P*-value \< 0.05), whereas male reproductive-biased genes had the highest average $\overset{-}{D}$~NS~, which was also significantly different from male head-biased and unbiased genes (Wilcoxon rank-sum test *P*-value \< 0.01 in both comparisons).  In contrast, female transcriptomes failed to show significant differences in any comparison ([supplementary fig.  S7](#sup1){ref-type="supplementary-material"}, [Supplementary Material](#sup1){ref-type="supplementary-material"} online). 

![ ---Boxplots of differentiation indexes measured among sex-biased and unbiased genes.  Differentiation was estimated as average allele frequency differences between *A. fraterculus* and *A. obliqua* using all ($\overset{-}{D}$~CDS~), nonsynonymous ($\overset{-}{D}$~NS~) and synonymous ($\overset{-}{D}$~S~) SNPs.  Sex-biased genes are grouped in genes with the same expression pattern in both species (both spp) and biased expression detected in a particular species (sp-specific).  Comparison of $\overset{-}{D}$~CDS~ (*A*), $\overset{-}{D}$~NS~ (*B*), and $\overset{-}{D}$~S~ (*C*) among male-, female-biased and unbiased genes expressed in reproductive tissues.  \*Holm-corrected *P*-value of Wilcoxon rank sum test \< 0.05. \*  just above the box indicates significant level in the comparison to unbiased genes.](evy009f4){#evy009-F4}

Analysis of approximately 4,000 orthologs between *A. fraterculus* and *A. obliqua* and approximately 3,000 among seven Tephritidae species revealed that male-biased genes in both nonspecies-specific and species-specific groups have significantly higher evolutionary rates than unbiased genes ([fig.  5](#evy009-F5){ref-type="fig"} and [supplementary table S6](#sup1){ref-type="supplementary-material"}, [Supplementary Material](#sup1){ref-type="supplementary-material"} online).  Likewise, female-biased unigenes also showed significantly higher evolutionary rates than unbiased, though the comparison involving the species-specific expression genes failed to reject the null hypothesis ([fig.  5](#evy009-F5){ref-type="fig"}).  Moreover, comparisons between tissues revealed that reproductive-biased genes, be it male or female, displayed higher rates of evolution than unbiased in both species ([fig.  5](#evy009-F5){ref-type="fig"}). 

![ ---Evolutionary rates for sex- and tissue-biased genes estimated based on pairwise comparison (*A. fraterculus* and *A. obliqua*) and seven Tephritidae species (*A. fraterculus*, *A. obliqua*, *Ceratitis capitata*, *Rhagoletis zephyria*, *Zeugodacus cucurbitae*, *Bactrocera dorsalis*, and *Bactrocera oleae*).  Sex-biased genes are grouped in genes with the same expression pattern in both species (both spp) and biased expression detected in a particular species (sp-specific).  Boxplots of Log~10~(Ka/Ks) from *A. fraterculus* and *A. obliqua* orthologs for sex-biased genes (*A*) and tissue-biased genes (*B*) and (*C*).  Boxplots of Log~10~ (*d*N/*d*S) from seven Tephritidae species orthologs for sex-biased genes (*D*) and tissue-biased genes (*E*) and (*F*).  \*Holm-corrected *P*-value of Wilcoxon rank sum test \< 0.05.  \*\*Holm-corrected *P*-value of Wilcoxon rank sum test \< 0.01. \*  or \*\* just above the box indicates significant level in comparisons with unbiased genes.](evy009f5){#evy009-F5}

Male-biased and male reproductive-biased genes displayed significantly greater proportion of genes evolving under positive selection than unbiased as evaluated by MKT and pairwise Ka/Ks ([table 1](#evy009-T1){ref-type="table"}).  Moreover, these contrasts showed a positive mean DoS, suggesting adaptive evolution, even though these set of genes exhibited higher values than unbiased, this difference was not statistically significant ([supplementary fig.  S8](#sup1){ref-type="supplementary-material"}, [Supplementary Material](#sup1){ref-type="supplementary-material"} online).  These results diverge from what was found by the branch-site tests, which indicated similar proportion of genes evolving under positive selection among biased and unbiased genes.  Interestingly, we noticed that 11% of the male-biased genes with signals of positive selection play important roles in *Drosophila*\'s reproduction ([table 2](#evy009-T2){ref-type="table"}), and it is possible that they may retain similar roles in *Anastrepha*.  However, we also found a greater proportion of genes with signals of relaxed selection among male-biased and reproductive-biased (both sexes) genes in comparison to unbiased genes, suggesting differences in selective constraints in these contrasts ([table 1](#evy009-T1){ref-type="table"}).  Table 1Patterns of Evolution for Sex- and Tissue-Biased and Unbiased GenesMcDonald--Kreitman TestPairwise Ka/KsBranch-Site TestRELAX*N****P* \< 0.05**[^a^](#tblfn2){ref-type="table-fn"}*N***Ka/Ks \> 1**[^b^](#tblfn3){ref-type="table-fn"}*N****P* \< 0.05**[^c^](#tblfn4){ref-type="table-fn"}*N****P* \< 0.05**[^d^](#tblfn5){ref-type="table-fn"}Reproductive Male-biased (both)155**16**\*\*488**23**\*27217272**70**\*\* Male-biased (specific)15**4**\*\*282101841218422 Female-biased (both)22011547137115 Female-biased (specific)361195**11**\*1366136**33**\*\* Unbiased4071532749424811202481431Male Reproductive-biased174**19**\*527**28**\*\*30314303**84**\*\* Head-biased2502579161616126 Unbiased4021733219225061272506432Female Reproductive-biased403354131406140**36**\* Head-biased371234102171421729 Unbiased4141731298224061112406416[^3][^4][^5][^6][^7][^8]Table 2Signals of Positive Selection in Sex-Biased Expressed Genes Potentially Associated with ReproductionAnnotation with *Drosophila melanogaster* DatabaseExpression PatternSignal of SelectionRoleReference*Neural Lazarillo*Male-biased[^a^](#tblfn8){ref-type="table-fn"}Ka/Ks \> 1Fertility and courtship behavior[@evy009-B87]*Takeout*Male-biased[^c^](#tblfn10){ref-type="table-fn"}Ka/Ks \> 1Courtship[@evy009-B17]*CG15406*Male-biased[^b^](#tblfn9){ref-type="table-fn"}BstInfluence female's remating[@evy009-B95]*Kelch-like family member 10*Male-biased[^c^](#tblfn10){ref-type="table-fn"}MKTSpermatogenesis[@evy009-B3]*Dynein intermediate chain at 61B*Male-biased[^c^](#tblfn10){ref-type="table-fn"}BstSpermatogenesis[@evy009-B26]*Hedgehog*Male-biased[^b^](#tblfn9){ref-type="table-fn"}BstMale's germ line maintenance[@evy009-B122]*Male fertility factor kl5*Male-biased[^c^](#tblfn10){ref-type="table-fn"}MKTSperm motility[@evy009-B12]*Lost boys*Male-biased[^c^](#tblfn10){ref-type="table-fn"}MKTSperm motility[@evy009-B113]*Tektin A*Male-biased[^c^](#tblfn10){ref-type="table-fn"}MKTSperm motility[@evy009-B20]*Egg-derived tyrosine phosphatase*Female-biased[^c^](#tblfn10){ref-type="table-fn"}BstOogenesis and embryogenesis[@evy009-B110]*CG14645*Female-biased[^b^](#tblfn9){ref-type="table-fn"}Ka/Ks \> 1Courtship[@evy009-B46]*CG14187*Female-biased[^a^](#tblfn8){ref-type="table-fn"}Ka/Ks \> 1Chorion protein[@evy009-B104][^9][^10][^11][^12]

Discussion
==========

The RNA-seq data generated high-quality de novo assemblies of *A. fraterculus* and *A. obliqua* as evaluated by length distribution and gene content metrics.  The N50 values of around 1,800 bp of these assemblies were in line with equivalent transcriptomes of other available tephritids ([@evy009-B43]; [@evy009-B68]; [@evy009-B89]).  Additionally, we found almost all conserved Arthropoda ortholog clusters in these transcriptome assemblies, suggesting a significant representation and completeness for the panel of genes expressed in head and reproductive tissues of both species.  A further indication of their completeness is that most of the CDSs were successfully annotated against proteins of the *D. melanogaster* database (∼70%) and the nr GenBank database (∼90%).  Furthermore, functional annotation using the distribution of GO and KOG categories for reproductive tissues of *A. fraterculus* and *A. obliqua* showed similar distributions, akin to what has been described for head tissues from these species ([@evy009-B84]). 

Transcriptome data studied here enabled identification of thousands of SNPs across *A. fraterculus* and *A. obliqua*.  Even though we identified hundreds of SNPs fixed, or nearly fixed, in one or the other species, the most common pattern observed for *A. fraterculus* and *A. obliqua* transcriptomes indicates that the species have diverged recently, because a great number of SNPs show little allele frequency difference across species.  However, these results should be interpreted with caution, because we estimated SNP allele frequency distributions from pools of individuals of a single population per species, and they may not represent the whole diversity across the species\' geographic distributions.  Microsatellite analyses across *A. fraterculus* Brazilian populations showed some evidence of differentiation, but over 90% of variation is intrapopulational ([@evy009-B63]).  Furthermore, there is evidence that these species differentiated with gene flow ([@evy009-B92]; Díaz F, Lima ALA, Nakamura AM, Fernandes F, Sobrinho I, de Brito RA, unpublished data), which would make variation in general common to several localities, rather than isolated, even across species boundaries.  If this pattern is common across the species\' distribution, it might indicate that the diversity distribution here inferred for *A. fraterculus* and *A. obliqua* may hold for the majority of SNPs identified. 

The allele frequency distributions are consistent with a scenario where the majority of the genome would be somewhat homogenous, interspersed by highly differentiated regions ([@evy009-B64]), such as what was found in two recently diverged species of *Gryllus* in the presence of gene flow ([@evy009-B2]).  This pattern would hold even if there were no selective forces involved and only drift would be driving the species apart ([@evy009-B16]).  Here, despite the reduced number of contigs with large allele frequency differences across species, we still detected at least 5% of the SNPs nearly fixed for different alleles in different species.  Interestingly, the distributions of *D* (allele frequency difference between species) inferred for synonymous, nonsynonymous, and noncoding SNPs are significantly different from one another.  These differences seem to be at least in part due to adaptive evolution as there is a significantly greater proportion of fixed differences across species which are associated with nonsynonymous substitutions, even when contrasted with noncoding substitutions.  This increased proportion of nonsynonymous substitutions in fixed differences between species also holds when we consider all substitutions present in a CDS.  The unigenes with high $\overset{-}{D}$~NS~ values may potentially be "islands of divergence," which are genomic regions that remain differentiated between species even in the presence of gene flow due to directional selection ([@evy009-B71]).  Even though there are other reasons why genomic islands of divergence may occur ([@evy009-B69]; [@evy009-B16]), many which cannot be tested for the data presented here because of the lack of a reference genome, these results point that at least a portion of the divergence between *A. fraterculus* and *A. obliqua* is due to regions affected by selection.  We found some male-biased expressed genes with signals of positive selection involved with male courtship and fertility ([table 2](#evy009-T2){ref-type="table"}), thus these genes may be related to the establishment of prezygotic barriers.  This observation agrees with studies on morphotypes of *A. fraterculus* complex species which have suggested that their reproductive isolation is mainly due to prezygotic barriers ([@evy009-B88]; [@evy009-B51]). 

The SNP allele frequency distributions allowed us not only to identify genes potentially involved with species differences but also to investigate general patterns of evolution for genes expressed in reproductive tissues across the two closely related species.  In general, SNPs in both species showed a large proportion of rare alleles ([supplementary fig.  S3](#sup1){ref-type="supplementary-material"}, [Supplementary Material](#sup1){ref-type="supplementary-material"} online), which might be due to demographic expansion and/or selective sweeps or weak purifying selection in particular genes ([@evy009-B30]; [@evy009-B27]).  We consider the former to be more likely considering that this pattern seems to be widespread across several genes, and that population expansion due to the increased distribution of host fruits with agriculture has been suggested to have happened to both studied species based on coalescent simulations (Díaz F, Lima ALA, Nakamura AM, Fernandes F, Sobrinho I, de Brito RA, unpublished data). 

In general, the expression profiles were similar between *A. fraterculus* and *A. obliqua* showing 20--30% sex-biased unigenes, a pattern similar to what was found in *Drosophila* species ([@evy009-B119]).  This may be due to stability of the gene expression control machinery because of evolutionary constraints ([@evy009-B119]; [@evy009-B39]).  Data from both species consistently indicate that the majority of sex-biased genes comes from reproductive tissues, which agrees with what was reported for a comparison of somatic tissues and gonad transcriptomes in *D. melanogaster* ([@evy009-B75]).  Furthermore, the approximately 0.5% of sex-biased genes displayed in head tissues of both species of *Anastrepha* contrast with the approximately 16% differentially expressed genes in *D. melanogaster* head ([@evy009-B13]).  However, in the latter there is a large difference in the genes with sexually diverged expression between central system and peripheral tissues ([@evy009-B34]), so this variation in expression pattern across head organs and structures could obscure the expression of sex-biased genes in the whole head.  Besides, our results show a higher number of upregulated genes and larger magnitudes of their fold changes in males than in females, which could be due to the existence of more male-biased genes or the differences in expression of female-biased are too small, which would require larger statistical power to detect these differences ([@evy009-B4]).  The comparison between tissues reveals that most tissue-biased expressed genes are in male reproductive tissues, possibly because most tissue-specific genes are expressed in testis, as indicated for *D. melanogaster* ([@evy009-B66]). 

Our results suggest that male-biased genes have higher ω such as it was found for several lineages ([@evy009-B105]; [@evy009-B106]; [@evy009-B37]), such as *Drosophila* species ([@evy009-B120]).  In addition, female-biased genes also evolve significantly faster than unbiased genes in *Anastrepha* species, which has also been described for some animal taxa such as birds ([@evy009-B61]), mosquitoes ([@evy009-B73]; [@evy009-B109]), and fishes ([@evy009-B111]).  We also found that sex-biased genes also tend to evolve more rapidly, particularly their nonsynonymous rates, but the rate of fixation of nonsynonymous mutations ($\overset{-}{D}$~NS~) was only significantly greater in males.  These findings suggest that male-biased genes may have been evolving under adaptive constraints in the *Anastrepha* species studied here.  However, the faster evolution of sex-biased genes may be explained by other factors than sexual selection such as tissue-specific expression, genetic drift, turnovers in expression patterns, and relaxed selective constraints ([@evy009-B62]; [@evy009-B67]; [@evy009-B33]; [@evy009-B60]).  We found evidence of both selection and relaxed constraints in these genes.  Positive selection was detected by greater proportion of genes with significant MKT and high rates of evolution (Ka/Ks \> 1) in male-biased genes when compared with unbiased genes.  However, branch-site tests displayed a similar proportion of genes evolving under positive selection in male-biased and unbiased genes.  It is possible that this is a consequence of the reduced number of substitutions in the short branch between the two recently diverged species that failed to reach the significance level in the branch-site test, but this may also be caused by higher rates of evolution in male-biased genes which would complicate orthology assignment even for species in the same genus ([@evy009-B24]).  In fact, approximately 50% of male-biased genes only showed orthologs in *Anastrepha* but not to other more distantly related species, preventing their analysis in the branch-site test and thus potentially producing a bias in the proportion of positively selected genes. 

Our data also suggest signals of relaxed selective constraints, as male-biased and species-specific female-biased genes displayed a greater proportion of genes with signals of relaxed selection than unbiased genes, which would imply that the values are converging to neutrality ([@evy009-B108]).  As there are many more sites possibly under purifying selection than under positive selection, for instance, we detected only two male-biased genes that displayed Ka/Ks \> 1, whereas 68 had Ka/Ks \< 1, this pattern would imply that the average ω would tend to increase, helping to explain their higher rates of evolution.  On the other hand, genes with tissue-specific expression are probably evolving under relaxed selection.  Comparison between head and reproductive transcriptomes from *A. fraterculus* and *A. obliqua* revealed that reproductive genes evolve faster than head genes in both sexes, showing similar patterns to *Drosophila* ([@evy009-B47]).  This outcome may be explained by broader patterns of expression in head-biased genes, but as there are no available data for other tissues in these species, we are not able to estimate their actual specificity.  Nevertheless, testis typically shows a greater proportion of tissue-specific genes ([@evy009-B5]; [@evy009-B25]; [@evy009-B66]; [@evy009-B111]), thus it is likely that several reproductive-biased genes would be tissue-specific, particularly for males.  These potentially tissue-specific genes would be more likely to evolve under positive or relaxed selection, whereas generally expressed genes seem to experience stronger evolutionary constraints, possibly due to pleiotropy ([@evy009-B62]; [@evy009-B38]; [@evy009-B56]), or because genes that are expressed in different tissues across the organism are more likely to be part of the constitutive set of essential genes for cell function, thus harboring lower nonsynonymous rates.  Finally, gene duplication can reduce functional constraints, so distinct levels of paralogy may also produce this outcome.  Nevertheless, the mammalian gene family content is equivalent between housekeeping and tissue-specific genes ([@evy009-B118]); hence, if this pattern is also true in insects, gene redundancy may not be a plausible explanation. 

Our analysis identified hundreds of SNPs associated with unigenes that showed fixed, or nearly fixed, differences between *A. fraterculus* and *A. obliqua*, which have been significantly more associated with nonsynonymous substitutions than to other substitutions and point to an important role for selection in their differentiation.  Even though we still lack a formal connection between sex-biased genes and speciation, the first "speciation gene" identified in *Drosophila* was the male-biased *Odysseus site homeo-box*, associated with postzygotic isolation mechanisms which produce sterility in male hybrids ([@evy009-B103]).  Furthermore, if a set of sex-biased genes evolve under sexual selection or sexual conflict which may lead to intraspecific intersexual divergence, these same differences may foster differentiation between species, especially when the genes involved are sex-biased in a species-specific manner. 

Conclusions
===========

Our work not only contributes to the current functional genomic knowledge on two of the most important fruit pests from the Neotropics by generating next-generation transcriptome data for reproductive tissues which have been hitherto unavailable, but also explored differences in expression patterns between sexes and tissues.  Although several studies are available on this matter for a wide variety of animals, particularly *Drosophila* ([@evy009-B59]; [@evy009-B112]; [@evy009-B78]; [@evy009-B37]; [@evy009-B18]), little had been known for Tephritidae species.  In this matter, our findings indicate that head tissues of *A. fraterculus* and *A. obliqua* exhibit few genes with sex-biased expression.  More importantly, sexual dimorphism in expression profiles of reproductive tissues revealed that sex-biased genes evolve faster than unbiased genes, especially in males, a pattern that was associated with signals of positive selection and relaxed constraints.  Our results shed some light on the evolution of sex- and tissue-biased genes expressed in reproductive and head tissue of *A. fraterculus* and *A. obliqua* which should be valuable to other species as well.  Furthermore, we found a set of sex-biased genes in reproductive tissues that may be candidates to be involved in the differentiation process of *A. fraterculus* and *A. obliqua*.  However, further studies that evaluate the populational variation of these genes are necessary to corroborate their role in the differentiation of these and other species of the *fraterculus* group. 

Supplementary Material
======================

[Supplementary data](#sup1){ref-type="supplementary-material"} are available at *Genome Biology and Evolution* online. 

Supplementary Material
======================

###### 

Click here for additional data file. 

The authors would like to thank the Associate Editor Prof. Judith Mank and three anonymous reviewers for providing invaluable comments and methodological advice that considerably improved the manuscript.  This work was supported by the São Paulo Research Foundation (FAPESP, Brazil) (2010/20455-4 to R.A.B. and 2012/17160-8 to E.B.C.); National Council of Scientific and Technological Development (CNPq, Brazil) (160002/2013-3 to C.C. and 150626/2012 to S.C.-E.); and Coordination for the Improvement of Higher Education Personnel (CAPES, Brazil) (PVE 056/13 to R.A.B.). 

[^1]: Associate editor: Judith Mank

[^2]: **Data deposition:** Sequenced reads derived from reproductive tissues of *Anastrepha fraterculus* and *Anastrepha obliqua* have been deposited in the NCBI Sequence Read Archive under the accession SRP128443.  The unigenes of *A. fraterculus* and *A. obliqua* assemblies are available upon request. 

[^3]: *N*, number of unigenes; both, same expression pattern in *A. fraterculus* and *A. obliqua*; specific, sex-biased expressed gene in either *A. fraterculus* or *A. obliqua*. 

[^4]: Significant departure from nonsynonymous and synonymous proportion of polymorphic and fixed SNPs using Fisher's exact test and NI \< 1. 

[^5]: Number of orthologs of *A. fraterculus* and *A. obliqua*. 

[^6]: Number of orthologs with significant LRTs between MA and MA1 using the *A. fraterculus* and *A. obliqua* ancestral branch as foreground. 

[^7]: Number of orthologs with *k* \< 1 (relaxed selection) and significant LRTs between null (*k *= 1) and alternative using the *A. fraterculus* and *A. obliqua* ancestral branch as foreground. 

[^8]: Fisher's exact test comparing biased with unbiased genes showing \**P* \<0.05 and \*\**P* \<0.01. 

[^9]: [Note]{.smallcaps}.---Bst, branch-site test. 

[^10]: Sex-biased expressed gene only in *A. fraterculus*. 

[^11]: Sex-biased expressed gene only in *A. obliqua*. 

[^12]: Sex-biased expressed gene in both species.