The Evolution of Gene Expression Underlying Vision Loss in Cave Animals.

Dissecting the evolutionary genetic processes underlying eye reduction and vision loss in obligate cave-dwelling organisms has been a long-standing challenge in evolutionary biology. Independent vision loss events in related subterranean organisms can provide critical insight into these processes as well as into the nature of convergent loss of complex traits. Advances in evolutionary developmental biology have illuminated the significant role of heritable gene expression variation in the evolution of new forms. Here, we analyze gene expression variation in adult eye tissue across the freshwater crayfish, representing four independent vision-loss events in caves. Species and individual expression patterns cluster by eye function rather than phylogeny, suggesting convergence in transcriptome evolution in independently blind animals. However, this clustering is not greater than what is observed in surface species with conserved eye function after accounting for phylogenetic expectations. Modeling expression evolution suggests that there is a common increase in evolutionary rates in the blind lineages, consistent with a relaxation of selective constraint maintaining optimal expression levels. This is evidence for a repeated loss of expression constraint in the transcriptomes of blind animals and that convergence occurs via a similar trajectory through genetic drift.

Introduction

The loss of complex traits is a ubiquitous evolutionary pattern and the loss of visual traits in cave dwelling organisms is one of the most striking and conspicuous examples (Porter and Crandall 2003). Unlike fauna of other dark habitats, such as the deep sea, who more often retain functional visual systems and interact with light through vertical migration or bioluminescence (Warrant and Locket 2004), many cave-dwelling animals spend their entire lifecycles in the complete absence of light (Poulson and White 1969; Culver and Pipan 2009). In response to this environmental pressure, vision loss in caves has evolved in parallel across widely divergent taxonomic groups. Comparing blind cave-dwelling organisms to their sighted close relatives provides an opportunity to dissect the evolutionary genomic processes involved in trait loss and parallel evolution.

Several nonmutually exclusive theories have persisted regarding the evolutionary processes involved in vision loss in aphotic environments (Culver 1982; Culver et al. 1995; Culver and Wilkens 2000). The earliest of these theories was that in the absence of light, functionally useless eyes are no longer maintained by selection and the genetic factors controlling the development and maintenance of vision can accumulate deleterious mutations without being removed by purifying selection (table 1—Loss through drift) (Barr 1968; Poulson and White 1969; Wilkens 1971; Wilkens 1988). Alternatively, it has been hypothesized that the cave environment imposes selective pressures that directly or indirectly favor the loss of eyes and vision (Jeffery 2005; Yamamoto et al. 2009; Borowsky and Cohen 2013; Moran et al. 2015) (table 1—Loss through parallel selection, Loss through divergent selection). DNA sequencing has allowed relevant data collection to directly test these hypotheses. Yet, molecular studies have revealed different genes in different systems with varying signatures of selection and there has yet to be strong consensus on the relative importance of these processes in shaping vision loss and adaptation to an aphotic environment (Crandall and Hillis 1997; Culver and Wilkens 2000; Leys et al. 2005; Friedrich et al. 2011; Carlini et al. 2013; Hinaux et al. 2013; Klaus et al. 2013; Niemiller et al. 2013).

Table 1.

Hypotheses Considered for the Role of Gene Expression Variation in the Evolution of Vision Loss.

Hypothesis	Predictions—Whole Transcriptome Expression Phylogeny	Predictions—Individual Orthogroup Variance/ Evolutionary Rates	Relevant Parameter OU Process
Null—Constant stabilizing selection	Relative tree lengths among blind and sighted species equal that of molecular phylogeny	No significant shift / difference in strength of stabilizing selection through time	αblind = αsighted
	—	No significant shift / difference in strength of drift	σ2blind = σ2sighted
	—	Equal variance within and among species	βblind = βsighted
	—	No significant shift in optimum	θblind = θsighted
Loss through drift	Relative increase or decrease in tree length for blind species is less than that of sighted species	Decrease in strength of stabilizing selection through time	αblind < αsighted
		Increase in strength of drift	σ2blind > σ2sighted
		Increase in within-species variance relative to among species variance	βblind > βsighted
		No parallel shift in optimum	θblind = θsighted
