Parallel molecular routes to cold adaptation in eight genera of New Zealand stick insects

The acquisition of physiological strategies to tolerate novel thermal conditions allows organisms to exploit new environments. As a result, thermal tolerance is a key determinant of the global distribution of biodiversity, yet the constraints on its evolution are not well understood. Here we investigate parallel evolution of cold tolerance in New Zealand stick insects, an endemic radiation containing three montane-occurring species. Using a phylogeny constructed from 274 orthologous genes, we show that stick insects have independently colonized montane environments at least twice. We compare supercooling point and survival of internal ice formation among ten species from eight genera, and identify both freeze tolerance and freeze avoidance in separate montane lineages. Freeze tolerance is also verified in both lowland and montane populations of a single, geographically widespread, species. Transcriptome sequencing following cold shock identifies a set of structural cuticular genes that are both differentially regulated and under positive sequence selection in each species. However, while cuticular proteins in general are associated with cold shock across the phylogeny, the specific genes at play differ among species. Thus, while processes related to cuticular structure are consistently associated with adaptation for cold, this may not be the consequence of shared ancestral genetic constraints.


28
These supplementary methods include additional information to the methods detailed in the 29 main article. Citation numbers refer to the main article, and additional citation are named in 30 text and listed alphabetically at the end of this section (Supplementary References). 31 ad libitum. Survival was assessed after seven days and used to indicate whether individuals 78 could survive a small amount of internal ice formation. 79 Because it can take several hours for ice formation to reach equilibrium in insects 51 , 80 we determined whether individuals of each species could survive 6 h of the progression of  To calculate the SCP for each locale/species, we used both the specific tests to were placed in a +21 ˚C incubator and held for 10-15 minutes prior to cooling to -5 ˚C at Sequence clean-up and assembly Prior to assembly, raw Illumina reads were filtered to retain only high quality sequences in 129 several steps. To remove sequences more likely to contain errors, all reads containing an 'N' 130 were discarded using ShortRead v1.16.3 55 implemented in R. Low quality (phred score < 30) 131 and Illumina primer sequences were then trimmed using Cutadapt v1.1 53 . Lastly, polyT 132 stretches on sequence ends were trimmed using PRINSEQ lite v 0.6 54 and remaining reads 133 shorter than 25bp were discarded. The quality of each individual library before and after 134 clean-up was verified using FastQC v 0.10.1 56 .

135
For each species, de novo assembly was conducted separately using Trinity  Tables 5). 143 To reduce an apparent overabundance of called splice variants, the resulting 144 transcripts were clustered using CD-HIT-EST 58 , which combined all contigs with > 95% 145 similarity. During library preparation, removal of ribosomal RNA could be both incomplete 146 and unequal among samples. To remove this bias among libraries, ribosomal sequences were 147 removed by Ribopicker online 59 . However, not all sequences were removed and in some 148 cases 12S, 16S and 28S were manually identified by blast and removed from the assembly.

149
The resulting number of contigs, as well as their N50 and N90, are listed in Supplemental 150 Table 5.
Assemblies were annotated by a BlastX 60 sequence search (e-value threshold 1e -10 ) against the National Centre for Biotechnology Information (NCBI) non-redundant (nr) 153 database. Approximately 10% of all contigs were annotated this way, but that number was 154 higher (~25%) when only contigs > 200 bp were retained. However, to maximize 155 comparisons among datasets, 100-200bp long contigs were retained in all analyses. Gene 156 ontology (GO) terms were assigned using Blast2GO 38 based on the BlastX results.

157
To determine the overall similarity among the transcripts sequenced across species, 158 pairwise blasts with 1e -10 e-value threshold were conducted among each of the Trinity 159 alignments. Based on this, all alignments shared more than 99% of their transcripts with one  Orthologs used in phylogenetic analyses were extracted from each species' Trinity assembly 166 using HaMStr 61 , which was kindly run by I. Ebersberger using a core ortholog set containing 167 2,423 orthologs built based on reference genomes in six arthropod taxa: Daphnia pulex, 168 Acyrthosiphon pisum, Pediculus humanus, Bombyx mori, Apis mellifera, and Drosophila 169 melanogaster. Recovery of orthologs was high across taxa; 2,336 of core orthologs were 170 detected in at least one dataset. Of these, 1,485 were found in all 11 datasets. However, many 171 of these genes were not full length across all taxa, so we retained only genes which were 172 missing less than 2% of their nucleotides (due to both gaps and missing data). The remaining 173 genes were aligned using MAFFT v.7.047beta 62 with a maximum of 1000 iterations and the -