Does effective population size affect rates of molecular evolution: Mitochondrial data for host/parasite species pairs in bees suggests not

Abstract Adaptive evolutionary theory argues that organisms with larger effective population size (N e) should have higher rates of adaptive evolution and therefore greater capacity to win evolutionary arm races. However, in some certain cases, species with much smaller N e may be able to survive besides their opponents for an extensive evolutionary time. Neutral theory predicts that accelerated rates of molecular evolution in organisms with exceedingly small N e are due to the effects of genetic drift and fixation of slightly deleterious mutations. We test this prediction in two obligate social parasite species and their respective host species from the bee tribe Allodapini. The parasites (genus Inquilina) have been locked into tight coevolutionary arm races with their exclusive hosts (genus Exoneura) for ~15 million years, even though Inquilina exhibit N e that are an order of magnitude smaller than their host. In this study, we compared rates of molecular evolution between host and parasite using nonsynonymous to synonymous substitution rate ratios (dN/dS) of eleven mitochondrial protein‐coding genes sequenced from transcriptomes. Tests of selection on mitochondrial genes indicated no significant differences between host and parasite dN/dS, with evidence for purifying selection acting on all mitochondrial genes of host and parasite species. Several potential factors which could weaken the inverse relationship between N e and rate of molecular evolution are discussed.