Loss through parallel selection	Relative decrease in tree length among blind species which is greater than that among sighted species	Increase or no change in strength of stabilizing selection through time	αblind ≥ αsighted
	—	Decrease or no change in strength of drift	σ2blind ≤ σ2sighted
	—	Increase or no change in among-species variance relative to within-species variance	βblind ≤ βsighted
	—	Significant parallel shift in optimum	θblind ≠ θsighted
Loss through divergent selection	Relative increase in tree length among blind species which is greater than that among sighted species	Increase or no change in strength of stabilizing selection through time	αblind ≥ αsighted
	—	Decrease or no change in strength of drift	σ2blind ≤ σ2sighted
	—	Increase in among-species variance relative to within-species variance	βblind < βsighted
	—	No parallel shift in optimum	θblind = θsighted

Our understanding of the evolutionary processes acting on vision in cave systems is influenced by the degree to which molecular changes occur in parallel in independently evolved blind animals; that is, whether convergent vision loss is the result of selection, constraint or random processes (Losos 2011). Parallel phenotypic evolution in a similar environment can be the result of a combination of the same and different mutations, genes, and gene functions, including expression patterns (Manceau et al. 2010; Pankey et al. 2014). Such variation could arise in an ancestral population and be independently fixed in different lineages, or independently arise in different lineages (Stern 2013). Some of the most compelling evidence in support of these different hypotheses in cave organisms comes from genetic crosses of independently evolved troglomorphic populations of the same species. Crosses of such populations of the Astyanax cave fish result in offspring with larger and better developed eyes and visual cells than either of the parent populations (Wilkens and Strecker 2003), providing strong evidence that different mutations are responsible for eye reduction in this system. In contrast, Protas et al. (2006) found evidence for different mutations acting on the same gene resulting in pigmentation loss in cave dwelling Astyanax. In a crustacean system, Protas et al. (2011) found different genetic loci responsible for loss and reduction phenotypes for both eyes and pigmentation in a cave versus surface population QTL comparison. These lines of evidence lend support to the hypothesis that complex trait loss in cave organisms occurs through different sets of de novo mutations arising in independent lineages.

It is increasingly being revealed that variation in gene regulation, rather than mutational differences in protein coding regions, is responsible for a significant proportion of phenotypic variance among closely related organisms (King and Wilson 1975; Carroll 2005). Indeed, gene expression variation is involved in local adaptation (Fraser 2013; Xu et al. 2015) and is a considerable source of variation upon which natural selection can act (Gilad et al. 2006; Whitehead and Crawford 2006b). Vision and eye development genes are highly conserved across taxonomic groups and are highly pleiotropic, meaning that complete loss of these genes or loss of function mutations to protein-coding regions are unlikely to be fixed in a population (but see Yang et al. 2016). Altered expression patterns of large-effect genes such as the Pax and Hedgehog gene families early in development are clearly important factors involved in the reduction of eyes in cave animals (Yamamoto et al. 2004). However, the role of evolutionary variation in gene expression across the transcriptome is in need of further investigation.

For an eye to remain functional, stabilizing selection should maintain gene expression levels at an optimum with variation within species due to a combination of environmental and genetic factors (Gilad et al. 2006; Bedford and Hartl 2009). The majority of differences in expression patterns in homologous tissues among species living under similar conditions should be due to genetic drift alone, unless there are considerable differences in function (table 1—Null) (Lande 1976; Yang et al. 2017). Within the constraints of stabilizing selection, it is expected that divergences in gene expression among species increase with phylogenetic distance, and that genes with greater variation within species should exhibit greater variation among species (Lande 1976; Whitehead and Crawford 2006b; Brawand et al. 2011; Musser and Wagner 2015; Rohlfs and Nielsen 2015). If vision loss in cave dwelling animals is driven primarily by directional selection on gene expression levels, one would expect optimal expression levels in the eyes of blind species to shift relative to their sighted counterparts and a subset of genes to decrease in variance within species (table 1—Loss through selection). If the stabilizing selection on gene expression is relaxed, one would expect expression levels to increase in variance within species relative to among species and the rate of expression evolution to increase (table 1—Loss through drift). If there is a general trend of a relaxation of stabilizing selection across the eye transcriptome, one would expect a greater proportion of genes to have accumulated selectively neutral variation resulting in no significant convergence in expression among blind species (table 1—Loss through drift). A directional selection scenario could increase the probability that changes in expression levels of the same genes would be selected for and possibly increase the overall similarity in expression patterns among blind species (table 1—Loss through parallel selection). Alternatively, independent vision-loss events could be the result of selection on different sets of genes or selection in different directions, resulting in divergent expression patterns (table 1—Loss through divergent selection).

