BACKGROUND
whole transcriptome shotgun sequencing, or rna-seq, is a high-throughput, next-generation sequencing tool that can efficiently identify tens of thousands of functional genes in an organism or specific tissue at a given time  <cit> .
this deep sequencing makes it a more attractive tool than microarrays for organisms lacking reference genomes, facilitating de novo transcriptome assembly .
its high coverage is desirable when profiling transcripts in tissues of unknown function, enabling researchers to generate and/or test multiple hypotheses at once , and in organisms potentially harboring symbiotic organisms that may or may not produce the transcripts of interest , as rna-seq can simultaneously identify genes from microbes and their vectors/hosts.
such a combination of low genome resource availability and enigmatic physiology exists in the stick and leaf insects , or phasmids.
though common in the pre-molecular biology era through the laboratory stick insect, carausius morosus
<cit> , phasmid research today is relatively limited.
few phasmids are pests of agricultural crops  <cit> , though they reach plague-like abundances in temperate forests  <cit>  and c. morosus is an invasive pest in several countries  <cit> .
all life stages of all species within the order feed exclusively on leaves  <cit> .
this obligate folivory is relatively rare: among insects it is known only from leaf beetles , while more basal “herbivores” such as grasshoppers and crickets  will quite readily scavenge vertebrate meat, engage in cannibalism, and even hunt and kill other insects  <cit> .
thus phasmids are an ideal system for studying the evolution of herbivory in the lower polyneoptera.
folivorous organisms benefit greatly from plant cell wall degrading enzymes , a group that includes cellulases, hemicellulases, lignases, pectinases, and xylanases  <cit> .
once thought to be limited to microbes, endogenous  pcwde production has since been found throughout the animalia.
in particular, cellulase  genes from the glycoside hydrolase family  <dig>  are now believed to have existed in the ancestor of all metazoan life  <cit>  as opposed to having been repeatedly acquired from microbes via horizontal gene transfer, as is thought to be origin of gh <dig> and gh <dig> cellulases in beetles  <cit> .
among insects, endogenous cellulases have been found in lower and higher termites, cockroaches, crickets, beetles , a firebrat  <cit> , a springtail  <cit> , and, recently, the phasmids.
high cellulase activity in the anterior midguts of two phasmid species, eurycantha calcarata and entoria okinawaensis, was detected, the responsible proteins isolated, and the genes encoding them sequenced.
sequence homology and antigency against an insect cellulase anti-serum supported an endogenous, insectan origin for the enzymes  <cit> .
such a process is slow and predicated on the translation of pcwde genes into enzymes active against laboratory substrates like carboxymethylcellulose or crystalline cellulose, whose specificity and selectivity are imperfect  <cit> .
whether or not phasmids contain other pcwdes, such as the cellobiases  that convert the products of cellulase into glucose monomers  <cit>  or the pectinases  that hydrolyze pectin into galacturonic acid monomers  <cit> , was unknown.
presence of such active enzymes could explain the obligate folivory of the phasmatodea and be a key factor in the order’s evolution  <cit> , which is itself a puzzle as the sister order to the phasmatodea is highly debated  <cit> .
microbiological assays of the phasmid gut suggest digestion in the order is symbiont-independent  <cit> , so any phasmid pcwde genes are likely endogenous, yet whether the genes show homology to insect or microbe genes depends on their own evolutionary origin  <cit> .
complicating the issue is the relative lack of genetic resources for phasmids or their most closely-related orders: orthoptera  <cit> , embioptera  <cit> , and notoptera/xenonomia  <cit> .
lastly, even if phasmids have pcwde genes, their expression is not a guarantee, nor are they necessarily expressed in the gut region where they are most active.
the phasmid midgut is physically differentiated into two sections: an anterior midgut  marked by circular pleating and folding, and a posterior midgut  studded irregularly by hollow bulbs with filamentous tubules called the “appendices of the midgut” that open into the midgut lumen  <cit> .
the appendices may have an excretory or secretory function, and the possibility exists that they produce digestive enzymes that are carried forward into the amg via countercurrent flow  <cit> .
in the face of all these unknowns, next-generation sequencing is the best resource for answering questions of phasmid digestive physiology efficiently and effectively.
here we used de novo transcriptome assembly to identify the genes expressed in the midguts of six species of phasmatodea from four families, while greatly increasing the publicly available genetic resources for the order.
we also used rna-seq on one exemplar species, peruphasma schultei  to quantitatively compare transcript expression between the amg and pmg, and produced a genomic dna library from the symbiont-free phasmid brain to confirm that identified transcripts were encoded by the insect itself.
our main goal was to identify the production organ of the phasmatodea endogenous cellulase, while simultaneously creating an inventory of expressed pcwde and other digestive genes in phasmids and generating hypotheses on their evolutionary origins and the putative functions of the midgut sections.
this study serves as a necessary preliminary for more targeted molecular work.
more broadly, our transcriptomes are useful for evolutionary analyses of non-cellulase pcwdes in insects and identifying potential genes with biotechnological applications such as in processing biofuel feedstock or improving its rheology  <cit> .
methods
insects and microscopy
insects used were peruphasma schultei , sipyloidea sipylus , aretaon asperrimus , and extatosoma tiaratum, medauroidea extradentata, and ramulus artemis  cultured at room temperature in the bohart museum of entomology, university of california, davis.
phasmids were fed an ad libitum diet of privet  for peruphasma, eucalyptus for extatosoma, and rosa sp.
for the others.
library prep and sequencing
the rna-seq study of peruphasma schultei made use of three biological replicates for both the anterior and the posterior midguts .
for each replicate, the guts of five fed, surface-sterilized, adult, male and female phasmids were removed under sterile conditions and emptied of their contents in several washes of 70% ethanol.
then the anterior and posterior sections were separately pooled and homogenized in trizol® reagent.
rna was extracted according the trizol-plus protocol, which includes an on-column dnaase digestion step.
total rna quality  was checked with the bioanalyzer  <dig>  libraries were made using the illumina truseq v <dig> kit according to the manufacturer’s instructions.
hundred base pair paired-end sequencing was performed on the hiseq  <dig> and the raw data uploaded to the ncbi sra database .
for quality control, low quality bases and adapter contamination were removed with the fastx toolkit  <cit>  and the cutadapt software packages  <cit> .
fastqc  <cit>  was used to check the final quality of reads prior to de novo assembly.
the number of reads generated for each biological replicate is shown in table  <dig> table  <dig>
total reads and trinity results for each transcriptomic or genomic library
peruphasma schultei
aretaon asperrimus
extatosoma tiaratum
medauroidea extradentata
ramulus artemis
sipyloidea sipylus
genomic
p. schultei
reads
total number of transcripts and components based on results of the trinity assembler  <cit>  with default parameters.
n <dig> statistic is a nucleotide length.
for gut transcriptomes of the other five species, the same method was used as for p. schultei with a few changes.
for each species, only one biological replicate was produced for both the anterior and the posterior midgut.
this library was made of pooled midguts .
rna was successfully extracted for all tissues with the exception of the s. sipylus pmg, for which the extraction failed and for which no additional specimens could be obtained.
rna-extraction and quality control were as for peruphasma, but libraries were made using the nebnext® ultra™ rna library prep kit for illumina kit according to the manufacturer’s instructions.
sequencing was on an illumina hiseq  <dig>  and the data uploaded to the ncbi sra database .
the numbers of reads produced for each sample are shown in table  <dig>
de novotranscriptome assembly
the trinity assembler with the default parameters was used to generate de novo transcriptomes for all species using quality controlled reads  <cit> .
tophat  was used for aligning reads to the transcriptome  <cit> .
htseq  <cit>  was used to quantify the number of reads aligning to each transcript.
gene and isoform abundances and expression levels from the p. schultei rna-seq data were quantified using rsem   <cit> .
this program was chosen over other programs as it does not rely on reference genomes, of which there are none for the phasmatodea.
for the p. schultei rna-seq, trinity assembly yielded  <dig>  transcripts .
differentially expressed genes were identified using ebseq, an r package that compares isoform expression across two or more biological conditions, in this case amg and pmg, using a bayesian heirarchical model  <cit> .
differentially expressed genes  were those with an adjusted p-value < <dig> .
transcriptome annotation and pcwde identification
due to the lack of closely related species with well-annotated genomes, or even consensus as to what is the most closely related order to the phasmatodea, we used several methods to annotate the assembled transcripts.
for all p. schultei transcripts and the top  <dig> most highly expressed transcripts for the other species, we used blast2go’s  <cit>  tblastx program to compare each sequence to the ncbi translated nucleotide collection  database, with an expect value threshold of e- <dig>  contigs with highly significant blast  <cit>  hits were mapped to the gene ontology  database and annotated using blast2go with an expect value threshold of e- <dig>  interpro annotations were performed using the blast2go remote connection to the interproebi server  <cit> .
go terms were modulated using annex and goslim, using the “generic” mapping  available in blast2go .
potential metabolic pathways represented in the transcriptome were identified using the kyoto encyclopedia of genes and genomes   <cit>  database via blast2go .
enrichment analysis  was used to find enriched go terms, with term filter value below  <dig> , term filter mode “fdr,” and two-tailed test options selected .
annotations were added to those provided by rsem.figure  <dig>
comparisons of the go terms expressed in the anterior and posterior midguts of the phasmatodea.
the top  <dig> most expressed transcripts are seen for each of the anterior  and posterior  midguts.
a) aretaon asperrimus.
b) extatosoma tiaratum.
c) medauroidea extradentata.
d) peruphasma schultei.
e) ramulus artemis.
go categories enriched for the most differentially expressed genes in each
peruphasma schultei
midgut segment.
values are relative to the overall transcriptome as per fisher’s exact test.
a) anterior midgut .
b) posterior midgut .
to specifically identify pcwde-encoding transcripts, we downloaded nucleotide sequences for representative pcwdes from the ncbi database, selecting known, endogenous insect proteins as well as fungal, bacterial, and protozoan proteins.
the query sequences from ncbi were blast-ed against the full transcriptomes after removing low-quality reads, with an expect value threshold of e- <dig>  only transcriptome isoforms that aligned to at least 75% of the representative gene downloaded from ncbi were included in later transcript number analyses.
amino acid alignment and phylogenetic analysis
for the putative pcwdes, the transcript sequences from the phasmids were converted to amino acid sequences using the expasy online translation tool  <cit> .
a representative sequence from each isogroup  was selected based on e-value and sim mean when compared to known enzymes in the ncbi database.
the number of isogroups and isotigs  within each group is listed in table  <dig>  and additional file 7: table s <dig> shows these sequence names.
other known protein sequences for these enzymes were collected from the ncbi database from a diversity of organisms including bacteria, fungi, plants, other insects, nematodes, and, when available, protists, chordates, and other invertebrates .
the amino acid sequences were aligned using muscle  <cit>  and manually curated using mesquite  <cit> .
further alignment and production of consensus sequences for clades were done using jalview  <cit> .
these alignments were then searched for the known conserved regions for the active/catalytic sites for each enzyme type, identified using the catalytic site atlas  <cit> , with blast alignment to confirm their presence in the phasmid transcripts.table  <dig>
number of pcwde isogroups and isotigs in the phasmid midgut
data from the full midgut transcriptomes with short sequences removed.
transcripts were identified as pcwdes based on amino acid alignment to known proteins.
muscle-aligned sequences curated on mesquite were converted to phylip format  <cit> .
for neighbor-joining trees, the phylip program “seqboot” was run to make multiple datasets for bootstrapping, and the results run through “protdist” and “neighbor”, then the trees combined with “consense”.
for parsimony trees, the “seqboot” datasets were run through “protpars” and “consense”.
for maximum likelihood trees, the muscle-alignment file as uploaded to the cipres portal   <cit>  and run on raxml-hpc <dig> on xsede for  <dig> bootstrapping runs  <cit> .
for bayesian analysis, mr.bayes  <dig> .
<dig>  <cit>  was run with the cipres datasets with  <dig> generations for bootstrapping.
consensus trees were viewed and prepared for figures using figtree  <dig> .
<dig>  <cit> .
the maximum likelihood tables were chosen as the figures for this manuscript.
testing for endogenous production of pcwdes
the possibility existed that some transcripts came from microbial symbionts or contaminants.
table s <dig> shows which pcwde encoding isogroups contained poly-adenylation signals, a feature predominantly of eukaryotic mrna that that, unlike bacterial rna, will pass in large amount through our cdna synthesis method, and whose presence is used to suggest endogenicity .
we also tested the translated transcripts for the presence of eukaryote-specific signal peptides using signalp  <dig>    <cit> , which is also evidence against bacterial origins for the transcripts.
such methods however cannot differentiate between enzymes produced by protozoan or fungal symbionts and those produced by insects, including insect-produced proteins whose genes were acquired from a eukaryote via horizontal gene transfer as has happened in beetles  <cit> .
to further show that particular pcwde genes are endogenous to the phasmid genome and not produced by gut symbionts or contaminants, we extracted dna from non-gut tissue for next generation sequencing.
finding a gene encoding a transcript protein in the genome is a strong and frequently used indicator of endogenicity .
by using non-alimentary tissue we avoided aspecific amplication of microbial genes and recovered insect dna alone, as has been done in other insects to show microbe-like genes are endogenously produced , including the first discovery of endogenously produced cellulase in an insect  <cit> .
we dissected out the brain under sterile conditions from one p. schultei individual and extracted dna using the chargeswitch® gdna mini tissue kit .
genomic illumina libraries for paired-end 100 bp sequencing were then produced using the illumina truseq kit and validated using the bioanalyzer  <dig>  sequencing was conducted on the illumina hiseq  <dig>  and the data uploaded to the ncbi sra database .
the numbers of reads generated for the sample are shown in table  <dig>  we tested if the p. schultei cellulase and pectinase genes were endogenous in origin by mapping all genomic reads from the brain tissue back to our de novo midgut transcriptome assembly.
we then took all the genomic reads that mapped to each pcwde gut transcriptome gene and blasted them to the entire gut transcriptome to narrow the list of reads down to those that uniquely map to a single pcwde gene .table  <dig>
peruphasma schultei
pectinase and cellulase genomic reads uniquely aligned to transcriptome isotigs
these reads are therefore of endogenous genes.
counts to  <dig> max.
as a final test, we blasted all pcwde transcripts to the draft genome for timema cristinae , which was under development but available at the webpage of the nosil lab of evolutionary biology at the university of sheffield, uk .
the timematidae are considered the sister group to all other phasmatodea  <cit> .
while the timema genome is not guaranteed to have the same genes as the species we examined, finding our transcripts within the timema genome would be strong evidence that the gene is both endogenous in and ancestral to phasmatodea.
RESULTS
phasmid midgut de novotranscriptome assemblies
from the extracted rna libraries of the pooled amgs or pmgs of each phasmid species we generated approx  <dig> million high quality, 100 bp, paired-end sequence reads , with the exception of sipyloidea sipylus , from which we were only able to successfully extract rna from the anterior midguts.
de novo assembly of each midgut section’s library with trinity  <cit>  produced ~114- <dig> thousand transcript contigs per species .
all reads and the final transcriptome for peruphasma schultei are available under bioproject accession prjna <dig>  and for all other phasmids under prjna <dig>
annotation of the p. schulteitranscriptomes
approximately  <dig>  of the  <dig> transcriptomic sequences had blast hits , most of which were homologous to sequences from other insects and arthropods .
more genes were homologous to the red flour beetle tribolium castaneum than to insects from more closely related orders, reflecting the relative dearth of available genetic information from insects in the polyneoptera clade and the relatively recent sequencing of tribolium.
the high percentage of sequences with no blast hits  is unsurprising given the lack of an annotated, sequenced phasmatodea genome.
the sequences may have represented noncoding regions, wrongly-assembled contigs, or novel genes whose significance is unknown.
differential gene expression across the phasmid midgut
rna-seq analysis of p. schultei suggested compartmentalization of gut function , as found in other insects  <cit>  and plant cell wall consuming organisms.
additional files  <dig> and  <dig> list all differentially expressed genes  between the p. schultei midgut sections, defined as genes with a posterior probability of differential expression  > <dig>   <cit> .
we also defined genes as being highly expressed if their expression levels were 10× higher than the mean for that midgut segment.
over  <dig> genes were differentially expressed in each gut section, with  <dig> genes expressed only in the amg and  <dig> expressed only in the pmg.
we found  <dig> highly expressed genes for the amg and  <dig> for the pmg .
analysis of the most highly expressed sequences in each midgut section for all species further suggested compartmentalization of digestion, or at least enzyme gene expression.
all species showed similar go category profiles for each midgut section , with nearly 50% reduction in hydrolase gene expression in the pmg relative to the amg.
enzyme-encoding transcripts that break down polymers at internal sites, such as serine proteases, lipases, and pcwdes  <cit>  were more abundant in the amg, as were carboxylesterases transcripts and sugar hydrolases.
transcripts abundant in the pmg encoded enzymes that break down dimers and monomers, such as dipeptidases, phospholipases, and trehalase, as well as some cell membrane receptor proteins and cytochrome p450s.
phasmatodea midgut pcwdes
among the sugar hydrolases were several isogroups of pcwdes including cellulases in the gh <dig> family and cellobiase, which together can digest cellulose polymers completely into sugar, and the pectinase endopolygalacturonase.
we also found transcripts encoding beta- <dig> -glucanase , a polysaccharide-degrading enzyme family known mostly from lepidoptera larval midguts and expressed in response to feeding on a diet containing bacteria  <cit> .
these four enzymes  were used in the manual annotations.
amino acid alignment and phylogenetic analysis
the phasmid cellulases aligned most closely with other polyneopteran cellulases — including the known, active, endogenous cellulases isolated from the phasmids eurycantha calcarata and entoria okinawaensis
<cit>  — as well as those of other invertebrates, tunicates, plants, and actinomycete bacteria, but not with nematode or beetle cellulases.
phasmid endoglucanases are of the gh <dig> family thought to be ancestral to all animal life  <cit> , as opposed to the gh <dig>   <dig>  or  <dig> cellulases found in nematodes and beetles  <cit> .
the phasmid transcripts either themselves included or were homologous to transcripts including the known active sites invariant in gh <dig> cellulases, based on work on thermobifida/thermomonospora fusca   <cit> : namely two conserved asp’s  functioning in catalytic base activity and a glu residue  that functions as the catalytic acid .
phylogenetic analysis could not resolve domain or phylum-level relationships among the sequences tested .
every cellulase transcript we isolated was homologous to two sequences from the timema genome.figure  <dig>
sections of cellulase amino acid sequence alignments containing conserved/active site residues.
catalytic sites are identified by the grey arrows  <cit> .
list of phasmatodea sequences in additional file  <dig>  and all others in additional file  <dig>  similar sequences from groups of related species were combined into consensus sequences.
red letters show identity with the overall consensus sequence based on conservation of physiochemical properties.
quality is the inverse likelihood of observing mutations based on the blosum <dig> matrix  <cit> .
strongylocentrotus purpuratus and flavobacterium branchiophilum are abbreviated.
the phasmid pectinase sequences aligned most closely with those of gamma proteobacteria, rather than other insects or eukaryotes.
the alignment also showed hemipteran and beetle pectinases as most similar to fungal pectinases, the latter homology already noted in the literature  <cit> .
nearly all phasmid pectinase enzymes contained the four conserved regions of the catalytic sites, based on work on erwinia carotovora   <cit> : asn226-thr227-asp <dig>  gly248-asp249-asp <dig>  gly274-his275-gly <dig>  and arg305-ile306-lys <dig> .
an exception is the transcripts from a. asperrimus, whose arg residue is replaced with a tyrosine .
phylogenetic analysis suggested the phasmid polygalacturonases are a monophyletic group within those of the gamma proteobacteria .figure  <dig>
section of pectinase amino acid sequence alignments containing conserved/active site residues.
see caption to figure  <dig>  catalytic residues identified with grey bars  <cit> .
thermoanaerobacterium thermosaccharolyticum is abbreviated.
maximum likelihood tree for the pectinases.
phylogenetic analysis of amino acid sequences made with raxml-hpc <dig> on the xsede system  <cit> .
numbers are bootstrap values .
branch widths based on bootstrap value, branch colors based on clade.
branch lengths based on the mean number of nucleotide substitutions per site .
all phasmatodea sequences  were a monophyletic group among the gamma proteobacteria.
thermoanaerobacterium thermosaccharolyticum is abbreviated.
the possibility exists that the pectinase transcripts came from gut bacteria or contaminants rather than phasmid genes.
however, many of the pectinase transcripts had poly-a tails, as did those of other pcwdes , and all the pcwde transcripts that were not truncated at the 5’ end had eukaryotic-specific signal peptides.
this includes transcripts that contained complete open reading frames  as well as those that were 3’ truncated.
each pectinase-encoding contig also had multiple  matching genomic reads from brain tissue that uniquely aligned to them .
the same matching genomic reads could be found for cellulases, which have been demonstrated to be endogenously produced in phasmids  <cit>  and are endogenously produced in many other insects  <cit>  and metazoans  <cit> .
however, none of the pectinase transcripts had homologues in the timema genome.
the phasmid beta-glucosidases/cellobiases were in the gh <dig> family and aligned most with those of other insects.
the phasmid transcripts mostly had the conserved residues of beta-glucosidases, including the catalytic sites, based on work with white clover, trifolium repens   <cit> : arg <dig>  his <dig>  asn <dig>  glu <dig>  asn <dig>  tyr <dig>  glu <dig>  and trp <dig> .
phylogenetic analysis suggested the phasmid cellobiases are nearly all monophyletic, except a strongly-supported clade consisting of one isogroup from each species but s. sipylus .
analysis could not determine interclass relationships among insect beta-glucosidases, but suggested the enzyme existed in the common ancestor of the insecta.
every beta-glucosidase transcript we isolated was homologous to four sequences from the timema genome.figure  <dig>
sections of cellobiase amino acid sequence alignments containing conserved/active site residues.
see caption to figure  <dig>  catalytic and conserved residues identified with grey arrows  <cit> .
maximum likelihood tree for the beta-glucosidases.
phylogenetic analysis of amino acid sequences made with raxml-hpc <dig> on the xsede system  <cit> .
numbers are bootstrap values .
branch widths based on bootstrap value, branch colors based on clade.
branch lengths based on the mean number of nucleotide substitutions per site .
the phasmatodea sequences  formed two strongly supported monophyletic group with weak relationships to other insect groups.
the phasmid beta- <dig> -glucanases were in the gh <dig> family and aligned most closely with other insect enzymes , however this is a relatively recently described enzyme with few recorded sequences in the literature or ncbi database.
this is the first known record of endogenous beta- <dig> -glucanase in the polyneoptera.
the phasmid beta- <dig> -glucanases could be divided into four clear, monophyletic groups  with no more than one representative isogroup per each of the six phasmid species.
each group differed in their homology with the known consensus pattern for catalytically active beta- <dig> -glucanases based on work on bacillus licheniformis
<cit> : e--d--x-e-x---x- .
one group of six sequences had 11/ <dig> amino acids conserved, a group of four had 10/ <dig>  another group of six had 8/ <dig>  and the final group consisting of one a. asperrimus sequence had 6/ <dig>  the last amino acid in the bacillus region was not conserved among any phasmids, nor is it conserved among the lepidoptera sequences, which were also 11/ <dig>  or many other organism sequences sampled.
every beta- <dig> -glucanase transcript we isolated had six to eight homologues in the timema genome.figure  <dig>
sections of beta- <dig> -glucanase amino acid sequence alignments containing conserved/active site residues.
see caption to figure  <dig>  conserved, twelve amino acid region identified with grey bar  <cit> .
the four groups of phasmatodea gene are separated by black lines.
maximum likelihood tree for the beta- <dig> -glucanases.
phylogenetic analysis of amino acid sequences made with raxml-hpc <dig> on the xsede system  <cit> .
numbers are bootstrap values .
branch widths based on bootstrap value, branch colors based on clade.
branch lengths based on the mean number of nucleotide substitutions per site .
the phasmatodea sequences formed four strongly supported groups.
discussion
using high coverage sequencing of rna expressed in their midguts, we were able to produce high quality transcriptomes of several phasmatodea species.
this new data doubles the genera of phasmids with publicly available genetic resources on the ncbi databases, increasing the amount of annotated genes available for future work not only on phasmatodea, but also on the polyneoptera in general.
covering six species in four families, while drawing from the draft genome of a seventh species in a fifth family, the data suggests the differential expression and enzyme gene diversity of the phasmid midgut sections is mostly conserved throughout the order.
our findings will serve as a reference set for studying phasmid digestion and a jumping point for future proteomic and biochemical assays.
the abundance of pcwde isogroups in phasmids is relatively high, and the diversity of pcwde types is comparable to those in certain leaf beetles  like phaedon cochleariae
<cit>  or wood-boring beetles  like anoplophora glabripennis
<cit> .
the current record is likely diabrotica virgifera virgifera  with seventy-eight genes putatively encoding proteins from the same four enzyme classes studied here  <cit> .
as phasmatodea and chrysomelidae are among the few insect groups to be exclusively folivorous, a possible correlation exists between that dietary niche and a diverse pcwde complement.
for de novo transcriptomes, assemblers such as trinity often cannot differentiate between homologous genes and isoforms or allelic variants of the same gene.
they can potentially overestimate the number of isotigs  within an isogroup   <cit> .
combined with relatively low genetic resource availability for closely related insects and relatively high representation of species like the aforementioned beetles, we cannot be certain whether phasmids express more or fewer pcwdes than the average herbivorous insect.
in addition to using programs like rsem designed to reduce such errors  <cit> , comparing the number of reads mapping to a locus on the genome can be used to infer the true isoform number and account for inflation  <cit> .
that most p. schultei transcriptome sequences  had more matching genome reads  than their corresponding isogroup has members  suggests that our contig numbers represent true isoforms within each isogroup, rather than an overestimation due to mis-assembly  <cit> .
these phasmid isoforms may reflect multiple gene copies or alternatively spliced genes, either case suggesting a diverse complement of proteins working together to fully digest multiple varieties of carbohydrate polymer.
a highly derived genetic capacity for plant cell wall breakdown  <cit> .
however, because some isotigs were truncated at the 3’ or 5’ end, the possibility exists that certain transcripts represent different ends of a single gene.
future work using race-pcr from primers based on the transcripts identified here would produce full-length cdna sequences that will determine which transcripts represent unique genes and which are fragments, such genes could then be expressed into insect cell cultures for use in downstream enzymatic activity assays  <cit> .
previous research has confirmed that the endogenous cellulase genes we demonstrated are most highly expressed in the anterior midgut are also most highly active in the anterior midgut  <cit> , making it the site of both cellulase translation and action.
the physical structure of the amg supports this hypothesis: the pleating and folding serves to greatly increase the available surface area of the amg while slowing down the transit speed of food, increasing the amount of time and space available for cellulase enzymes to hydrolyze ingested plant material.
cellulase activity falls to nearly nothing in the pmg, tracking with cellulase gene expression.
if we extend the results of cellulases to those of the other pcwdes, then we hypothesize that phasmid digestive enzymes are active in the same region of the gut where they are expressed, making the pleated amg the site of primary plant cell wall and polymer digestion and the pmg the site of secondary digestion of smaller oligomers at most.
we also hypothesize that phasmids can fully digest cellulose into glucose, as they have the two enzymesnecessary to do so, and also actively degrade pectin into galacturonic acid.
such digestive abilities could explain how phasmids survive on otherwise uncommon, obligately folivorous diets: by fully breaking down plant cell walls into assimilatable nutrients rather than just degrading the walls to access the nutrient-rich cytoplasm within.
as transcriptomics only demonstrates gene expression, not translation or activity, these hypotheses cannot be confirmed with this data alone.
however, the fact that phasmids have active cellulases  <cit>  and the presence of the relevant catalytic residues on the phasmid cellulase, pectinase, and cellobiase transcripts  and some beta- <dig> -glucanase transcripts  suggests the transcripts code for functional enzymes, supporting the hypothesis that these enzymes are indeed actively degrading plant cell walls.
we have thus provided the necessary preliminary work justifying biochemical and proteomic assays into cellobiase, pectinase, and beta-glucanase activity in the phasmid gut.
phasmid pectinases are all endopolygalacturonases in the gh <dig> group, known in insects only from the beetles and the hemiptera, but they show homology and align to those from gamma proteobacteria .
pectinase genes were also absent in the timema genome.
our pectinase transcripts may have come from a bacterial symbiont.
the successful mapping of all p. schultei pectinase and cellulase transcripts to genes in the p. schultei genome, the presence of eukaryote-specific poly-a tails and signal peptides on phasmid transcripts with complete open reading frames, previous studies with p. schultei and r. artemis suggesting their digestion is symbiont independent  <cit> , and the absence of characteristic paunches for microbial fermentation  <cit> , all tentatively suggest these pectinases are encoded in the phasmids own genome and not produced by gut microbes.
however, the possibility remains that the samples were contaminated by a non-symbiotic microbe: either a rare bacteria that poly-adenylates its rna or a fungal symbiont that acquired a bacterial gene via horizontal transfer.
a more parsimonious hypothesis is that the phasmid pectinase gene was acquired through horizontal transfer from a bacterial ancestor, much as the beetle pectinases are thought to have been acquired though horizontal transfer from an ascomycete fungus  <cit> , or as leaf beetle xylanases may have also transferred from a gamma proteobacteria  <cit> .
the absence of the genes in timema would suggest either that the transfer event occurred after the split between the timematidae and the other phasmatodea, or that the pectinase genes are ancestral to both and lost in timema.
lastly, the similarity between phasmid and bacterial sequences could simply be an artefact of the over-representation of microbial and dearth of animal pectinases in the ncbi database at this time  <cit> .
using long range or race pcr to clone entire genes from phasmid genomic dna and get introns would conclusively demonstrate whether or not the pectinases are endogenously produced or not, and such work on all six species studied here is underway.
whether pectinase genes exist in other polyneoptera remains to be seen, but would help determine when the horizontal transfer event could have occurred.
failure to find endogenous pectinases in closely related insects would mean the transfer occurred in an early phasmatodea ancestor: a development that would have expanded the digestive abilities of the order and may have played a significant role in their evolution of obligate folivory.
pcwde diversity could also be correlated to the development of the longer and larger body sizes of the euphasmatodea.
broader, multi-phyla, phylogenetic analysis for pectinolytic enzyme genes in the animalia can answer this question  <cit> .
phasmid cellobiases are gh <dig>  as are those of other insects who produce them endogenously like the higher termites  <cit> .
such species can break down cellobiose independently, unlike the lower termites that have symbiotic microorganisms to produce their cellobiases for them.
the beta- <dig> -glucanases are gh <dig>  similar to those found in lepidoptera  <cit> , however the prevalence of this recently described gene in animals has not been sufficiently examined.
a greater sampling of this gene's presence in other animal, fungal, and bacterial species, as well as biochemical studies to determine the conserved catalytic residues for the protein, are needed before ordinal-level hypotheses can be made for the enzyme's evolutionary history.
so far, the animal enzymes appear most closely related, and we hypothesize that at least three distinct beta- <dig> -glucanase gene families existed in an early phasmatodea ancestor, including the ancestor of the timematidae.
our work promotes phasmids as a potentially high-value source of novel pcwdes for biotechnological applications  <cit> .
cellulases and pectinases are highly sought after by the biofuel industry to degrade feedstock into the monomers later converted into fuel, or to improve the flow rate of the material by reducing the amount of solid matter  <cit> .
pectinases are also used in the production of coffee, tea, and juice, and in waste-water treatment  <cit> .
phasmid pcwdes could be introduced into bacteria or fungi like trichoderma reesei, for industrial-scale enzyme production or direct use in bioreactors for wastewater treatment or biofuel production  <cit> .
our de novo midgut transcriptomes enabled us to survey all expressed pcwdes of the phasmatodea at once and identify conserved catalytic domains, justifying downstream translation and activity level analyses.
a benefit of this is system is the increased speed and efficiency compared to the converse  <cit> : running chemical assays to identify enzyme activity, using proteomics to identify the amino acid sequence of isolated enzymes, and working backwards from there to design a primer for the enzyme-encoding gene and hope it exists within the target organism’s genome itself and not a symbiont or contaminant  <cit> .
another benefit is that transcriptomics can reveal genes useful for phylogenetic analysis but that are not translated or whose proteins are modified post-translation such that the standard biochemical tests do not detect their function.
an associated drawback is that expressed genes are not necessarily translated into active proteins, nor are they necessarily active at the site of expression  <cit> .
however, when a reference genome is not available, a transcriptome can provide large sets of potential genes for study and, combined with genomic data, can determine whether or not they are endogenous to the target organism, which cannot be determined by homology alone.
de novo transcriptome assembly combined with rna-seq is a powerful tool for suggesting putative functions for unknown tissues in understudied organisms and directions for future study.
CONCLUSIONS
the folivorous phasmatodea are an ideal system to study the evolution of obligate herbivory, yet a paucity of genetic resources and poorly understood basic biology impede such work.
using rna-seq, we demonstrated a diversity of plant cell wall degrading enzymes expressed differentially in the anterior section of the phasmid midgut.
of these, the cellulases, cellobiases, and beta- <dig> -glucanases are likely all encoded in the insect’s own genome, as are the pectinases, though we could not definitively rule out a microbial source for the latter.
such an abundance of endogenous enzymes was not expected from the polyneoptera, raising important questions on their evolutionary history.
the efficiency by which our de novo transcriptomes generated new genomic resources and hypotheses for future research on polyneopteran digestion demonstrate the power of such methods to analyze organisms lacking sequenced genomes.
our findings strongly encourage expanding the searches for pcwdes, most notably the pectinases and beta- <dig> -glucanases, into other, lower polyneopteran insects.
availability of supporting data
all reads and sequence files described in the manuscript are available under bioproject accessions prjna <dig> for p. schultei and prjna <dig> for the other phasmids.
electronic supplementary material
additional file 1: table s1: kegg table for the most highly expressed genes of the aretaon asperrimus midgut.
top  <dig> most highly expressed genes in each of the anterior midgut  and posterior midgut  used.
additional file 2: table s2: kegg table for the most highly expressed genes of the extatosoma tiaratum midgut.
top  <dig> most highly expressed genes in each of the anterior midgut  and posterior midgut  used.
additional file 3: table s3: kegg table for the most highly expressed genes of the medauroidea extradentata midgut.
top  <dig> most highly expressed genes in each of the anterior midgut  and posterior midgut  used.
additional file 4: table s4: kegg table for the most highly expressed genes of the peruphasma schultei midgut.
top  <dig> most highly expressed genes in each of the anterior midgut  and posterior midgut  used.
additional file 5: table s5: kegg table for the most highly expressed genes of the ramulus artemis midgut.
top  <dig> most highly expressed genes in each of the anterior midgut  and posterior midgut  used.
additional file 6: table s6: kegg table for the most highly expressed genes of the sipyloidea sipylus midgut.
top  <dig> most highly expressed genes in each of the anterior midgut  only.
additional file 7: table s7: representative isotigs  for each pcwde isogroup for the phasmatodea.
data is from the full midgut transcriptomes with short sequences removed.
transcripts were identified as pcwdes based on amino acid alignment to known proteins from the ncbi database.
* = the sequence contained a complete open-reading frame.
other sequences could have been truncated at the 5’ or 3’ or both.
† = the sequence or another in that isogroup  had a poly-a tail.
additional file 8: table s8: – species and ncbi accession no's for cellulase  proteins compared to phasmid proteins.
sequences with *were not included in the alignment figure  due to poor or absent overlap with other sequences at that region.
table s <dig>  – species and ncbi accession no's for pectinase  proteins compared to phasmid proteins.
sequences with *were not included in the alignment figure  due to poor or absent overlap with other sequences at that region.
table s <dig>  – species and ncbi accession no's for cellobiase  proteins compared to phasmid proteins.
sequences with *were not included in the alignment figure  due to poor or absent overlap with other sequences at that region.
table s <dig>  – species and ncbi accession no's for beta- <dig> -glucanase proteins compared to phasmid proteins.
sequences with *were not included in the alignment figure figure 9) due to poor or absent overlap with other sequences at that region.
additional file 9: table s9: number of peruphasma schultei midgut transcriptome sequences with successful blast search, mapping, and annotation.
additional file 10: figure s1: species distribution for top-hit blast results of p. schultei midgut transcriptome.
additional file 11: table s12: the most differentially expressed genes  in the p. schultei amg.
includes genes found in both or only one tissue type.
means measured in rpkm.
ppde = posterior probability of differential expression.
annotations made with blast2go.
additional file 12: table s13: the most differentially expressed genes  in the p. schultei pmg.
includes genes found in both or only one tissue type.
means measured in rpkm.
ppde = posterior probability of differential expression.
annotations made with blast2go.
additional file 13: table s14: the most highly expressed genes of the p. schultei amg.
identified as genes with expression levels ten times greater than the mean for that section.
means measured in rpkm.
annotations made with blast2go.
additional file 14: table s15: the most highly expressed genes of the p. schultei pmg.
identified as genes with expression levels ten times greater than the mean for that section.
means measured in rpkm.
annotations made with blast2go.
abbreviations
amganterior midgut
ciprescyberinfrastructure for phylogenetic research
degdifferentially expressed genes
ecenzyme commission
ghglycoside hydrolase
gogene ontology
heghighly expressed gene
keggkyoto encyclopedia of genes and genomes
ncbinational center for biotechnology information
pcwdeplant cell wall degrading enzyme
pdbprotein data bank
pmgposterior midgut
ppdeposterior probability of differential expression
rpkmreads per kilobase per million
rsemrna-seq by expectation maximization
srasequence read archive
xsedeextreme science and engineering discovery environment.
competing interests
the authors declare that they have no competing interests.
authors’ contributions
authors ms and brj conceived the study and contributed equally to the work.
lsk and brj provided resources and expertise.
ms and lsk reared the insects.
wcj, ja, and brj carried out the rna seq and de novo transcriptome assembly.
ms, wcj, ja, and brj analyzed and interpreted data.
all authors read and approved the final manuscript.

 whole transcriptome shotgun sequencing, or rna-seq, is a high-throughput, next-generation sequencing tool that can efficiently identify tens of thousands of functional genes in an organism or specific tissue at a given time. this deep sequencing makes it a more attractive tool than microarrays for organisms lacking reference genomes, facilitating de novo transcriptome assembly. its high coverage is desirable when profiling transcripts in tissues of unknown function, enabling researchers to generate and/or test multiple hypotheses at once, and in organisms potentially harboring symbiotic organisms that may or may not produce the transcripts of interest, as rna-seq can simultaneously identify genes from microbes and their vectors/hosts. such a combination of low genome resource availability and enigmatic physiology exists in the stick and leaf insects, or phasmids. though common in the pre-molecular biology era through the laboratory stick insect, carausius morosus, phasmid research today is relatively limited. few phasmids are pests of agricultural crops, though they reach plague-like abundances in temperate forests and c. morosus is an invasive pest in several countries. all life stages of all species within the order feed exclusively on leaves. this obligate folivory is relatively rare: among insects it is known only from leaf beetles, while more basal “herbivores” such as grasshoppers and crickets will quite readily scavenge vertebrate meat, engage in cannibalism, and even hunt and kill other insects. thus phasmids are an ideal system for studying the evolution of herbivory in the lower polyneoptera. folivorous organisms benefit greatly from plant cell wall degrading enzymes, a group that includes cellulases, hemicellulases, lignases, pectinases, and xylanases. once thought to be limited to microbes, endogenous pcwde production has since been found throughout the animalia. in particular, cellulase genes from the glycoside hydrolase family are now believed to have existed in the ancestor of all metazoan life as opposed to having been repeatedly acquired from microbes via horizontal gene transfer, as is thought to be origin of gh and gh cellulases in beetles. among insects, endogenous cellulases have been found in lower and higher termites, cockroaches, crickets, beetles, a firebrat, a springtail, and, recently, the phasmids. high cellulase activity in the anterior midguts of two phasmid species, eurycantha calcarata and entoria okinawaensis, was detected, the responsible proteins isolated, and the genes encoding them sequenced. sequence homology and antigency against an insect cellulase anti-serum supported an endogenous, insectan origin for the enzymes. such a process is slow and predicated on the translation of pcwde genes into enzymes active against laboratory substrates like carboxymethylcellulose or crystalline cellulose, whose specificity and selectivity are imperfect. whether or not phasmids contain other pcwdes, such as the cellobiases that convert the products of cellulase into glucose monomers or the pectinases that hydrolyze pectin into galacturonic acid monomers, was unknown. presence of such active enzymes could explain the obligate folivory of the phasmatodea and be a key factor in the order’s evolution, which is itself a puzzle as the sister order to the phasmatodea is highly debated. microbiological assays of the phasmid gut suggest digestion in the order is symbiont-independent, so any phasmid pcwde genes are likely endogenous, yet whether the genes show homology to insect or microbe genes depends on their own evolutionary origin. complicating the issue is the relative lack of genetic resources for phasmids or their most closely-related orders: orthoptera, embioptera, and notoptera/xenonomia. lastly, even if phasmids have pcwde genes, their expression is not a guarantee, nor are they necessarily expressed in the gut region where they are most active. the phasmid midgut is physically differentiated into two sections: an anterior midgut marked by circular pleating and folding, and a posterior midgut studded irregularly by hollow bulbs with filamentous tubules called the “appendices of the midgut” that open into the midgut lumen. the appendices may have an excretory or secretory function, and the possibility exists that they produce digestive enzymes that are carried forward into the amg via countercurrent flow. in the face of all these unknowns, next-generation sequencing is the best resource for answering questions of phasmid digestive physiology efficiently and effectively. here we used de novo transcriptome assembly to identify the genes expressed in the midguts of six species of phasmatodea from four families, while greatly increasing the publicly available genetic resources for the order. we also used rna-seq on one exemplar species, peruphasma schultei to quantitatively compare transcript expression between the amg and pmg, and produced a genomic dna library from the symbiont-free phasmid brain to confirm that identified transcripts were encoded by the insect itself. our main goal was to identify the production organ of the phasmatodea endogenous cellulase, while simultaneously creating an inventory of expressed pcwde and other digestive genes in phasmids and generating hypotheses on their evolutionary origins and the putative functions of the midgut sections. this study serves as a necessary preliminary for more targeted molecular work. more broadly, our transcriptomes are useful for evolutionary analyses of non-cellulase pcwdes in insects and identifying potential genes with biotechnological applications such as in processing biofuel feedstock or improving its rheology. methods insects and microscopy insects used were peruphasma schultei, sipyloidea sipylus, aretaon asperrimus, and extatosoma tiaratum, medauroidea extradentata, and ramulus artemis cultured at room temperature in the bohart museum of entomology, university of california, davis. phasmids were fed an ad libitum diet of privet for peruphasma, eucalyptus for extatosoma, and rosa sp. for the others. library prep and sequencing the rna-seq study of peruphasma schultei made use of three biological replicates for both the anterior and the posterior midguts. for each replicate, the guts of five fed, surface-sterilized, adult, male and female phasmids were removed under sterile conditions and emptied of their contents in several washes of 70% ethanol. then the anterior and posterior sections were separately pooled and homogenized in trizol® reagent. rna was extracted according the trizol-plus protocol, which includes an on-column dnaase digestion step. total rna quality was checked with the bioanalyzer libraries were made using the illumina truseq v kit according to the manufacturer’s instructions. hundred base pair paired-end sequencing was performed on the hiseq and the raw data uploaded to the ncbi sra database. for quality control, low quality bases and adapter contamination were removed with the fastx toolkit and the cutadapt software packages. fastqc was used to check the final quality of reads prior to de novo assembly. the number of reads generated for each biological replicate is shown in table  table total reads and trinity results for each transcriptomic or genomic library peruphasma schultei aretaon asperrimus extatosoma tiaratum medauroidea extradentata ramulus artemis sipyloidea sipylus genomic p. schultei reads total number of transcripts and components based on results of the trinity assembler with default parameters. n statistic is a nucleotide length. for gut transcriptomes of the other five species, the same method was used as for p. schultei with a few changes. for each species, only one biological replicate was produced for both the anterior and the posterior midgut. this library was made of pooled midguts. rna was successfully extracted for all tissues with the exception of the s. sipylus pmg, for which the extraction failed and for which no additional specimens could be obtained. rna-extraction and quality control were as for peruphasma, but libraries were made using the nebnext® ultra™ rna library prep kit for illumina kit according to the manufacturer’s instructions. sequencing was on an illumina hiseq and the data uploaded to the ncbi sra database. the numbers of reads produced for each sample are shown in table  de novotranscriptome assembly the trinity assembler with the default parameters was used to generate de novo transcriptomes for all species using quality controlled reads. tophat was used for aligning reads to the transcriptome. htseq was used to quantify the number of reads aligning to each transcript. gene and isoform abundances and expression levels from the p. schultei rna-seq data were quantified using rsem. this program was chosen over other programs as it does not rely on reference genomes, of which there are none for the phasmatodea. for the p. schultei rna-seq, trinity assembly yielded transcripts. differentially expressed genes were identified using ebseq, an r package that compares isoform expression across two or more biological conditions, in this case amg and pmg, using a bayesian heirarchical model. differentially expressed genes were those with an adjusted p-value <. transcriptome annotation and pcwde identification due to the lack of closely related species with well-annotated genomes, or even consensus as to what is the most closely related order to the phasmatodea, we used several methods to annotate the assembled transcripts. for all p. schultei transcripts and the top most highly expressed transcripts for the other species, we used blast2go’s tblastx program to compare each sequence to the ncbi translated nucleotide collection database, with an expect value threshold of e- contigs with highly significant blast hits were mapped to the gene ontology database and annotated using blast2go with an expect value threshold of e- interpro annotations were performed using the blast2go remote connection to the interproebi server. go terms were modulated using annex and goslim, using the “generic” mapping available in blast2go. potential metabolic pathways represented in the transcriptome were identified using the kyoto encyclopedia of genes and genomes database via blast2go. enrichment analysis was used to find enriched go terms, with term filter value below, term filter mode “fdr,” and two-tailed test options selected. annotations were added to those provided by rsem.figure comparisons of the go terms expressed in the anterior and posterior midguts of the phasmatodea. the top most expressed transcripts are seen for each of the anterior and posterior midguts. a) aretaon asperrimus. b) extatosoma tiaratum. c) medauroidea extradentata. d) peruphasma schultei. e) ramulus artemis. go categories enriched for the most differentially expressed genes in each peruphasma schultei midgut segment. values are relative to the overall transcriptome as per fisher’s exact test. a) anterior midgut. b) posterior midgut. to specifically identify pcwde-encoding transcripts, we downloaded nucleotide sequences for representative pcwdes from the ncbi database, selecting known, endogenous insect proteins as well as fungal, bacterial, and protozoan proteins. the query sequences from ncbi were blast-ed against the full transcriptomes after removing low-quality reads, with an expect value threshold of e- only transcriptome isoforms that aligned to at least 75% of the representative gene downloaded from ncbi were included in later transcript number analyses. amino acid alignment and phylogenetic analysis for the putative pcwdes, the transcript sequences from the phasmids were converted to amino acid sequences using the expasy online translation tool. a representative sequence from each isogroup was selected based on e-value and sim mean when compared to known enzymes in the ncbi database. the number of isogroups and isotigs within each group is listed in table  and additional file 7: table s shows these sequence names. other known protein sequences for these enzymes were collected from the ncbi database from a diversity of organisms including bacteria, fungi, plants, other insects, nematodes, and, when available, protists, chordates, and other invertebrates. the amino acid sequences were aligned using muscle and manually curated using mesquite. further alignment and production of consensus sequences for clades were done using jalview. these alignments were then searched for the known conserved regions for the active/catalytic sites for each enzyme type, identified using the catalytic site atlas, with blast alignment to confirm their presence in the phasmid transcripts.table number of pcwde isogroups and isotigs in the phasmid midgut data from the full midgut transcriptomes with short sequences removed. transcripts were identified as pcwdes based on amino acid alignment to known proteins. muscle-aligned sequences curated on mesquite were converted to phylip format. for neighbor-joining trees, the phylip program “seqboot” was run to make multiple datasets for bootstrapping, and the results run through “protdist” and “neighbor”, then the trees combined with “consense”. for parsimony trees, the “seqboot” datasets were run through “protpars” and “consense”. for maximum likelihood trees, the muscle-alignment file as uploaded to the cipres portal and run on raxml-hpc on xsede for bootstrapping runs. for bayesian analysis, mr.bayes. was run with the cipres datasets with generations for bootstrapping. consensus trees were viewed and prepared for figures using figtree.. the maximum likelihood tables were chosen as the figures for this manuscript. testing for endogenous production of pcwdes the possibility existed that some transcripts came from microbial symbionts or contaminants. table s shows which pcwde encoding isogroups contained poly-adenylation signals, a feature predominantly of eukaryotic mrna that that, unlike bacterial rna, will pass in large amount through our cdna synthesis method, and whose presence is used to suggest endogenicity. we also tested the translated transcripts for the presence of eukaryote-specific signal peptides using signalp, which is also evidence against bacterial origins for the transcripts. such methods however cannot differentiate between enzymes produced by protozoan or fungal symbionts and those produced by insects, including insect-produced proteins whose genes were acquired from a eukaryote via horizontal gene transfer as has happened in beetles. to further show that particular pcwde genes are endogenous to the phasmid genome and not produced by gut symbionts or contaminants, we extracted dna from non-gut tissue for next generation sequencing. finding a gene encoding a transcript protein in the genome is a strong and frequently used indicator of endogenicity. by using non-alimentary tissue we avoided aspecific amplication of microbial genes and recovered insect dna alone, as has been done in other insects to show microbe-like genes are endogenously produced, including the first discovery of endogenously produced cellulase in an insect. we dissected out the brain under sterile conditions from one p. schultei individual and extracted dna using the chargeswitch® gdna mini tissue kit. genomic illumina libraries for paired-end 100 bp sequencing were then produced using the illumina truseq kit and validated using the bioanalyzer sequencing was conducted on the illumina hiseq and the data uploaded to the ncbi sra database. the numbers of reads generated for the sample are shown in table  we tested if the p. schultei cellulase and pectinase genes were endogenous in origin by mapping all genomic reads from the brain tissue back to our de novo midgut transcriptome assembly. we then took all the genomic reads that mapped to each pcwde gut transcriptome gene and blasted them to the entire gut transcriptome to narrow the list of reads down to those that uniquely map to a single pcwde gene.table peruphasma schultei pectinase and cellulase genomic reads uniquely aligned to transcriptome isotigs these reads are therefore of endogenous genes. counts to max. as a final test, we blasted all pcwde transcripts to the draft genome for timema cristinae, which was under development but available at the webpage of the nosil lab of evolutionary biology at the university of sheffield, uk. the timematidae are considered the sister group to all other phasmatodea. while the timema genome is not guaranteed to have the same genes as the species we examined, finding our transcripts within the timema genome would be strong evidence that the gene is both endogenous in and ancestral to phasmatodea.  phasmid midgut de novotranscriptome assemblies from the extracted rna libraries of the pooled amgs or pmgs of each phasmid species we generated approx million high quality, 100 bp, paired-end sequence reads, with the exception of sipyloidea sipylus, from which we were only able to successfully extract rna from the anterior midguts. de novo assembly of each midgut section’s library with trinity produced ~114- thousand transcript contigs per species. all reads and the final transcriptome for peruphasma schultei are available under bioproject accession prjna and for all other phasmids under prjna annotation of the p. schulteitranscriptomes approximately of the transcriptomic sequences had blast hits, most of which were homologous to sequences from other insects and arthropods. more genes were homologous to the red flour beetle tribolium castaneum than to insects from more closely related orders, reflecting the relative dearth of available genetic information from insects in the polyneoptera clade and the relatively recent sequencing of tribolium. the high percentage of sequences with no blast hits is unsurprising given the lack of an annotated, sequenced phasmatodea genome. the sequences may have represented noncoding regions, wrongly-assembled contigs, or novel genes whose significance is unknown. differential gene expression across the phasmid midgut rna-seq analysis of p. schultei suggested compartmentalization of gut function, as found in other insects and plant cell wall consuming organisms. additional files and list all differentially expressed genes between the p. schultei midgut sections, defined as genes with a posterior probability of differential expression >. we also defined genes as being highly expressed if their expression levels were 10× higher than the mean for that midgut segment. over genes were differentially expressed in each gut section, with genes expressed only in the amg and expressed only in the pmg. we found highly expressed genes for the amg and for the pmg. analysis of the most highly expressed sequences in each midgut section for all species further suggested compartmentalization of digestion, or at least enzyme gene expression. all species showed similar go category profiles for each midgut section, with nearly 50% reduction in hydrolase gene expression in the pmg relative to the amg. enzyme-encoding transcripts that break down polymers at internal sites, such as serine proteases, lipases, and pcwdes were more abundant in the amg, as were carboxylesterases transcripts and sugar hydrolases. transcripts abundant in the pmg encoded enzymes that break down dimers and monomers, such as dipeptidases, phospholipases, and trehalase, as well as some cell membrane receptor proteins and cytochrome p450s. phasmatodea midgut pcwdes among the sugar hydrolases were several isogroups of pcwdes including cellulases in the gh family and cellobiase, which together can digest cellulose polymers completely into sugar, and the pectinase endopolygalacturonase. we also found transcripts encoding beta- -glucanase, a polysaccharide-degrading enzyme family known mostly from lepidoptera larval midguts and expressed in response to feeding on a diet containing bacteria. these four enzymes were used in the manual annotations. amino acid alignment and phylogenetic analysis the phasmid cellulases aligned most closely with other polyneopteran cellulases — including the known, active, endogenous cellulases isolated from the phasmids eurycantha calcarata and entoria okinawaensis — as well as those of other invertebrates, tunicates, plants, and actinomycete bacteria, but not with nematode or beetle cellulases. phasmid endoglucanases are of the gh family thought to be ancestral to all animal life, as opposed to the gh or cellulases found in nematodes and beetles. the phasmid transcripts either themselves included or were homologous to transcripts including the known active sites invariant in gh cellulases, based on work on thermobifida/thermomonospora fusca: namely two conserved asp’s functioning in catalytic base activity and a glu residue that functions as the catalytic acid. phylogenetic analysis could not resolve domain or phylum-level relationships among the sequences tested. every cellulase transcript we isolated was homologous to two sequences from the timema genome.figure sections of cellulase amino acid sequence alignments containing conserved/active site residues. catalytic sites are identified by the grey arrows. list of phasmatodea sequences in additional file and all others in additional file similar sequences from groups of related species were combined into consensus sequences. red letters show identity with the overall consensus sequence based on conservation of physiochemical properties. quality is the inverse likelihood of observing mutations based on the blosum matrix. strongylocentrotus purpuratus and flavobacterium branchiophilum are abbreviated. the phasmid pectinase sequences aligned most closely with those of gamma proteobacteria, rather than other insects or eukaryotes. the alignment also showed hemipteran and beetle pectinases as most similar to fungal pectinases, the latter homology already noted in the literature. nearly all phasmid pectinase enzymes contained the four conserved regions of the catalytic sites, based on work on erwinia carotovora: asn226-thr227-asp gly248-asp249-asp gly274-his275-gly and arg305-ile306-lys. an exception is the transcripts from a. asperrimus, whose arg residue is replaced with a tyrosine. phylogenetic analysis suggested the phasmid polygalacturonases are a monophyletic group within those of the gamma proteobacteria.figure section of pectinase amino acid sequence alignments containing conserved/active site residues. see caption to figure  catalytic residues identified with grey bars. thermoanaerobacterium thermosaccharolyticum is abbreviated. maximum likelihood tree for the pectinases. phylogenetic analysis of amino acid sequences made with raxml-hpc on the xsede system. numbers are bootstrap values. branch widths based on bootstrap value, branch colors based on clade. branch lengths based on the mean number of nucleotide substitutions per site. all phasmatodea sequences were a monophyletic group among the gamma proteobacteria. thermoanaerobacterium thermosaccharolyticum is abbreviated. the possibility exists that the pectinase transcripts came from gut bacteria or contaminants rather than phasmid genes. however, many of the pectinase transcripts had poly-a tails, as did those of other pcwdes, and all the pcwde transcripts that were not truncated at the 5’ end had eukaryotic-specific signal peptides. this includes transcripts that contained complete open reading frames as well as those that were 3’ truncated. each pectinase-encoding contig also had multiple matching genomic reads from brain tissue that uniquely aligned to them. the same matching genomic reads could be found for cellulases, which have been demonstrated to be endogenously produced in phasmids and are endogenously produced in many other insects and metazoans. however, none of the pectinase transcripts had homologues in the timema genome. the phasmid beta-glucosidases/cellobiases were in the gh family and aligned most with those of other insects. the phasmid transcripts mostly had the conserved residues of beta-glucosidases, including the catalytic sites, based on work with white clover, trifolium repens: arg his asn glu asn tyr glu and trp. phylogenetic analysis suggested the phasmid cellobiases are nearly all monophyletic, except a strongly-supported clade consisting of one isogroup from each species but s. sipylus. analysis could not determine interclass relationships among insect beta-glucosidases, but suggested the enzyme existed in the common ancestor of the insecta. every beta-glucosidase transcript we isolated was homologous to four sequences from the timema genome.figure sections of cellobiase amino acid sequence alignments containing conserved/active site residues. see caption to figure  catalytic and conserved residues identified with grey arrows. maximum likelihood tree for the beta-glucosidases. phylogenetic analysis of amino acid sequences made with raxml-hpc on the xsede system. numbers are bootstrap values. branch widths based on bootstrap value, branch colors based on clade. branch lengths based on the mean number of nucleotide substitutions per site. the phasmatodea sequences formed two strongly supported monophyletic group with weak relationships to other insect groups. the phasmid beta- -glucanases were in the gh family and aligned most closely with other insect enzymes, however this is a relatively recently described enzyme with few recorded sequences in the literature or ncbi database. this is the first known record of endogenous beta- -glucanase in the polyneoptera. the phasmid beta- -glucanases could be divided into four clear, monophyletic groups with no more than one representative isogroup per each of the six phasmid species. each group differed in their homology with the known consensus pattern for catalytically active beta- -glucanases based on work on bacillus licheniformis: e--d--x-e-x---x-. one group of six sequences had 11/ amino acids conserved, a group of four had 10/ another group of six had 8/ and the final group consisting of one a. asperrimus sequence had 6/ the last amino acid in the bacillus region was not conserved among any phasmids, nor is it conserved among the lepidoptera sequences, which were also 11/ or many other organism sequences sampled. every beta- -glucanase transcript we isolated had six to eight homologues in the timema genome.figure sections of beta- -glucanase amino acid sequence alignments containing conserved/active site residues. see caption to figure  conserved, twelve amino acid region identified with grey bar. the four groups of phasmatodea gene are separated by black lines. maximum likelihood tree for the beta- -glucanases. phylogenetic analysis of amino acid sequences made with raxml-hpc on the xsede system. numbers are bootstrap values. branch widths based on bootstrap value, branch colors based on clade. branch lengths based on the mean number of nucleotide substitutions per site. the phasmatodea sequences formed four strongly supported groups. discussion using high coverage sequencing of rna expressed in their midguts, we were able to produce high quality transcriptomes of several phasmatodea species. this new data doubles the genera of phasmids with publicly available genetic resources on the ncbi databases, increasing the amount of annotated genes available for future work not only on phasmatodea, but also on the polyneoptera in general. covering six species in four families, while drawing from the draft genome of a seventh species in a fifth family, the data suggests the differential expression and enzyme gene diversity of the phasmid midgut sections is mostly conserved throughout the order. our findings will serve as a reference set for studying phasmid digestion and a jumping point for future proteomic and biochemical assays. the abundance of pcwde isogroups in phasmids is relatively high, and the diversity of pcwde types is comparable to those in certain leaf beetles like phaedon cochleariae or wood-boring beetles like anoplophora glabripennis. the current record is likely diabrotica virgifera virgifera with seventy-eight genes putatively encoding proteins from the same four enzyme classes studied here. as phasmatodea and chrysomelidae are among the few insect groups to be exclusively folivorous, a possible correlation exists between that dietary niche and a diverse pcwde complement. for de novo transcriptomes, assemblers such as trinity often cannot differentiate between homologous genes and isoforms or allelic variants of the same gene. they can potentially overestimate the number of isotigs within an isogroup. combined with relatively low genetic resource availability for closely related insects and relatively high representation of species like the aforementioned beetles, we cannot be certain whether phasmids express more or fewer pcwdes than the average herbivorous insect. in addition to using programs like rsem designed to reduce such errors, comparing the number of reads mapping to a locus on the genome can be used to infer the true isoform number and account for inflation. that most p. schultei transcriptome sequences had more matching genome reads than their corresponding isogroup has members suggests that our contig numbers represent true isoforms within each isogroup, rather than an overestimation due to mis-assembly. these phasmid isoforms may reflect multiple gene copies or alternatively spliced genes, either case suggesting a diverse complement of proteins working together to fully digest multiple varieties of carbohydrate polymer. a highly derived genetic capacity for plant cell wall breakdown. however, because some isotigs were truncated at the 3’ or 5’ end, the possibility exists that certain transcripts represent different ends of a single gene. future work using race-pcr from primers based on the transcripts identified here would produce full-length cdna sequences that will determine which transcripts represent unique genes and which are fragments, such genes could then be expressed into insect cell cultures for use in downstream enzymatic activity assays. previous research has confirmed that the endogenous cellulase genes we demonstrated are most highly expressed in the anterior midgut are also most highly active in the anterior midgut, making it the site of both cellulase translation and action. the physical structure of the amg supports this hypothesis: the pleating and folding serves to greatly increase the available surface area of the amg while slowing down the transit speed of food, increasing the amount of time and space available for cellulase enzymes to hydrolyze ingested plant material. cellulase activity falls to nearly nothing in the pmg, tracking with cellulase gene expression. if we extend the results of cellulases to those of the other pcwdes, then we hypothesize that phasmid digestive enzymes are active in the same region of the gut where they are expressed, making the pleated amg the site of primary plant cell wall and polymer digestion and the pmg the site of secondary digestion of smaller oligomers at most. we also hypothesize that phasmids can fully digest cellulose into glucose, as they have the two enzymesnecessary to do so, and also actively degrade pectin into galacturonic acid. such digestive abilities could explain how phasmids survive on otherwise uncommon, obligately folivorous diets: by fully breaking down plant cell walls into assimilatable nutrients rather than just degrading the walls to access the nutrient-rich cytoplasm within. as transcriptomics only demonstrates gene expression, not translation or activity, these hypotheses cannot be confirmed with this data alone. however, the fact that phasmids have active cellulases and the presence of the relevant catalytic residues on the phasmid cellulase, pectinase, and cellobiase transcripts and some beta- -glucanase transcripts suggests the transcripts code for functional enzymes, supporting the hypothesis that these enzymes are indeed actively degrading plant cell walls. we have thus provided the necessary preliminary work justifying biochemical and proteomic assays into cellobiase, pectinase, and beta-glucanase activity in the phasmid gut. phasmid pectinases are all endopolygalacturonases in the gh group, known in insects only from the beetles and the hemiptera, but they show homology and align to those from gamma proteobacteria. pectinase genes were also absent in the timema genome. our pectinase transcripts may have come from a bacterial symbiont. the successful mapping of all p. schultei pectinase and cellulase transcripts to genes in the p. schultei genome, the presence of eukaryote-specific poly-a tails and signal peptides on phasmid transcripts with complete open reading frames, previous studies with p. schultei and r. artemis suggesting their digestion is symbiont independent, and the absence of characteristic paunches for microbial fermentation, all tentatively suggest these pectinases are encoded in the phasmids own genome and not produced by gut microbes. however, the possibility remains that the samples were contaminated by a non-symbiotic microbe: either a rare bacteria that poly-adenylates its rna or a fungal symbiont that acquired a bacterial gene via horizontal transfer. a more parsimonious hypothesis is that the phasmid pectinase gene was acquired through horizontal transfer from a bacterial ancestor, much as the beetle pectinases are thought to have been acquired though horizontal transfer from an ascomycete fungus, or as leaf beetle xylanases may have also transferred from a gamma proteobacteria. the absence of the genes in timema would suggest either that the transfer event occurred after the split between the timematidae and the other phasmatodea, or that the pectinase genes are ancestral to both and lost in timema. lastly, the similarity between phasmid and bacterial sequences could simply be an artefact of the over-representation of microbial and dearth of animal pectinases in the ncbi database at this time. using long range or race pcr to clone entire genes from phasmid genomic dna and get introns would conclusively demonstrate whether or not the pectinases are endogenously produced or not, and such work on all six species studied here is underway. whether pectinase genes exist in other polyneoptera remains to be seen, but would help determine when the horizontal transfer event could have occurred. failure to find endogenous pectinases in closely related insects would mean the transfer occurred in an early phasmatodea ancestor: a development that would have expanded the digestive abilities of the order and may have played a significant role in their evolution of obligate folivory. pcwde diversity could also be correlated to the development of the longer and larger body sizes of the euphasmatodea. broader, multi-phyla, phylogenetic analysis for pectinolytic enzyme genes in the animalia can answer this question. phasmid cellobiases are gh as are those of other insects who produce them endogenously like the higher termites. such species can break down cellobiose independently, unlike the lower termites that have symbiotic microorganisms to produce their cellobiases for them. the beta- -glucanases are gh similar to those found in lepidoptera, however the prevalence of this recently described gene in animals has not been sufficiently examined. a greater sampling of this gene's presence in other animal, fungal, and bacterial species, as well as biochemical studies to determine the conserved catalytic residues for the protein, are needed before ordinal-level hypotheses can be made for the enzyme's evolutionary history. so far, the animal enzymes appear most closely related, and we hypothesize that at least three distinct beta- -glucanase gene families existed in an early phasmatodea ancestor, including the ancestor of the timematidae. our work promotes phasmids as a potentially high-value source of novel pcwdes for biotechnological applications. cellulases and pectinases are highly sought after by the biofuel industry to degrade feedstock into the monomers later converted into fuel, or to improve the flow rate of the material by reducing the amount of solid matter. pectinases are also used in the production of coffee, tea, and juice, and in waste-water treatment. phasmid pcwdes could be introduced into bacteria or fungi like trichoderma reesei, for industrial-scale enzyme production or direct use in bioreactors for wastewater treatment or biofuel production. our de novo midgut transcriptomes enabled us to survey all expressed pcwdes of the phasmatodea at once and identify conserved catalytic domains, justifying downstream translation and activity level analyses. a benefit of this is system is the increased speed and efficiency compared to the converse: running chemical assays to identify enzyme activity, using proteomics to identify the amino acid sequence of isolated enzymes, and working backwards from there to design a primer for the enzyme-encoding gene and hope it exists within the target organism’s genome itself and not a symbiont or contaminant. another benefit is that transcriptomics can reveal genes useful for phylogenetic analysis but that are not translated or whose proteins are modified post-translation such that the standard biochemical tests do not detect their function. an associated drawback is that expressed genes are not necessarily translated into active proteins, nor are they necessarily active at the site of expression. however, when a reference genome is not available, a transcriptome can provide large sets of potential genes for study and, combined with genomic data, can determine whether or not they are endogenous to the target organism, which cannot be determined by homology alone. de novo transcriptome assembly combined with rna-seq is a powerful tool for suggesting putative functions for unknown tissues in understudied organisms and directions for future study.  the folivorous phasmatodea are an ideal system to study the evolution of obligate herbivory, yet a paucity of genetic resources and poorly understood basic biology impede such work. using rna-seq, we demonstrated a diversity of plant cell wall degrading enzymes expressed differentially in the anterior section of the phasmid midgut. of these, the cellulases, cellobiases, and beta- -glucanases are likely all encoded in the insect’s own genome, as are the pectinases, though we could not definitively rule out a microbial source for the latter. such an abundance of endogenous enzymes was not expected from the polyneoptera, raising important questions on their evolutionary history. the efficiency by which our de novo transcriptomes generated new genomic resources and hypotheses for future research on polyneopteran digestion demonstrate the power of such methods to analyze organisms lacking sequenced genomes. our findings strongly encourage expanding the searches for pcwdes, most notably the pectinases and beta- -glucanases, into other, lower polyneopteran insects. availability of supporting data all reads and sequence files described in the manuscript are available under bioproject accessions prjna for p. schultei and prjna for the other phasmids. electronic supplementary material additional file 1: table s1: kegg table for the most highly expressed genes of the aretaon asperrimus midgut. top most highly expressed genes in each of the anterior midgut and posterior midgut used. additional file 2: table s2: kegg table for the most highly expressed genes of the extatosoma tiaratum midgut. top most highly expressed genes in each of the anterior midgut and posterior midgut used. additional file 3: table s3: kegg table for the most highly expressed genes of the medauroidea extradentata midgut. top most highly expressed genes in each of the anterior midgut and posterior midgut used. additional file 4: table s4: kegg table for the most highly expressed genes of the peruphasma schultei midgut. top most highly expressed genes in each of the anterior midgut and posterior midgut used. additional file 5: table s5: kegg table for the most highly expressed genes of the ramulus artemis midgut. top most highly expressed genes in each of the anterior midgut and posterior midgut used. additional file 6: table s6: kegg table for the most highly expressed genes of the sipyloidea sipylus midgut. top most highly expressed genes in each of the anterior midgut only. additional file 7: table s7: representative isotigs for each pcwde isogroup for the phasmatodea. data is from the full midgut transcriptomes with short sequences removed. transcripts were identified as pcwdes based on amino acid alignment to known proteins from the ncbi database. * = the sequence contained a complete open-reading frame. other sequences could have been truncated at the 5’ or 3’ or both. † = the sequence or another in that isogroup had a poly-a tail. additional file 8: table s8: – species and ncbi accession no's for cellulase proteins compared to phasmid proteins. sequences with *were not included in the alignment figure due to poor or absent overlap with other sequences at that region. table s – species and ncbi accession no's for pectinase proteins compared to phasmid proteins. sequences with *were not included in the alignment figure due to poor or absent overlap with other sequences at that region. table s – species and ncbi accession no's for cellobiase proteins compared to phasmid proteins. sequences with *were not included in the alignment figure due to poor or absent overlap with other sequences at that region. table s – species and ncbi accession no's for beta- -glucanase proteins compared to phasmid proteins. sequences with *were not included in the alignment figure figure 9) due to poor or absent overlap with other sequences at that region. additional file 9: table s9: number of peruphasma schultei midgut transcriptome sequences with successful blast search, mapping, and annotation. additional file 10: figure s1: species distribution for top-hit blast results of p. schultei midgut transcriptome. additional file 11: table s12: the most differentially expressed genes in the p. schultei amg. includes genes found in both or only one tissue type. means measured in rpkm. ppde = posterior probability of differential expression. annotations made with blast2go. additional file 12: table s13: the most differentially expressed genes in the p. schultei pmg. includes genes found in both or only one tissue type. means measured in rpkm. ppde = posterior probability of differential expression. annotations made with blast2go. additional file 13: table s14: the most highly expressed genes of the p. schultei amg. identified as genes with expression levels ten times greater than the mean for that section. means measured in rpkm. annotations made with blast2go. additional file 14: table s15: the most highly expressed genes of the p. schultei pmg. identified as genes with expression levels ten times greater than the mean for that section. means measured in rpkm. annotations made with blast2go. abbreviations amganterior midgut ciprescyberinfrastructure for phylogenetic research degdifferentially expressed genes ecenzyme commission ghglycoside hydrolase gogene ontology heghighly expressed gene keggkyoto encyclopedia of genes and genomes ncbinational center for biotechnology information pcwdeplant cell wall degrading enzyme pdbprotein data bank pmgposterior midgut ppdeposterior probability of differential expression rpkmreads per kilobase per million rsemrna-seq by expectation maximization srasequence read archive xsedeextreme science and engineering discovery environment. competing interests the authors declare that they have no competing interests. authors’ contributions authors ms and brj conceived the study and contributed equally to the work. lsk and brj provided resources and expertise. ms and lsk reared the insects. wcj, ja, and brj carried out the rna seq and de novo transcriptome assembly. ms, wcj, ja, and brj analyzed and interpreted data. all authors read and approved the final manuscript. 
