We generated maximum likelihood trees for 109 genes from a metazoan Light Interaction Toolkit (LIT 1.0; Additional file 1: Table S1; Additional file 2). From our efforts, we noted that many LIT genes do not have orthologous relationships across Metazoa. For example, we find that LIT genes with similar functions in distantly related taxa (e.g. arthropods and vertebrates) are often paralogs, not orthologs, due to lineage-specific gene duplications. Although evidence suggests that orthologs tend to be more similar functionally than paralogs, this does not hold true in the case of all gene families 26 , 27 . Thus, we conclude that tree-building is a useful approach for inferring the function of certain LIT genes, yet we still urge researchers to use caution when assigning functions to uncharacterized genes in the absence of functional tests and based on phylogenetic relationships alone. Further, how one selects the sequences used to build a tree may have a significant influence on the inferences drawn about the ancestral function and patterns of diversification of gene families. Thus, we advocate for an objective, repeatable approach to choosing sequences for gene trees (such as the one we employ here), especially when using phylogenetic relationships between these sequences to make inferences about the functions of newly sequenced genes.

We generated 28 transcriptomes for light-interacting tissues from a range of cephalopod mollusks, arthropods, and cubozoan cnidarians (Table 1). We sequenced transcriptomes from these taxa because they lacked genomic or transcriptomic resources, but are well-suited for answering certain questions about the function, development, and evolution of eyes and other light-interacting structures. The sizes of our transcriptomes varied (Additional file 3: Table S2). For example, the number of total bases in our transcriptomes averaged 2,903,000 ± 1,185,000 (mean ± std deviation) with a low of 89,000 bp (for an eye from the cephalopod Vampyroteuthis infernalis) and a high of 5,185,000 bp (for a ventral eye from the horseshoe crab Limulus polyphemus). The number of isotigs per transcriptome ranged from 168 (V. infernalis eye dataset) to 5,447 (for tissue from the eyes and head of the isopod Asellus aquaticus) and the mean sizes of isotigs ranged from 530 bp (V. infernalis eye dataset) to 1,397 bp (L. polyphemus ventral eye dataset).

To rapidly identify potential LIT 1.0 genes in our genetic datasets, we implemented PIA in Galaxy. Unlike past efforts at using phylogenetic methods to annotate transcriptomes, PIA does not require the re-calculation of gene trees every time a new sequence is to be analyzed. The output from PIA is a gene tree or a set of gene trees in Newick format that can be viewed using existing tools in Galaxy, such as tab2trees from the Osiris package 11 . These trees include sequences identified from predicted protein databases associated with 29 fully-sequenced genomes (Additional file 4: Table S3 and Additional file 5: Supplementary References for Table S3). The trees also include sequences marked as either Landmarks or Queries (Figure 1). Landmarks are genes (usually from model organisms) whose functions and/or patterns of expression have been characterized relatively well. Sequences marked "LANDMARK1" – which are highlighted with red squares when the trees are viewed using our tab2trees tool – are well-characterized LIT genes. Sequences marked "LANDMARK2" are also well-characterized genes, but are those that have functions different than the LIT genes that we are seeking. For example, we included certain non-opsin GPCRs in our trees for opsins, but we labeled them as LANDMARK2s because we have evidence that they are not involved in the detection of light. Queries marked “QUERY” – which are noted by yellow circles in the output from the tab2trees tool – represent potential LIT genes that PIA has identified from a particular genetic dataset. Promising queries from a transcriptome (i.e. ones that may represent orthologs of LIT genes) will tend to fall on short branches in phylogenetic positions that are sensible given established relationships between species (for an example, see the query tagged “Gprk1 hit UN0029 ORF1” in Figure 1). The output from PIA may also include query sequences that are close relatives, but not orthologs, of LIT genes (for examples, see the queries tagged “Gprk1 hit UN1121 ORF1” and “Gprk1 hit UN2338 ORF1” in Figure 1). Thus, we urge users of PIA to inspect carefully where queries fall on their respective gene trees and to make inferences about function accordingly.