In this study, we take a comparative transcriptomics approach to assess the relative roles of positive and relaxed selection on gene expression resulting in the parallel loss of vision in cave-dwelling crayfish, a group with over 45 described blind cave-dwelling species in several phylogenetically independent groups (fig. 1; Stern et al. 2017). To do so, we tackle the following aims: 1) Assess the degree of expression convergence in the eyes of independently blind species both across the transcriptome and in individual gene families; 2) Model the strength of stabilizing selection and drift on expression levels in blind lineages versus their sighted relatives; and 3) Test alternative hypotheses concerning the mechanism(s) of loss of vision in cave crayfish.

. 1.

Molecular phylogeny of species and individuals sampled with RNA-seq, adapted from Stern et al. (2017) and heatmap of pairwise Pearson correlation coefficients between individuals.

Results

Transcriptomes Cluster by Eye Function, Not Phylogeny

Under the null expectation that gene expression differences among species are largely due to drift constrained by constant stabilizing selection, gene expression divergences in homologous eye tissue should match the phylogeny of the species (Gilad et al. 2006; Whitehead and Crawford 2006a; Brawand et al. 2011; Musser and Wagner 2015; Yang et al. 2017). Significant clustering or divergence among blind species relative to sighted ones, considering expectations from the molecular phylogeny, would support the “Loss through parallel selection” or “Loss through divergent selection” hypotheses, respectively (table 1). The neighbor-joining expression phylogeny based on Pearson correlations from 3,560 orthologous transcripts did not match the molecular phylogeny, rejecting the null expectation of constant stabilizing selection alone acting on expression levels in this group (figs. 1 and 2). Although the majority of the groupings in the expression phylogeny have low bootstrap support, the groupings with the strongest support are the reciprocally monophyletic clustering of the sighted and blind species (fig. 2). The cave-dwelling sighted species fall outside of the surface-dwelling sighted species, but still group strongly with the rest of the sighted species. Although this suggests that habitat has an effect on expression patterns, the signal of vision loss is stronger than that of habitat alone. Expression divergence as measured by proportion of total tree length was higher among the blind species than the sighted ones (Wilcoxon Rank-Sum Test P < 0.001). However, both the blind and sighted species grouped together more tightly than is expected by the molecular phylogeny and the decrease in expected tree length is greater among the sighted species than the blind species (fig. 1ilcoxon Rank-Sum Test P < 0.001). This suggests a weak convergence in expression patterns in the blind species that is not greater than expected from random processes as species with conserved eye function exhibit less expression divergence than the blind species do. This is in support of the “Loss Through Drift” hypothesis (table 1).

. 2.

Analyses based on neighbor-joining trees of pairwise Pearson correlation coefficients of expression levels. (a) Neighbor-joining tree based on Pearson distances using expression means for each species. Node support based on bootstrap (N = 10,000) replicates. (b) Box-plots of proportion of total tree length using bootstrap replicate expression phylogenies (N = 10,000). Horizontal lines indicate the proportions of total tree length for each group using the maximum-likelihood molecular phylogeny. Vertical bars highlight the median differences between the proportions of total tree length in the expression phylogeny and the molecular phylogeny. The decrease in tree length in the expression phylogeny relative to the molecular phylogeny is greater in the sighted species (Wilcoxon Rank-Sum Test).