| INTRODUC TI ON
There has been very substantial and sustained interest in how the evolutionary rates of species vary with selection pressure and effective population size (N e ): a function of the number of individuals that contribute alleles to each generation (Kimura & Ohta, 1971;Woolfit, 2009;Woolfit & Bromham, 2003). One particular aspect of this interest has focused on evolutionary "arms races" between species that are locked into competition, such as hosts and their parasites, but which differ in their capacities for evolutionary rates (Kimura & Ohta, 1971;Lanfear et al., 2014;Woolfit, 2009). Very generally, evolutionary rates are influenced by selection and genetic drift. Both selection and drift can be strongly influenced by N e .
Larger populations provide more opportunities for mutations to enter each generation and hence be subjected to selection, and this can equate to more favorable mutations being available for selection to promote, leading to high rates of molecular evolution. On the other hand, smaller populations are expected to be subjected to higher rates of genetic drift which should also lead to relatively high evolutionary rates for neutral mutations and even deleterious mutations. These two different considerations of population size lead to specific predictions (Woolfit, 2009;Woolfit & Bromham, 2003), which we now briefly discuss.
For populations or species with relatively small N e , Ohta's (1973Ohta's ( , 1992 nearly neutral theory predicts that most nonsynonymous nucleotide substitutions will fall into a "nearly neutral" state -wherein extremely deleterious mutations should still be removed by selection and highly favorable mutations should be rare. Under this scenario, the nonsynonymous substitution rate (d N ) for "slightly" deleterious mutations increases because of enhanced genetic drift while d N for the slightly advantageous mutations decreases because small population size will lower the number of favorable mutations entering a gene pool; and also, slightly advantageous mutations are likely to be lost due to the lowered efficacy of positive selection. Nevertheless, given that slightly detrimental mutations comprise a substantial proportion of mutations in lineages with small N e , previous studies (Eyre-Walker et al., 2002;James et al., 2016;Woolfit, 2009) have suggested that the ratio of non-synonymous to synonymous substitution rates (referred to as ω or rate of molecular evolution) should be greater in lineages with relatively small N e (Ohta, 1992).
The ratio ω (d N /d S ) can provide an indication of the mode and strength of selection acting on protein-coding genes (Nielsen, 2005;Yang & Bielawski, 2000), where ω = 1 signifies neutral evolution, ω ~ 0 is suggestive of strong selective constraints, ω < 1 indicates purifying selection and ω > 1 indicates strong positive selection.
Empirical measurements of ω can therefore provide insights into the history of selection relative to population sizes of particular lineages and focal genes (Wagner, 2002).

| Practical consequences of varying effective population sizes
The above issues regarding ω and N e can potentially become very problematic for species that are locked into co-evolutionary arms races with other species that have much larger effective population sizes, and consequently potentially greater rates of adaptive evolution (Shokri Bousjein et al., 2016, 2017. Such a situation could arise in obligate parasite-host associations (species dyads) if parasites have much smaller N e than their hosts, where we might expect hosts to have higher rates of adaptive evolution and parasites to accumulate slightly deleterious mutations more rapidly, as suggested for some species of socially parasitic inquiline bees (Shokri Bousjein et al., 2016).
Several studies have compared rates of molecular evolution in organisms assumed to have different effective population sizes, but they have not produced consistent results. Spradling et al. (2001) compared the rate of cytochrome b evolution in 21 rodent species and found an inverse relationship between effective population size and the rate of molecular evolution. Johnson and Seger (2001) showed an increase in evolutionary rates of island avian species, which are supposed to have smaller effective population sizes compared to those occurring on the mainland. Furthermore, Woolfit and Bromham (2003) argued that long-term reductions in population sizes of endosymbiotic microorganisms compared to their free-living relatives caused an increase in the nonsynonymous substitution rates of the 16S rRNA gene. Bromham and Leys (2005) conducted a comparative analysis on social parasites of bees, wasps, and ants, using concatenated mtDNA, nuclear and ribosomal genes, and found that social parasites tend to have faster rates of nonsynonymous substitution than their social hosts, which they attributed to the effect of smaller effective population of parasites on the rate of molecular evolution. In contrast, Erler et al. (2014) found similar rates of evolution among almost all defense-related genes (antimicrobial peptide genes) when comparing host and socially parasitic bumblebees. A study by Helbing and Lattorff (2016) revealed that three antiviral siRNA genes evolved faster in host bumblebees compared to their respective parasitic species. Furthermore, Fouks and Lattorff (2016) discovered similar protein evolutionary rates for nuclear gene EF-1α between social Bombus terrestris and its parasitic species B. vestalis.
To further examine this issue, we compared rates of mitochondrial evolution using two allodapine bee host species and their obligate social parasite bee species, otherwise known as allodapine inquilines. These inquilines very rarely infest more than 5% of host colonies (Smith et al., 2013;Smith & Schwarz, 2006a, 2009).
They are locked in tight co-evolutionary arms races with their host because they are obligate parasites and host-specific (Smith et al., 2013;Smith & Schwarz, 2009). These parasitic species spend their almost entire life cycle within the nest of the host species and have extreme adaptations to social parasitism, including strongly reduced mouth parts and pollen-collecting scopae. Given these morphological variations, they are completely dependent on their host's colony for brood rearing (Michener, 1965(Michener, , 1970b(Michener, , 1971a(Michener, , 1975(Michener, , 1983.
Allodapine host and parasite clades are mostly sibling lineages (Smith et al., 2013) and their life-history traits such as body size and generation times are similar (Michener, 1965;Smith & Schwarz, 2006a, 2006b, making this group a good model system for examining rates of evolution in both hosts and parasites. The Australian inquiline genus Inquilina Michener forms the sister clade to its host genus Exoneura Smith (Chenoweth & Schwarz, 2011;Smith et al., 2013). A previous study estimated the relative N e of Exoneura and its parasite Inquilina based on the mean number of reproductive host and parasite females per nest and found that N e of parasite species is at least an order of magnitude lower than its host (Shokri Bousjein et al., 2016). Despite this, previous phylogenetic analyses of allodapine parasites by Smith et al. (2013) showed that they have been able to persist for long periods of evolutionary time (about 15 million years ago from their initial divergence) and are presumed to have followed their hosts through multiple speciation events.
In this study, we examine whether the much smaller effective population sizes of allodapine parasite species in the genus Inquilina give rise to faster rates of molecular evolution (dN/dS) compared to their Exoneura host species; as theoretically predicted (Woolfit, 2009;Woolfit & Bromham, 2003). We target mitochondrial genes to examine this hypothesis because mitochondrial DNA is highly variable in many animal species because of its elevated mutation rate, which stands as ideal marker for the study of evolutionary events among relatively close phylogenetic species, such as allodapine hostparasite dyads (Galtier et al., 2009).
Given that we have focused on the rate of molecular evolution (ω) as measures of the efficiency of selection, our a priori expectation is that selection should be relaxed on mitochondrial genes of parasite species compared to their hosts because of the predicted enhanced effects of genetic drift on species with small N e (Ohta, 1973;Weber & Diggins, 1990;Woolfit, 2009). In other words, the ratio of non-synonymous (slightly deleterious substitutions) to synonymous substitution (neutral substitutions) rates of mtDNA genes (referred to d N /d S or ω or rate of molecular evolution) should be greater in inquilina species with much smaller N e compared to Exoneura species.
Sampling was undertaken in December 2013, from the Gembrook region in the Dandenong Ranges of Victoria, Australia.
Nests containing host and parasite species were collected from dead and fallen fronds of the tree fern, Cyathea australis, early in the morning when bees were not active ( Figure 1). All collected nests were immediately stored in insulated boxes on ice and transported back to Flinders University for nest dissection. Ethics approval was not required for this study.

| RNA preparation and highthroughput sequencing
Specimens were snap frozen on dry ice and the head and metasoma of each species were dissected on dry ice in the laboratory and immediately preserved in RNAlater ® (Weber et al., 2010) to prevent RNA degradation. In order to replicate sampling procedures, head and metasomal tissues were analyzed separately, and we pooled tissues

| Bioinformatics
2.3.1 | Quality control, assembly, and gene orthology Raw sequence data was quality controlled using FASTQC (Babraham Institute) and trimmed using CUTADAPT (Martin, 2011) to remove Nextra adapters, SMARTER PCR Primers, reads containing suspected poly-Adenine and poly-Thymine tails, low quality reads (phred scores <30), and sequences shorter than 25 bp. The Trinity platform was used for de novo assembly of transcript sequences using a default k-mer of 25 (Haas et al., 2013;Tierney et al., 2015).
The quality of assembled contigs was then assessed by mapping them back to raw reads using the BOWTIE 2 aligner and the quantity of proper-paired reads was calculated. Read coverage statistics were also retrieved using BEDTOOLS (v.2.22.0). We used two approaches to determine and verify orthologous genes, sequence similarity and phylogenetic reconstruction (Tekaia, 2016;Tierney et al., 2015).

Sequence similarity inference
We matched assembled transcripts using the BLASTX algorithm against seven NCBI reference species (mitochondrial proteincoding genes), with a 10 −6 E-value cut-off. Reference species used in this study included the African honey bee Apis mellifera scutel- Colletidae; Table S1). We then ran two reciprocal best similarity hit methods (tBLASTn and BLASTn), with a 10 −6 E-value threshold to verify putative orthologs. We considered stringent criteria to call reciprocal best hits as orthologs if we found: alignment length ≥50%; protein identity ≥30%; and nucleotide similarity ≥50% (Tommaso et al., 2011). However, a few best hit contigs didn't match to any of the reference species mtDNA genes due to their short length. We therefore performed BLAST2BLASTN alignment between best-hit contigs that were perfectly matched to reference species mtDNA genes (subject sequence) and un-matched best-hit contigs (query sequence); as per Tierney et al. (2015).

Phylogenetic reconstructions
For each species, a consensus sequence was generated in Geneious v.9.1.4 (Kearse et al., 2012) from best hit contigs of head and metasomal tissue using pairwise alignment, Geneious algorithm, setting global alignments with free end gaps for alignment type, 65% similarity for cost matrix, 9 for gap open penalty, 3 for gap extension penalty and allowing determine sequence directions automatically. In the case of heterozygote sites between best hit contigs of head and metasomal tissues, the IUPAC ambiguity codes were used. Multiple sequence alignments of consensus sequences and reference orthologous genes were created using Translation alignment and MAFFT/ Geneious algorithms. Amino acid sequences were translated using the standard invertebrate mitochondria code and examined visually for unusual features.
A matrix of nucleotide alignments was then analyzed using a Bayesian inference method implemented in MrBayes v.3.2 (Huelsenbeck & Ronquist, 2001). Two different approaches were used to reconstruct Bayesian gene trees. Gene trees for each mtDNA gene were reconstructed using both host/parasite focal species and the reference species. We also constructed a phylogenetic tree by concatenating all mtDNA genes.
For MCMC analyses of each individual mtDNA genes, we partitioned nucleotide sequences by codon positions and applied a GTR + I + Γ nucleotide substitution model (the least restrictive model was chosen to avoid potential errors related to incorrect a priori model selection) for each codon partition. Nucleotide sequences were partitioned by genes for MCMC analyses of concatenated mtDNA genes. Analyses were run two times, with each run comprising 10 million generations, using three heated chains and one cold chain and with variable rate permitted, and sampling every 1000th iteration.
Likelihood plots and standard deviation of split frequencies (below 0.01) were used to verify stationarity distribution and length of run.
Parameter trace files of each run were also examined in TRACER v.1.6 (Rambaut et al., 2014) where the first 10% of trees were discarded as burn-in FigTree v1.3.1 (Rambaut & Drummond, 2010) was then used to visualize the Bayesian Inference (BI) phylogenetic tree and Hylaeus dilatatus and Colletes gigas were used as outgroups to root the tree.

| Detection of selection pressure
Given that purifying selection is the most prevalent form of selection in essential genes (e.g., mtDNA genes) (Comas et al., 2010;Jordan et al., 2002;Monteiro et al., 2010), we assumed this was also true in our focal bees mtDNA genes unless we found signals of BUSTED assigns what percentages of the sites evolve with negative selection (ω < 1).

BUSTED splits branches into foreground and background
partitions and fits a codon model with three rate categories as ω 1 ≤ ω 2 ≤ 1 ≤ ω 3 for each partition. It then estimates the proportion of sites per partition belonging to each ω class. This model is referred to as unconstrained model or alternative model. If the unconstrained model shows evidence of positive selection (ω 3 > 1) and proportion of sites assigned to that category is non-zero, BUSTED fits the constrained model (null model) where positive selection on the foreground branches is disallowed by constraining ω 3 = 1. A likelihood ratio test is then used to compare the unconstrained model with the constrained one. If the null hypothesis is rejected, it indicates that at least one site, at least some of the time, is under positive selection on the foreground branches.

| Comparison of mitochondrial evolutionary rate
Hyphy v.2.5 was utilized to compare rates of molecular evolution of mtDNA genes between each host and its associated parasite using a portion of inferred concatenated Bayesian tree (containing two focal hosts and their respective parasite species) and the likelihood function created by Hyphy.
Evolutionary comparisons were carried out at two levels: 1. We first estimated branch-by Likelihood ratio tests (LRT) were then used to explore evidence of branch-by-branch rate heterogeneity if H 0 was rejected.
2. We then compared the rate of molecular evolution when codons within a branch were assumed to have different evolutionary rates. A hypothesis testing framework RELAX 3 (a part of the Hyphy software package) (Wertheim et al., 2015) was used to compare rates of molecular evolution between host and parasite species using the efficacy of selection and quantifying the ω ratio.  Figure S2, Treatment 1). For the second treatment, the rate of evolution was compared between the common ancestral lineages for each of the host and parasite clades. For this, stem lineage of host and parasite clades were considered as reference and test branches respectively (electronic supplementary material, Figure S2, Treatment 2). Each parasite (test branch) was also compared with its own host (reference branch) for the third and fourth treatments (electronic supplementary material, Figure S2, Treatments 3 and 4). These treatments were first examined on individual gene because the rate of evolution is not equal between different mtDNA genes (Rand, 2001). However, given that mtDNA genes are completely linked and therefore, are not independent replicates of evolution (Mitterboeck & Adamowicz, 2013), we concatenated all mitochondrial genes into a single alignment and repeated all treatments for concatenated mtDNA genes ( Figure S2).

| Transcriptome assembly statistics
In total, Illumina sequencing is generated between 27.3 and 37.7 million paired-end reads, of which 16.6-29.2 million remained after quality trimming. Total assembled transcripts varied from 28,372-82,102 per species, and contig N50 ranged from 552-1526 bp.
Overall, 30.4-54.3 million reads in the assembled files were aligned to the post-trimming transcripts and 52.83%-77.62% of reads were determined as proper pairs. Detailed summary statistics on reads and assembly can be found in Table S2.

| Phylogenetic tree inference
Gene trees based on individual genes were sometimes poorly supported ( Figure S1) which is not unexpected given the short sequences for some genes. We therefore present here a phylogenetic tree of concatenated all mtDNA genes (Figure 2). The concatenated mitochondrial sequences retained for analyses were 8948 bp long.

| Test of selection
Although tree wide global rate of selection (Table 3) and all ω values yielded from RELAX analyses ( cepted where no statistically significant evidence of positive selection (ω > 1) on any host and parasite mitochondrial genes was found (Table 2). In addition, values of point estimates of ω calculated by BUSTED revealed strong evidence of purifying selection acting on mtDNA genes (reported as ω 1 and ω 2 in Table 2).

| Estimating branch-by-branch variation in rates
Likelihood ratio tests for each gene resulted in support for global rates, which posits the same ω for all host and parasite branches where all codon positions are assumed to have equal rate of evolution (Table 3).

| Estimating codon-by-codon variation in rates
We continued our comparisons by allowing variation in rates across all codon positions of each branch based on inferred equal ω among host and parasite branches. Note: The mtDNA best hit contigs derived from BLASTx which met the criteria of tBLASTn and BLASTn (alignment length ≥50%; protein identity ≥30% and nucleotide similarity ≥50%) are highlighted. ND-4L best hit contigs of both host and one inquiline species didn't meet those criteria and obtained from BLAST2BLASTN alignment. Genebank accession numbers were provided below the size of each mtDNA genes.

TA B L E 1 Characteristics of the mitochondrial genes of focal species
In all treatments, we found no evidence that purifying selection of parasite species is less efficient than host species. In other words, our analyses did not support predictions that rates of molecular evolution in mitochondrial genes are higher in parasite species than their hosts (Table 4).

| DISCUSS ION
Our initial aim was to examine whether the evolution processes of mitochondrial genes are faster in allodapine obligate social parasite species which exhibit much smaller effective population sizes than their hosts. Because N e has been predicted to affect the pattern and rates of molecular evolution (Woolfit, 2009;Woolfit & Bromham, 2003), a priori we expected socially parasitic species have faster rates of molecular evolution than their hosts due to their greatly reduced effective population sizes.
ω values calculated by BUSTED analyses indicated that all host and parasite mtDNA genes have undergone purifying selection, which is the major mode of selection acting on mitochondrial genes (Castellana et al., 2011;Soares et al., 2013;Stewart et al., 2008).
However, we found no evidence that the efficiency of purifying selection is reduced in parasite lineages due to possible increased genetic drift associated with much smaller N e . In other words, the RELAX outcome provided no support to our hypothesis that the rate of molecular evolution of mtDNA genes is higher in parasite species than their hosts in any of the treatments.
Our finding here is inconsistent with some previous studies, which found inverse relationships between N e and rates of molecular evolution (Bromham & Leys, 2005;DeSalle & Templeton, 1988;Johnson & Seger, 2001;Spradling et al., 2001;Weinreich, 2001;Woolfit & Bromham, 2003;Wu & Li, 1985). This discrepancy might be partly due to some deficiencies in previous studies, but also to different approaches used for comparisons. For instance, Wu and F I G U R E 2 Mitochondrial genes phylogeny. A consensus tree (from MrBayes analyses) derived from 11 concatenated protein-coding mtDNA genes, with posterior probability node support. Parasites and hosts are indicated by colored branches: pink for hosts and blue for parasites. The lines to the right of the terminal branches link each host to its associated parasite Li (1985) showed that globin genes evolve at higher rates in rodents with smaller N e compared to primate lineages. However, each analysis in their study was based on a single comparison, which might have limited the breadth of their conclusions. In addition, Spradling et al. (2001) found differences in rates of Cyt-b among rodent species with different N e s, but for the rate comparison, they didn't use phylogenetically independent comparison methods (e.g., as developed by Felsenstein, 1985) which led to inflated type I and type II errors (Harvey & Rambaut, 1998). Johnson and Seger (2001) also confirmed an opposite relationship between N e and rates of molecular evolution using island and mainland bird species; however, they used a small and taxonomically restricted dataset for their comparison. Later, Wright et al. (2009) used a much larger and more varied dataset and found no significant difference between island and mainland species.
Although Woolfit and Bromham (2003) found increased rates of 16S rRNA evolution in endosymbiotic bacteria and fungi with small effective population sizes, they found that A+T base composition Note: Significant positive diversifying selection on each branch was tested using the likelihood ratio test by comparing a constrained model in which positive selection is disallowed on foreground branches and unconstrained model which allows positive selection at a proportion of sites on foreground and background branches. Accordingly, none of the mtDNA genes were found to be under significant positive selection. The statics for constrained model was not produced for mtDNA genes including CO1, CO2, CO3, Cyt-b, ND-1, ND-4L, and ND-5 because unconstrained model didn't find any evidence of positive selection (ω 3 > 1) and proportion of sites assigned to that category was also zero. For these genes constrained model is the same as the unconstrained model. The small-sample Akaike Information Criterion (AICc, Sugiura, 1978) was used to compare the goodness of fit of two models. ω 1 and ω 2 reported here are calculated under the accepted constrained model. Note: Tree wide global rate posits ω does not vary from branch to branch in the tree and provides a crude measure of the overall strength of selection acting on each mtDNA gene. Local model allows a separate ω in every branch of the tree. The likelihood ratio test was used to compare global model (null hypothesis) and local model (alternative hypothesis) on each mtDNA gene and resulted in a very strong (p ≥ 0.05 in all cases) support in favor of the global rate. DF = Degrees of freedom. Note: The RELAX test was used to examine whether selection is relaxed or intensified on a subset of test branches compared with a subset of reference branches in a predefined tree. The relaxation coefficient, k, is used to estimate selection intensity. In the null model, the selection intensity is constrained to 1 for all branches, whereas in the alternative model, k is allowed to differ between reference and test groups. Acceptance or rejection of the alternative model is tested using a likelihood-ratio test, but Akaike Information Criterion was also included as measures of fit of the null model and the alternative model. Four different treatments were examined on host and parasite branches for RELAX analyses which are as follows. Treatment 1: the combined parasite clade was selected as a test branch and the combined host clade was treated as a reference branch. Treatment 2: common ancestral lineage of parasite was tested against common ancestral lineage of host. Treatment 3: I. schwarzi as a test branch and E. robusta as a reference branch. Treatment 4: I. excavata as a test branch and E. angophorae as a reference branch. All treatments were also tested on concatenated mitochondrial genes. Two treatments were examined on ND-3 gene as it wasn't found in I. schwarzi. No significant differences in purifying selection efficiency were found between host and parasite species in any treatments. ω reported here is calculated under the accepted null model. of the same gene, as another indicator of increased genetic drift, is not significantly different between compared taxa. Furthermore, Bromham and Leys (2005) used hosts and parasites of bumble bees, allodapine bees and ants for their comparisons and found consistently higher rates of molecular evolution in parasite lineages than hosts. However, they used concatenated mtDNA, nuclear and ribosomal genes for almost all their comparisons, which may conflate quite different evolutionary processes (e.g., recombinant nuclear genes with non-recombinant mt genes). In the case of the allodapine bees in their study, they did not compare inquiline species with their corresponding hosts, while in our study we explicitly compared rates between each host and its specific parasite.

TA B L E 3 Tree wide global rate of selection versus local rates
Our finding in this study raises a critical question which is: why the comparisons failed to support the predicted differences in efficiency of positive and purifying selection (using RELAX method) between hosts and their social parasites? We now put forward several possibilities that might explain this outcome.
One possible interpretation is that N e of inquilines is not reduced enough to affect the rates of molecular evolution despite their comparatively lower N e compared to the host. For example, while it has been found that N e of allodapine inquilines is about an order of magnitude lower than their hosts based on incidences of parasitization (Shokri Bousjein et al., 2016), it is possible that its absolute N e is still large enough to allow effective purifying selection. In addition, N e could also be influenced by variation in the number of offspring per individual, sex ratio (Hedrick & Parker, 1997;Shokri Bousjein et al., 2017), reproductive skew among individuals, and colony productivity. These factors are connected to each other in a complex way in allodapine bees (Schwarz et al., 2007). For instance, reproductive skew can be affected by relatedness among nest mates (Harradine et al., 2012;Langer et al., 2004), however, skew is also linked with colony productivity (Schwarz, 1986), and both of these co-vary with sex allocation (Schwarz, 1994). Because of these complicating factors, a precise estimation of N e is problematic. Nevertheless, genetic methods such as using the d N /d S ratio are an attempt to obtain an accurate measure of N e ; however, these methods are impractical when there have been historical changes in population size or selection on coding genes (Gregory & Witt, 2008;Whitney & Garland, 2010).

An alternative issue concerns how great the effect of any change
in N e should be on the rates of molecular evolution; so, while it is clear that N e may potentially have considerable effects on the rates of evolution, estimation of the magnitude of that effect is not simple (Bromham & Leys, 2005;Woolfit, 2009). Previous studies (e.g., Bachtrog, 2008;Woolfit, 2009) (1977) suggested that the effect of a change in population size on the rate of evolution is expected to be quite large because a substantial proportion of mutations fall in the range from 1/N e L to 1/N e S. By contrast, Kimura (1979) suggested fewer mutations will be in this range, therefore the difference in evolutionary rates between lineages with different N e will also be low. However, Silander et al. (2007) argued that the distribution of selective coefficients of slightly deleterious mutations is dynamic between taxa and even within a species. Initially, slightly deleterious mutations were often taken into account to determine the difference in rates of evolution because they are more common in small population, while slightly advantageous mutations are assumed to be rare in such lineages (Woolfit, 2009). However, Charlesworth and Eyre-Walker (2007) found that slightly advantageous mutations that are under positive selection are also relatively common in small populations. This subsequently led Woolfit (2009)  Furthermore, in this study we targeted only mtDNA genes to examine the hypothesis and found no difference between evolutionary rates of mtDNA genes among host and parasite species in any of the treatments. It is possible that sole use of mtDNA genes might not be adequate to examine the hypothesis. This might be partly due to the fact that mitochondrial genes have been less likely involved in the adaptation process given their role in basic metabolic functions (respiration; Galtier et al., 2009). In addition, given that the mtDNA genes are maternally inherited and have reduced effective population size, they therefore do not represent the true genomic inheritance of organisms (Ballard & Whitlock, 2004). Our findings suggested that nuclear genes must be included in order to discover potential changes in selection pressure between the host and parasite species.
In conclusion, our analyses suggest that relatively small N e of parasite lineages has seemingly no strong impact on rates of molecular evolution of their protein-coding mtDNA genes compared to their hosts. This does not match the patterns found by Bromham and Leys (2005) who included Exoneura and Inquilina in their study and it also conflicts with some broader theoretical considerations. Our results indicate a need to consider a variety of theoretical bases for comparative rates of evolution for host/parasite relationships, and these may include the evolution of highstake arms races into more benign relationships such as symbiosis or mildly deleterious associations. At present, the evolution of more benign species interactions has not been explored in terms of rates of molecular evolution and we argue that our data calls for such an examination in future studies.

ACK N OWLED G M ENTS
We are grateful to Mohammad Javidkar for his help with field collections, to Kathy Saint for her help in lab experiments, to Dr. Terry Bertozzi for his advice on NGS data analyses and to eResearch SA for bioinformatic support. This research was supported by two Australian grants from the Holsworth Wildlife Research Endowment (3000005577) and Nature Foundation, SA (3000006475).

CO N FLI C T O F I NTE R E S T
The authors have no conflict of interest to declare. Writing -review & editing (equal).

DATA AVA I L A B I L I T Y S TAT E M E N T
All mtDNA gene sequences were deposited in GenBank (See accession numbers in Table 1).