Building additional gene trees will allow researchers to use PIA to search for LIT genes that are not included in our initial list or to search for new sets of genes. For example, there has been much recent interest in the sets of genes that underlie the process of biomineralization in animals 28 , 29 . Extensive databases of these genes have been published for invertebrate taxa (e.g. mollusks) that are not traditional model systems 30 . By building trees for these sets of genes and applying our methods for PIA, researchers will be able to survey new transcriptomes rapidly for genes that may be involved in biomineralization.

Across our 28 newly-sequenced transcriptomes, PIA identified potential orthologs of 69 of the 109 genes included in LIT 1.0 (Figure 2). We recovered certain genes from our transcriptomes far more often than others. Genes from LIT 1.0 that were expressed in ten or more of our transcriptomes included several components of the rhabdomeric phototransduction pathway, which is employed by the photoreceptors found in the eyes of many invertebrates 3 , 15 . These components include: Arrestin (Arr), Gq alpha (Galpha49B), Gq beta (Gbeta76C), protein kinase C (inaC), r-opsin (ninaE), phospholipase C (norpA), and transient receptor potential protein (trp). Ten or more of our transcriptomes also contained the enzyme aminolevulinate synthase (Alas), a component of the heme synthesis pathway 31 , as well as an aldehyde dehydrogenase (Aldh) that is related to the Ω-crystallins expressed in the lenses of the camera eyes of cephalopods 32 , 33 and the mirror-based eyes of scallops 34 - 36 . Amino acid and nucleotide sequences for the potential LIT genes that we identified from our 28 transcriptomes, along with the corresponding gene trees, are available on our publicly-accessible web server (http://galaxy-dev.cnsi.ucsb.edu/pia/) under the Shared Data tab.

The transcriptomes that we generated for cephalopods contained between two (Vampyroteuthis infernalis) and sixteen (Euprymna scolopes) LIT genes (Figure 2). The majority of these genes represent components of the rhabdomeric phototransduction pathway, which is known to confer light-sensitivity to photoreceptors from the retinas of cephalopods 37 - 40 . These components include r-opsins, subunits of the hetero-trimeric Gq protein with which r-opsins interact, phospholipase C, and the ion channel TRP. We also recovered a number of lens crystallins, including relatives of the S- and Ω-crystallins indentified previously from the lenses of cephalopods 18 . Lastly, we found orthologs of LIT genes that are associated with two vision-related features that may be unique to the eyes of cephalopods. First, we found orthologs of the retinoid-binding protein RALBP, which is involved in regenerating the chromophores employed by the visual pigments of cephalopods 41 , 42 . Second, we found sequences that may represent reflectins, which are proteins that contribute to the biological mirrors that cephalopods use to camouflage their eyes 43 , 44 .

Our transcriptomes from arthropods contained between three (Procambarus franzi – a crayfish) and 28 (Ischnura ramburii – a damselfly) potential orthologs of genes from LIT 1.0 (Figure 2). As in cephalopods, many of the genes we identified in arthropods represent components of the rhabdomeric phototransduction pathway. We also identified a number of genes that are associated with the synthesis of pterins and ommochromes, types of pigment found previously in the compound eyes of certain arthropods 45 . Additionally, several transcriptomes contained genes (e.g. pinta, ninaG, ninaD, ninaB, and santamaria) related to those that help synthesize the chromophores employed by the visual pigments in the eyes of the fruit fly Drosophila melanogaster 46 , 47 .

Finally, our transcriptomes for the adult rhopalia and planula larvae of the cubozoan cnidarian Tripedalia cystophora contained thirteen and seven light-interacting genes, respectively (Figure 2). The majority of these genes are related to those associated with phototransduction in bilaterians. For example, we found a previously characterized opsin in our transcriptome for rhopalia from adult T. cystophora 48 , as well as a Gs alpha subunit that is associated with light-detection in other cnidarians 49 . We also found evidence of J-crystallins, which are lens crystallins unique to the camera-type eyes of cubozoans 50 .

It is important to note that the absence of a particular gene from a transcriptome is not necessarily informative. Even if a transcriptome is “complete”, it is only complete for a particular piece of tissue, from a particular animal, at a particular time. Thus, we have tried to draw general conclusions and points of future interest from the genes that we identified from our transcriptomes. We hope researchers will apply the approach that we have developed here to re-visit these light-interacting tissues and sequence more deeply and more broadly so that statistically meaningful comparisons of gene expression may be drawn between them.

We sequenced 28 transcriptomes from 20 invertebrate taxa that lack genomic resources, but are well-suited for answering questions about the function, development, and evolution of eyes and other light-interacting structures (Table 1). For example, we generated transcriptomes from RNA expressed by the eyes and skin of certain cephalopod mollusks (squid and octopus). These animals may have the most complex light-influenced behaviors of any invertebrate 51 , 52 , but it appears that the eyes of cephalopods tend to contain only a single spectral class of photoreceptor ( 53 ; though see 54 as an exception). Additional physiological complexity may be suggested by the results of high throughput sequencing. It is also possible that certain visually-influenced behaviors in cephalopods – such as dynamic camouflage – may be influenced by molecular components that are expressed outside of their eyes. For example, past work suggests that certain cephalopods express LIT genes in their light-producing photophores 55 and in certain dermal cells 56 .

We also sequenced transcriptomes for a range of arthropods. We chose to study stomatopods (mantis shrimp) because they have an unsurpassed ability to distinguish different aspects of light. Certain species are maximally sensitive to twelve distinct wavelength peaks and some species can identify both linearly and circularly polarized light 57 - 60 . Similarly, we chose to study odonates (damselflies and dragonflies) because they have physiologically complex eyes 61 and display a diversity of visually-influenced behaviors 62 - 64 . To study the degeneration of eyes in arthropods from subterranean environments, we examined certain species of isopods and crayfish in which closely related species or populations live either above or below ground. Specifically, we sequenced tissues from the eye-bearing, surface-dwelling isopod Caecidotea forbesi and its eyeless, cave-dwelling congeneric C. bricrenata . We also sequenced transcriptomes for different populations of the isopod Asellus aquaticus , which has a surface-dwelling form and multiple cave-dwelling populations with typical cave morphologies like degenerated eyes 65 , 66 . Likewise, we generated transcriptome data from a pair of surface ( Procambarus alleni ) and cave ( P. franzi ) freshwater crayfish. Crayfish have previously been the focus of molecular evolutionary studies of opsin in cave/surface comparisons 67 . To study the evolution of sexually dimorphic eyes, we generated a transcriptome for the RNA expressed by developing eyes from the ostracod Euphilomedes carcharodonta , a species in which males have compound eyes, but females do not 68 , 69 . Other species in this family of ostracods exhibit a similar, but independently evolved eye dimorphism, suggesting that these ostracods may be a promising system for the study of sex-specific convergent phenotypic evolution 70 .

Lastly, we sequenced transcriptomes for Tripedalia cystophora , a cubozoan cnidarian (box jellyfish). Cubozoans are the only cnidarians with camera-type eyes and, for that reason, have been the subject of numerous studies of visual neurobiology 71 - 74 , morphology 75 , 76 , and behavior 77 , 78 . Transcriptomic resources will aid these efforts. Further, as cnidarians, cubozoans may help us understand the evolutionary origins of the metazoan phototransduction cascade 79 - 81 .

We extracted RNA from our tissue samples using either the organic solvent TRIzol (Invitrogen) or the Nucleospin RNA XS kit (Macherey-Nagel), in both cases following manufacturer’s protocol. In cases where we used TRIzol, we removed trace DNA with the Ambion TURBO DNA-free kit (Invitrogen). In all cases, we quantified RNA yield with a Qubit Fluorometer (Invitrogen), following manufacturer’s protocol. To generate cDNA from RNA, we used the SMARTer cDNA synthesis kit (Clontech). To reduce sequencing artifacts due to poly-T tracts, we used modified 3’-primers for first-strand synthesis: 5’-AAG CAG TGG TAT CAA CGC AGA GTA CTTTTTTCTTTTTT-3’. For second strand synthesis, we used the protocol outlined in the SMARTer cDNA kits and a number of cycles determined by a series of optimization procedures. We then purified the amplified cDNA using one volume per sample of phenol:chloroform:isoamyl (25:24:1 v/v/v) and standard protocols. Finally, we sequenced cDNA using the Roche 454 platform. Here, we followed manufacturer’s instructions and employed partial runs with a manifold to separate samples. To assemble our transcriptomes, we used GS De novo Assembler v2.3 (“Newbler”; 454 Life Sciences/Roche Branford, CT USA) set to default threshold options, and using the –vt option to remove adapters. Following assembly, we used LUCY 82 , 83 to trim low-quality nucleotide reads and delete any assembled contigs below 100 bp in length. Next, we ran isotigs from Newbler through the program iAssembler 84 to combine redundant isotigs, then ran the resulting sequences through the program ‘Get ORFs’ 12 , 85 , ignoring any sequences less than 30 amino acids in length, to produce the predicted protein sequences that we used in our PIA analyses. Assembled sequences and ORFs for our 28 transcriptomes are available on the Bitbucket public repository (http://bitbucket.org/osiris_phylogenetics/pia) and on a publicly-accessible web server (http://galaxy-dev.cnsi.ucsb.edu/pia/).

We assembled the LIT 1.0 by reviewing past research into the molecular components that underlie the function and development of light-interacting structures in metazoans (Additional file 1: Table S1). Specifically, the LIT 1.0 contains molecular components of rhabdomeric- and ciliary- type phototransduction 1 , 15 , 16 , 86 , transcription factors involved in the specification and development of photoreceptors and eyes 3 , 20 , 87 , genes involved in the synthesis and regeneration of the chromophore retinal 46 , 47 , 88 , lens crystallins 18 , 19 , 36 , 50 , reflectins 44 , components of the circadian clock pathway 15 , and the enzymes that transport and produce pigments such as melanins 89 , 90 , pterins 91 , ommochromes 15 , 17 , and hemes 31 . Genes from LIT 1.0 are an appropriate test case for PIA because the specific functions and expression patterns of many of these genes are well-characterized. Also, certain fundamental aspects of light detection – such as opsin-based phototransduction – appear to involve molecular components that are conserved broadly across metazoan phyla. Thus, we can make well-informed inferences about the functions of new sequences from transcriptomes based on their phylogenetic relatedness to LIT genes that have been characterized previously.

After assembling our list of genes for LIT 1.0, we used functionally characterized exemplars of each of these genes (i.e. those from model systems such as fly or mouse; see Additional file 1: Table S1) and the blastp algorithm to search the predicted protein databases associated with 29 fully-sequenced genomes, including those from 24 metazoans, two choanoflagellates, and three fungi (see Additional file 1: Table S1 for search settings and Additional file 4: Table S3 for details on the predicted protein databases that we searched). After removing duplicate genes with Similar Sequence Remover 11 , we aligned BLAST hits for each gene using MAFFT 24 , 25 . We removed genes on long branches using Long Branch Remover 11 and built trees with RAxML assuming WAG as the protein model, and using 100 bootstrap pseudoreplicates followed by maximum likelihood search for the best tree 92 , 93 .

We used PIA to search our 28 new transcriptomes for potential orthologs of genes from LIT 1.0. First, we searched translated versions of our transcriptomes using blastp and the same queries that we used to identify sequences for our pre-calculated gene trees (Additional file 1: Table S1). We used stringent settings for blastp, specifically taking the top three hits that surpassed an e-value cut-off of 1e-20. Future users of PIA have the option of retaining different numbers of top hits and adjusting the e-value cut-off point as they see fit. Next, we used MAFFT to align the hits from our BLAST searches against the sequences that we used to calculate trees for our 109 LIT genes. PIA provides the option of aligning sequences using MUSCLE 23 , MAFFT 24 , 25 , or MAFFT-profile 94 . MAFFT-profile is the fastest of the three options because it does not re-align all sequences, but rather aligns the new sequences to an existing alignment. Finally, PIA uses an Evolutionary Placement Algorithm (implemented in RAxML; 26 , 27 ) to place the potential LIT genes on to our pre-calculated gene trees using Maximum Likelihood. Briefly, EPA places new genes on each branch of a pre-calculated phylogeny and calculates a likelihood score. The placement with the best score is retained. This is much faster than recalculating the entire gene tree.