Individual Expression Shifts Do Not Drive Clustering of the Blind Transcriptomes

To test hypotheses about convergence and selection on expression levels, we modeled the evolution of expression for each orthogroup as a quantitative trait evolving under an Ornstein–Uhlenbeck (OU) process (Hansen 1997). This approach has been found to adequately model gene expression evolution (Bedford and Hartl 2009; Rohlfs et al. 2014) and can account for stochastic changes to expression levels (σ2), often interpreted as genetic drift, as well as a “rubber band” force (α) attracting expression levels toward an optimum (θ), often interpreted as the strength of a constant stabilizing selection through time (Lande 1976; Hansen 1997; Beaulieu et al. 2012). Comparing per-orthogroup likelihoods of single expression optima models to two optima models revealed 93 orthogroups (2.6%) with a significant expression shift in blind lineages (FDR-adjusted P < 0.05; supplementary table S2, Supplementary Material online). Interestingly, removing these 93 orthogroups actually increased the relative clustering of blind species in the expression tree (Wilcoxon Rank-Sum Test, P < 0.001; supplementary fig. S1, Supplementary Material online), suggesting that these genes with “convergent” expression patterns are not driving overall similarity of blind transcriptomes.

Expression Patterns Are Consistent with the Loss Through Drift Hypothesis

To test if vision loss is associated with weaker stabilizing selection and stronger drift (table 1, “Loss through drift”), we tested for differences in the α and σ2 parameters of the OU model between blind and sighted lineages. Allowing the α parameter to vary, 1,433 orthogroups (40.3%) were estimated to have a significant shift in stabilizing selection with 93.6% of these having lower α estimates in the blind lineages. Estimates of α across the transcriptome were significantly lower in blind lineages than sighted lineages (Paired Wilcoxon Rank-Sum Test P < 0.001; table 2). Allowing the σ2 parameter to vary between blind and sighted revealed a similar, but less pronounced signal on a per-gene basis. 55 orthogroups were found to have a significant shift in drift (σ2) estimates with 56% of these having greater σ2 estimates in blind lineages. Across the transcriptome, the strength of drift was found to be significantly stronger in the blind lineages than the sighted lineages (Paired Wilcoxon Rank-Sum Test P < 0.001; table 2). Removing orthogroups with significant shifts in α or σ2 decreased the relative clustering of blind species, suggesting that this increase in variance and evolutionary rate is what is driving the transcriptome clustering based on pairwise distances (supplementary figs. S2 and S3, Supplementary Material online).

Table 2.

Median Parameter Estimates Across Orthogroups (N = 3560) Comparing Expression Evolution Rates and Variance in Blind and Sighted (or cave and surface) Lineages.

	All Species			Sighted Species Only			Cave Species Only
	Blind (N = 8)	Sighted (N = 6)	P Value	Cave (N = 2)	Surface (N = 4)	P Value	Blind (N = 8)	Sighted (N = 2)	P Value
α	0.513	0.675	<0.001	2.17	1.59	<0.001	0.846	1.57	<0.001
σ2	0.374	0.134	<0.001	0.056	0.033	0.212	0.233	0.050	<0.001
β	2.88	1.93	<0.001	—	—	—	—	—	—

P-values are from paired Wilcoxon Rank-Sum Tests.

Finally, we estimated expression variance within versus among species using the Expression Variance and Evolution (EVE) model (Rohlfs and Nielsen 2015) considering the 12 species with more than one individual sampled. This approach extends the above OU model by estimating the proportion of variance within versus among species in the β parameter. Assuming that the majority of gene expression levels are under constant stabilizing selection (table 1—Null), maintaining a linear relationship of variances within and among species (Lande 1976; Whitehead and Crawford 2006b), a decrease in β estimates from ancestral sighted to derived blind lineages would suggest directional selection with lineage-specific expression levels (table 1—“Loss through divergent selection”). Significant increases in β estimates from sighted to blind lineages would suggest a relaxation of stabilizing selection maintaining expression levels at an optimum within species (table 1—“Loss through drift”).