the latter’s assembly yielded > transcripts, with over transcripts uniquely or more highly expressed in specific midgut sections. the differential gene expression between the two midgut sections provides the first molecular hints as to their function in living phasmids. if these expressed genes are translated into active enzymes, then phasmids can theoretically break plant cell walls into their monomer components independently of microbial symbionts. sequence alignments revealed catalytic sites on most pcwde transcripts. they represent an ideal system for investigating the adaptations necessary for obligate folivory, including plant cell wall degrading enzymes. we de novo assembled transcriptomes for the anterior and posterior midguts of six diverse phasmatodea species, with rna-seq on one exemplar species, peruphasma schultei. while most phasmid pcwde genes showed homology with those of other insects, the pectinases were homologous to bacterial genes. stick and leaf insects are an exclusively leaf-feeding order of insects with no record of omnivory, unlike other “herbivorous” polyneoptera. these genes were also found in genomic dna from phasmid brain tissue, suggesting endogenous production. two to three dozen pcwde encoding gene families, including cellulases and pectinases, were differentially expressed in the anterior midgut.
 stick and leaf insects are an exclusively leaf-feeding order of insects with no record of omnivory, unlike other “herbivorous” polyneoptera. they represent an ideal system for investigating the adaptations necessary for obligate folivory, including plant cell wall degrading enzymes. however, their physiology and internal anatomy is poorly understood, with limited genomic resources available.  we de novo assembled transcriptomes for the anterior and posterior midguts of six diverse phasmatodea species, with rna-seq on one exemplar species, peruphasma schultei. the latter’s assembly yielded > transcripts, with over transcripts uniquely or more highly expressed in specific midgut sections. two to three dozen pcwde encoding gene families, including cellulases and pectinases, were differentially expressed in the anterior midgut. these genes were also found in genomic dna from phasmid brain tissue, suggesting endogenous production. sequence alignments revealed catalytic sites on most pcwde transcripts. while most phasmid pcwde genes showed homology with those of other insects, the pectinases were homologous to bacterial genes.  we identified a large and diverse pcwde repertoire endogenous to the phasmids. if these expressed genes are translated into active enzymes, then phasmids can theoretically break plant cell walls into their monomer components independently of microbial symbionts. the differential gene expression between the two midgut sections provides the first molecular hints as to their function in living phasmids. our work expands the resources available for industrial applications of animal-derived pcwdes, and facilitates evolutionary analysis of lower polyneopteran digestive enzymes, including the pectinases whose origin in phasmatodea may have been a horizontal transfer event from bacteria. electronic supplementary material the online version of this article contains supplementary material, which is available to authorized users. issue-copyright-statement© the author 2014 