Estimates of this ratio of within to among species variance were significantly higher when considering blind species alone versus sighted species alone (Paired Wilcoxon Rank-Sum Test, P < 0.001; table 2) indicating that there is a general increase in within-versus-among species variance estimates in expression levels in the blind species relative to the sighted species. The estimated shared β across all orthogroups was 1.52 times higher in blind than sighted lineages. Together, these results indicate that vision loss results in an increase in expression variance both within and among species and this is likely due to an increase in the strength of drift as a result of decreased stabilizing selection on expression levels across the transcriptome relative to sighted species, supporting the “Loss through drift” hypothesis (table 1).

Increased Expression Variance Is Due to Evolved Vision Loss

To test whether this signal of decreased stabilizing selection and increased drift is due to vision loss or to cave dwelling alone, we capitalized on our sampling of two cave-dwelling sighted species and repeated the above analyses comparing sighted cave versus sighted surface species, as well as blind versus sighted species within the cave habitat (i.e., excluding the surface species). The above pattern held true when considering blind versus sighted cave-dwelling species, but not when considering sighted cave versus sighted surface species (table 2). In fact, the strength of stabilizing selection was found to be higher in sighted cave species than sighted surface species across the transcriptome (table 2). This strongly suggests that the decreased strength of stabilizing selection and increased strength of drift is a vision-loss signal and not an environmental signal.

Discussion

We present evidence that parallel evolution of vision loss does result in convergent transcriptomes, but not in the sense of many of the same sets of genes being similarly up or down regulated. Rather, the gene expression similarity in blind lineages (and divergence from sighted lineages) is the result of similar increases in expression variance, which may be the result of relaxation of selective constraint on expression levels transcriptome-wide. Together, these lines of evidence are most consistent with the “Loss through drift” hypothesis (table 1).

The molecular genetic basis of eye reduction and vision loss in cave organisms is a highly complex evolutionary process that may never be reduced to a universal answer across taxa (Protas and Jeffery 2012), that is, there may be many paths to becoming blind even among close relatives. Such an outcome has serious implications for our understanding of the genomics of parallel and convergent evolution, as well as the evolution of trait loss. Across several studied systems, it seems clear that different processes underlie trait reduction versus augmentation in subterranean ecosystems (Protas et al. 2006, 2011; Hinaux et al. 2013; Klaus et al. 2013; McGaugh et al. 2014). In this study, we analyze parallel loss of function of a homologous complex trait and although there is convergence in expression patterns, it is not as strong as the similarity among tissues with conserved function (fig. 2). This is contrasted by what has been found in nonhomologous gains of function, for example, the remarkably similar gene expression patterns in independently evolved cephalopod photophores (Pankey et al. 2014). Other studies have pointed to different genomic patterns and processes responsible for loss of vision in distantly related lineages (Meng et al. 2013; Protas and Jeffery 2012). Here, we show in a group of closely related organisms that while expression shifts in different genes likely underlie vision loss in different lineages, a common mechanism of increased rates of expression evolution through relaxation of constraint underlies the parallel loss of vision.

The arthropod eye is a complex organ with a history of modular evolution (Oakley 2003). Even among blind cave crayfish, individual eye tissues have degenerated in different ways and likely at different evolutionary time points (Hobbs et al. 1977). For example, some species of the genus Procambarus, who are hypothesized to have more recently invaded caves, retain rudimentary ommatidia (Mejia-Ortiz and Hartnoll 2005; Wilkens and Strecker 2017). More ancient cave species in the genus Orconectes lack ommatidia completely and have highly disorganized retina and other eye structures (Cooper et al. 2001). All tissues collected for this study were collected in the light in order minimize effects due to variable environmental conditions across species, which may have resulted in light-induced gene expression even in the blind species. Therefore, although the blind species sampled here are united by their nonfunctional eyes, the lack of strong transcriptome convergence may be due to differential degradation of eye structures or variable responses to light exposure.

The increase in estimated within-species variation, as well as the decrease in estimated stabilizing selection pulling expression levels toward an optimum provide multiple lines of evidence for a relaxation of selection on gene expression levels across the transcriptome associated with the loss of vision. Increased expression noise (i.e., expression variance in isogenic cells) has been found to provide a selective advantage in fluctuating environments (Zhang et al. 2009). Although our recovered pattern could be the result of increased environmental variability within or across caves compared with surface habitats, considering the relatively homogeneous and buffered environment of cave habitats, we posit that relaxation of constraint maintaining optimal expression levels for eye function is a likely explanation for our recovered patterns. Additionally, the lack of significant convergence in expression patterns in blind lineages lends further support for the importance of random processes influencing patterns of gene expression. This is consistent with experimental evidence that random mutations to cis-regulatory elements can both increase and decrease expression levels (Metzger et al. 2016). There is also experimental evidence that selection tends to favor mutations that decrease expression variance (Metzger et al. 2015). Therefore, in the absence of selection maintaining functional eyes in blind cave species, mutations increasing expression variance may be allowed to persist without fitness consequences.

The importance of positive selection driving vision loss in subterranean animals has been debated for many decades and our study adds a meaningful piece to the puzzle. Given that ∼2% of the orthogroups analyzed were found to be commonly differentially expressed across vision loss events, it is possible that mutations in regions regulating the expression of these genes are primarily responsible for the loss of vision, even though the similar expression of these genes do not drive the overall transcriptome clustering of blind species. Indeed, several of these orthogroups exhibiting convergent expression shifts have functions likely related to visual phototransduction and photoreceptor specification (supplementary table S2, Supplementary Material online) and some genes related to the phototransduction pathway have been found to have signatures of positive selection on both the molecular sequences and expression levels in this system (Stern and Crandall 2018). It is also possible that the initial genomic mutation(s) disrupting expression levels in the eyes of blind animals were driven to fixation by positive selection. Future studies analyzing genomic sequences, especially regulatory regions, will be critical in determining the exact mutations responsible for vision loss in this system. Nevertheless, we do observe a likely reduction of selective constraint across the transcriptome that is not necessarily observed when considering individual coding sequences (Crandall and Hillis 1997; Carlini et al. 2013).

Our study focused on comparing transcriptome-wide gene expression in the eyes of adult individuals. Future studies could compare expression levels both across species and across tissues or developmental time-points. Such approaches could pinpoint the developmental and tissue-specific components of gene expression, which are of course important factors in the evolution of vision loss in cave animals. Experimental studies may also be able to determine the different genomic mechanisms underlying vision loss in different crayfish lineages, whereas our phylogenetic approach can only detect significant commonalities. For the purposes of this study, we focused only on orthogroups that were found in all species in our analysis to reduce biases from missing data. This approach neglects consideration of orthogroups that may not be expressed at a detectable level in certain species due to biological or technical reasons or genes that have been lost through evolution. Additionally, we considered the expression level of each orthogroup as an evolving unit. Insight could be gleaned from considering the phylogenetic history of each gene family, which often does not match the species phylogeny and includes gene duplications and losses. The methods we used here consider the expression level of each orthogroup as an independently evolving trait, which neglects the covariance among traits. Multivariate phylogenetic comparative methods are not yet applicable to genome-wide datasets, especially when the number of traits significantly exceeds the number of species, and methods development in this area will be of great interest (Dunn et al. 2013; Adams and Collyer 2017).

Studies of convergent and parallel evolution are continually finding that different genomic changes can be responsible for the evolution of similar phenotypes (Elmer and Meyer 2011; Stern 2013; Storz 2016). In other words, even as evolution may find similar solutions to an environmental challenge, the underlying genomic architecture can be remarkably different and dependent on the genomic background on which new traits are built and reduced. Although convergent evolution is often taken as some of the clearest evidence for adaptation, convergence can occur due to evolutionary constraint or stochastic processes (Gould and Lewontin 1979; Gould 1989; Losos 2011). Our study supports this hypothesis, as the same sets of genes do not seem to be responsible for the evolution of vision loss in independent events, at least in terms of mean expression levels. Nonetheless, the pervasive signature of increased expression variation within species and relaxed stabilizing selection across these independent events does suggest that the same evolutionary genomic processes, rather than the same loci, can result in the parallel trait loss.

Materials and Methods

Study System, Sample Collection, and RNA Sequencing

Of the ∼450 crayfish species in the North American family Cambaridae, whose crown age is estimated to be at ∼60 Ma (Bracken-Grissom et al. 2014), approximately 45 are obligate cave-dwellers (Crandall and De Grave 2017; Stern et al. 2017). Blind cave crayfish exhibit evolutionarily fixed vision loss with highly degraded and disorganized eye structures (Crandall and Cronin 1997; Cooper et al. 2001). We targeted eight blind crayfish species with large, stable populations from four lineages with independent vision-loss events (fig. 1 and supplementary table S1, Supplementary Material online; Stern et al. 2017). Six sighted species with close phylogenetic relationships to the blind species were also collected. Two of the sighted species were collected from caves (C. nerterius and C. tenebrosus) and may be recent cave inhabitants (Finlay et al. 2006). These two species have pigmented, fully developed eyes that do exhibit some reduction in size. Specimens were not acclimated to common conditions before tissue dissection, but to minimize environmental variation, eye tissue was dissected in light environments (either in the field or the lab) and stored in RNAlater and kept on ice or at −20 °C until RNA extraction. RNA was extracted from both eyes of 1–3 individuals of each species using a Trizol, Qiagen RNeasy hybrid extraction protocol. Strand-specific RNA-seq libraries were prepared for each individual sample using the Nugen Ovation Universal RNA-seq kit with probes designed to deplete rRNA using consensus sequences of 18S and 28S sequence alignments of crayfish and related species from GenBank. Libraries were sequenced on a NextSeq 500 at the George Washington University, Milken School of Public Health, generating 75 bp paired-end reads of fragments with a mean size of 200 bp. For the two species with only one individual sampled, we generated additional single-end 75 bp reads in order to aid in transcriptome assembly. An average of 25.1 million paired-end reads were generated per library (supplementary table S1, Supplementary Material online; NCBI SRA BioProject PRJNA464169, accession SRR7124921-SRR7124956). Trinity assembly statistics are available in Supplementary Material (supplementary table S1, Supplementary Material online).

Transcriptome Assembly, Expression Estimation and Orthogroup Identification

Raw reads were processed using BBDuk in the BBtools package (https://jgi.doe.gov/data-and-tools/bbtools; last Accessed May 12, 2017) to filter adapter sequences, low complexity sequences and low quality (Q < 10) bases using a sliding window (ktrim = r k = 23 mink = 11 hdist = 1 qtrim = w trimq = 10 minlen = 36 entropy = 0.01 entropywindow = 50 entropyk = 5 tbo). To filter sequences that may have originated from microbial contamination prior to assembly, we mapped reads against a database of reference and representative bacterial, archeal, and fungal genomes from the NCBI RefSeq using Bowtie 2 and removed those that mapped concordantly (Langmead and Salzberg 2012). We also filtered reads that mapped to crayfish 28S and 18S rRNA sequences to remove those that were not depleted during library preparation. Clean reads for each species were pooled and assembled using Trinity v2.4 using the appropriate strand-specific option for these libraries (–SS_lib_type FR) (Grabherr et al. 2011). Open reading frames were predicted using Transdecoder 3.0 retaining amino acid sequences with significant blastp hits to the proteomes of Eurytemora affinis (Eyun et al. 2017), Hyallela azteca (Poyton et al. 2018), and Daphnia pulexColbourne et al. 2011 as well as HMMER v3.1b1 (http://hmmer.org; Accessed March 7, 2014) hits to the Pfam protein domain database (Haas et al. 2013; Finn et al. 2016).

To facilitate comparisons of expression levels across species, we identified putative “orthogroups,” that is, sets of orthologs and paralogs that have descended from a single copy sequence in the ancestor of the species being analyzed, using Orthofinder (Emms and Kelly 2015). Redundancy in each transcriptome was reduced by clustering amino-acid sequences with CD-HIT to 99.5% similarity (Fu et al. 2012). The longest protein sequence from each orthogroup was searched against the UniProtKB/Swiss-Prot protein database for approximate annotations. To minimize biases resulting from incomplete sampling and uneven sequencing depth, we analyzed orthogroups that contained at least one transcript from all 14 species, resulting in 3,560 orthogroups. These orthogroups contained an average of 47.89 sequences or 3.42 transcripts per species. 89.9% of these had hits with a bitscore above 40 in the UniProt/Swiss-Prot database.

Expression levels were estimated for each orthogroup using RSEM (Li and Dewey 2011) by mapping trimmed reads against the assembled transcripts for each species with Bowtie 2 (Langmead and Salzberg 2012). Where multiple transcripts were present for a species in an orthogroup, expected counts were summed across transcripts. This approach takes into account the presence of incompletely assembled transcripts that are included in an orthogroup and the fact that RSEM effectively distributes read counts across transcripts that share multi-mapping reads. Orthogroup expression levels were normalized for sequencing depth and transcript length by transforming to transcripts per million (TPM) (Wagner et al. 2012). Cross-sample normalization was accomplished by using the trimmed mean of M-values (TMM) method (Robinson and Oshlack 2010). TMM normalized TPM values were square root transformed before all comparative analyses to correct for heteroscedacity in expression levels (Musser and Wagner 2015). Rarefaction curves were used to confirm transcriptome saturation using the R package vegan (Dixon 2003).

Comparative Analyses

Neighbor-joining trees were calculated based on pairwise Pearson distances using the R package ape (Saitou and Nei 1987; Paradis 2004). 10,000 bootstrap replicates of the expression data matrix were used to assess node support. To assess the degrees of clustering among the blind and sighted individuals, we calculated the proportions of total branch lengths in the bootstrapped neighbor-joining trees and compared this to the expected relative branch lengths for each group using the molecular phylogeny of the species from Stern et al. (2017) using the R package picante (Kembel et al. 2010).

The R package OUwie was used to estimate relevant parameters of the Ornstein–Uhlenbeck process for each orthogroup (Beaulieu et al. 2012). A phylogeny of the 14 species used in this analysis was pruned from the time-calibrated maximum likelihood phylogeny from (Stern et al. 2017) which was estimated using a molecular data matrix with 466 taxa, ∼70% of the described diversity (Crandall and De Grave 2017). The root of the phylogeny was assumed to be in the “sighted” state. Within-species variance was included as measurement error. Significant differences in log-likelihoods were assessed using likelihood ratio tests (LRT) and assuming a chi-square distribution with one degree of freedom. P-values were corrected for multiple testing using the Benjamini–Hochberg procedure (Benjamini and Hochberg 1995). When considering differences in parameter estimates of the models allowing rates to vary across the tree, we used a parametric bootstrap to avoid relying on maximum likelihood point estimates. For each orthogroup, we used the maximum-likelihood parameter estimates to simulate 1,000 expression matrices using the R package OUwie (Beaulieu et al. 2012). These were used to reestimate parameters and generate confidence intervals for parameter estimates. Medians of these distributions were used to compare transcriptome-wide parameter estimates with a paired Wilcoxon Rank-Sum Test.

We made maximum-likelihood estimates of the parameters of the EVE model using a phylogeny with only the blind species and one with only the sighted species. These two trees have the same root age and number of species; therefore, parameter estimates can be compared without transformation. To test for a general difference in within-versus-among species variation across the transcriptome, we compared median β estimates using a parametric bootstrap of ML estimates in blind and sighted species using a paired Wilcoxon Rank-Sum Test.

Data and Code Availability

RNA-seq data has been deposited to the NCBI Short Read Archive under BioProject PRJNA464169, numbers SRR7124921-SRR7124956. Code to perform the analyses is available on https://github.com/TheDBStern/interspecific_rnaseq (v0.1.0 DOI: 10.5281/zenodo.1243204).

Supplementary Material

Supplementary data are available at Molecular Biology and Evolution online.

Click here for additional data file.