Phylogenetic relationships among botryllid species

200 Anchored Hybrid Enrichment (AHE) loci were sequenced in 55 specimens. The mean length of the loci is 705 bp and the median length is 296 bp, indicating that the distribution of locus length is skewed to the right. The range of locus length is between 121 and 5172 bp. The total length of the alignment is 141,107 bp.

Figure 1 is the tree generated by ASTRAL, and Supplementary Fig. ^S3 is the tree generated by RAxML. The topologies of the two trees are identical. The Botryllus genus is paraphyletic with respect to the Botrylloides genus. Focusing on the Botryllus genus first, there are three Botryllus-only clades (a, c, and d in Fig. 1), Clade a comprising Botryllus sp. (Bp1), and two specimens from the Western Pacific (the Philippines and Australia: Bsp1 and Bsp2). Clade c includes a Botryllus specimen from Florida (Bsp3) and a specimen of Botryllus horridus from Japan (Bh1). Clade d is formed by Botryllus gaiae (Bga1, Bga2), and Botryllus schlosseri from Italy (Bsp4), and is a sister group to Clade c. A specimen we collected in Bocas del Drago (Bsp5), Panama (in the Bocas del Toro archipelago), is the only Botryllus in a clade (e) that includes all the Botrylloides species. We will refer to this specimen as Bocas del Drago.

Open in a separate window

Figure 1

Maximum Quartet Support Species Tree of relationships among botryllid species, using the program Accurate Species TRee ALgorithm (ASTRAL). Nodes with less than 50% posterior probability support have been collapsed. The scale bar indicates two coalescent units. Red text=Genus Botryllus, Blue text=Genus Botrylloides. Colored circles refer to geographic regions (Orange=Western Pacific/Hawaii, Green=Eastern Pacific, Purple=Western Atlantic/Caribbean, Pink=English Channel/Mediterranean).

Moving on to the Botrylloides portion of the phylogeny, Clade g contains Botrylloides fuscus (Bf1–Bf3 in Clade h) as a sister taxon to a Botrylloides giganteus/Botrylloides violaceus group (Clade j). Clade r contains three specimens (Bsp7–9) of a species we are calling Rabbit Key in this manuscript, based on its collection location offshore of Rabbit Key in the Florida Keys. This species has also been found offshore of Barnes Key in the Florida Keys (the location of the sample in the phylogenomic trees). Clade r is a sister clade to Clade q, containing a Botrylloides sp. from the Bahamas (Bsp6). Clade r+q is a sister group to Clade p, containing specimens (Bsp10–13) from a species found in Bocas del Toro, Panama, which we will refer to as Bocas del Toro. Clades p, q, and r form Clade o, which is a sister clade to Clade n, a species from the Verde Island Passage in the Philippines which is represented by Bsp16–Bsp23. Clade s includes a specimen from the Pacific (Sydney, Australia: Bsp14) as a sister group to a clade labeled Botrylloides diegensis (v)+a clade labeled Botrylloides niger (u). The clade labeled Botrylloides diegensis (v) contains specimens with three different names: Botrylloides diegensis from California (Bd4 and Bd5), Botrylloides leachii from New Zealand (Bd6–8), and Botrylloides praelongus from Japan (Bd1–3). For brevity’s sake, we labeled Botrylloides diegensis/leachii/praelongus samples as Botrylloides diegensis in Fig. 1, Supplementary Figs. ^S3 and ^S4.

Sample identification: molecular techniques

We assigned samples to species by sequencing the mitochondrial cytochrome oxidase I (mtCOI) gene (Table ​(Table1).1). DNA was extracted using the Nucleospin Tissue Kit (Macherey Nagel). DNA was initially extracted from pieces of whole colony for each sample. In cases where extracted DNA failed to amplify, DNA was then extracted from zooids that were dissected from the colony. PCR amplification was performed using either OneTaq DNA Polymerase (New England Biolabs) or Phusion High-Fidelity DNA Polymerase (New England Biolabs). OneTaq reactions were as follows: 20 µl total reaction volume with 2 mM MgCl₂, 0.2 mM dNTPs, 2 µl of 10× buffer, 0.2 mM of each primer, and 0.16 U of OneTaq. Phusion reactions were as follows: 20 µl total reaction volume with 4 µl HF buffer, 0.2 mM dNTPs, 0.6 µl of 100% DMSO, 0.2 U of Phusion. The amount of water and template DNA was individually determined for each PCR reaction, based on the concentration of template DNA in the sample: at least 30 ng of DNA was added to each PCR reaction.

Each DNA sample was amplified with one of two PCR primer pairs: Tun_forward/Tun_reverse2⁵⁸, or LCO1490/HC02198⁵⁹. Tun primers were only used with OneTaq polymerase, using this protocol: 94 °C for 1 min, 60× (94 °C for 10 s, 50 °C for 30 s, 72 °C for 50 s), 72 °C for 10 min. Folmer primers were only used with Phusion polymerase, using this protocol: 98 °C for 30 s, 35× (98 °C for 10 s, 48 °C for 30 s, 72 °C for 30 s), 72 °C for 5 min. PCR products were incubated with 1 µl each of Exonuclease I (New England Biolabs) and Antarctic Phosphatase (New England Biolabs) at 37 °C for 1 h, followed by 90 °C for 10 min. The PCR products were sequenced at the University of Kentucky's HealthCare Genomics Core Laboratory using an ABI-3730 automated sequencer (Applied Biosystems). Forward and reverse sequences were edited and combined into a consensus sequence using Codon Codes Aligner (Codon Code Corporation).

We compared each sequence to botryllid sequences available on GenBank. If the sequence we obtained had 98–100% identity to a sequence identified on GenBank using blastn, we considered our sample to be the same species as the sample on GenBank. Because GenBank sequences can be mis-assigned, we only used GenBank identifications in which the submitting author had independently verified the taxonomic assignment of the sample using morphological characters. In many cases, the mtCOI sequence had no close match on GenBank.

Sample identification: morphological techniques

Taxa were also assigned to species by examination of morphological characters, when formalin preserved samples were available, using descriptions from the literature^{11,30–36,38,40,43–48,60–65}. The 31 analyzed morphological characters are summarized as follows: arrangement of systems, position of ovaries and testes, testes morphology, number of stigmatal rows, completeness of the second stigmatal row, arrangement of stigmata, shape of intestine, location of anterior edge of intestinal loop, location of anus, number of stomach folds, appearance of the stomach folds, shape of the stomach, shape and size of the pyloric caecum, number and size orders of the oral tentacles, distribution of pigment cells in the zooid, zooid length, colony color when living and after fixation, and tunic thickness^30–32. If the species could not be identified with complete certainty from the literature, we used a subset of the 31 characters: arrangement of systems, position of ovaries and testes (if ovaries and testes were present), appearance of the stomach folds, shape of the stomach, and size of the pyloric caecum. At least 30 zooids were examined from each colony.

When possible, we examined the morphology of the same colony from which the mtCOI sequencing was obtained. If a formalin preserved sample of the original colony was not available, we were often able to examine the morphology of another colony with an identical mtCOI sequence. We did not assign a specimen to a specific morphologically-described species if we could not obtain a formalin sample of the original colony or of a colony with an identical mtCOI sequence.

In the course of our species assignment using morphology, several specimens did not match the descriptions of any species previously found in the geographic region in which the specimen was collected. Either these samples represent new species, or they are known species that have only recently been identified from the location where they were collected. To determine where these species fit in relation to the morphologies of the other Botryllus/Botrylloides species, we searched the literature for morphological information on the 53 described botryllid species. We compiled data on the 31 morphological characters described above. The genus names (Botrylloides or Botryllus) and the morphological data come from the type description. In a few cases, the type description is lacking data on the majority of these characters. If a re-description was available, it was used to supplement the type description. Because our examinations were more thorough than type descriptions (often type descriptions do not provide data on individual zooids or individual colonies), we averaged our quantitative data across the 30+zooids we examined in each colony to obtain a single value for each colony. We then averaged quantitative data across multiple colonies to obtain a single value for each character for each species, to match the data available in the literature. A brief description of each considered morphological character is available in Supplementary Fig. ^S2, and the entire data matrix is available in Supplementary Table ^S1. Using the morphological data from the 53 described species and the 4 undescribed species, we then conducted a Principal Components Analysis using PCAmixdata⁶⁶ as implemented in R version 3.6.1.

To accompany the PCA, we constructed a phylogenetic tree of the 57 species using the 31 morphological characters using MrBayes 3.2.2⁶⁷ on the CIPRES (Cyberinfrastructure for Phylogenic Research) Science Gateway⁶⁸. The GTR+G model of nucleotide substitution was applied to all data sets (Nset=6). Each analysis was run for 10 million generations, with sampling every 1000 generations. The first 2000 trees were eliminated as burn-in.

Anchored hybrid enrichment (AHE) locus identification and probe design

Our aim was to develop a resource for collecting hundreds of orthologous loci across the botryllid ascidians using Anchored Hybrid Enrichment (AHE)⁶⁹. The pre-existing genomic resources included an assembled genome of Botryllus schlosseri⁷⁰, and two assembled transcriptomes: Botryllus schlosseri⁷¹, and Botrylloides leachii⁷², recently re-assigned to Botrylloides diegensis in Ref.²⁴. In order to better represent the high diversity of the botryllid group, we collected low-coverage, whole genome data assemblies for seven additional species (details are given in Supplementary Table ^S2). DNA extracts for these seven species were sent to the Center for Anchored Phylogenomics (http://www.anchoredphylogeny.com) for processing. In brief, after the quality/quantity of DNA was assessed using Qubit, Illumina libraries with single 8 bp indexes were prepared following⁷³, with modifications described in Ref.⁷⁴. Libraries were pooled and sequenced on two Illumina HiSeq2500 lanes with a paired-end 150 bp protocol. A total of 125 Gb of data was collected yielding 25–65×coverage per species. Reads were filtered for quality using the Cassava high chastity filter, demultiplexed with no mismatches tolerated, and merged to remove sequence adapters⁷⁵ prior to downstream processing.

In order to identify suitable conserved targets for AHE, we performed reciprocal blast on local machines at the Center for Anchored Phylogenomics using the two assembled transcriptomes (blastn). Using the results from the blast searches, we identified 482 preliminary targets with matching transcripts, which we aligned using MAFFT v7.023b⁷⁶. Alignments were manually inspected in Geneious (vR9, Biomatters Ltd., Kearse et al. 2012), then trimmed to regions that were well-aligned. For the remainder of the locus development/identification, we followed the protocol outlined in Ref.⁷⁷. More specifically, we isolated the Botryllus schlosseri (transcriptome) sequences from the aforementioned alignments, and using those as a reference scanned the Botryllus schlosseri genome for the AHE regions. Regions of 10,000 bp containing a 17 of 20 initial spaced k-mer match, followed by a 55 of 100 confirmation match to one of the references were kept. K-mers are all of a sequence’s subsequences of length=k. For example, the sequence GCTA would have the following k-mers: G, C, T, A, GC, CT, TA, GCT, CTA, and GCTA. K-mers from the Botryllus schlosseri transcriptome were used to search the Botryllus schlosseri genome for AHE regions, and matches were based on spaced seeds as described in Ref.⁷⁸. We then aligned (using MAFFT), the best matching genome sequence for each locus to the two transcriptome-derived sequences for that locus. Using Geneious (vR9, Biomatters Ltd.), we identified well-aligned regions of each three-sequence alignment and trimmed the alignments accordingly. The three-sequence alignment contained only two species: Botryllus schlosseri and Botrylloides leachii (recently re-assigned as Botrylloides diegensis)²⁴.

In order to incorporate whole genome sequencing (WGS) data from the seven additional species, we utilized sequences from Botrylloides leachii and Botryllus schlosseri in the alignments as references. Each WGS read was checked against the reference database and reads with a preliminary 17 of 20 initial spaced k-mer match, followed by a final 55 of 100 bp consecutive match were retained, then aligned by locus to form seeds for an extension assembly that allowed flanking regions to be recovered (see Ref.⁷⁷ for details and scripts). In order to construct the final alignments, the (up to) 10 sequences for each locus were aligned in MAFFT, then trimmed to well-aligned regions after inspection in Geneious (vR9, Biomatters Ltd.). In order to avoid problems associated with missing data in downstream projects⁷⁹, loci represented by less than 50% of the sequences in the alignment were removed from downstream analysis. When alignments from two loci were found to be overlapping (i.e. containing some of the same 20-mers), one locus was removed to ensure that each locus was a unique target. Lastly, we checked for repetitive elements by profiling the k-mers found in the alignments with respect to their occurrence in the WGS reads. Regions with a substantially elevated k-mer coverage were masked. A total of 200 AHE targets resulted from the process. Supplementary Table ^S3 contains the size of each locus, and genomic position of each locus in the Botryllus schlosseri genome, as determined from the best blastn match to the Botryllus schlosseri genome assembly using the locus sequence as a query. Finally, in-silico probes were tiled uniformly across the 10 sequences for each locus at 3.5×coverage depth. A total of 54,350 probes covered the 200 AHE targets (total target size~139 kb) that resulted from the process. These loci were successfully amplified in Symplegma brakenhielmi, to provide an outgroup for the phylogenomic tree. These loci will therefore be useful for Symplegma, which is the sister group to the Botrylloides/Botryllus clade⁸⁰. The utility of these loci beyond the genera Botrylloides, Botryllus and Symplegma has not been investigated.

DNA extraction and library preparation methods

DNA for all samples in Table ​Table11 was extracted using the E.Z.N.A DNA isolation kit (Omega BioTek), and an additional isopropanol precipitation was performed to further purify the DNA. The quantity and quality of the DNA extractions were determined using Qubit and 2% TAE agarose gels. The extracted DNA was fragmented into 300–500 bp pieces using a Covaris E220 focused-ultrasonicator with microTUBES (Covaris). Then, library preparation and indexing were performed on a Beckman-Coulter Biomek FXp liquid-handling robot, using a protocol based on Ref.⁷³. Anchored hybrid enrichment was performed using a custom SureSelect kit (Agilent Technologies) targeting loci designed from the whole genome alignment. Sequencing data were generated on an Illumina HiSeq2500 platform at the Center for Anchored Phylogenomics at Florida State University (www.anchoredphylogeny.com), as in Ref.⁷⁹. Sequencing was performed in the Translational Science Laboratory in the College of Medicine at Florida State University.

Raw read alignment

Sequence reads were demultiplexed with no mismatches tolerated and filtered for quality using the Illumina CASAVA pipeline with a high chastity setting. Overlapping reads were identified and merged using the approach described by Ref.⁷⁵. This process removes sequence adapters and corrects sequencing errors in overlapping regions. Reads were then assembled using the quasi-de novo approach described by Ref.⁷⁷. This assembly approach uses divergent references to identify sequences coming from conserved regions to which reads can be mapped. The mapped reads are in turn used as references when the assembly is extended into less conserved regions (see Ref.⁷⁷ for details). Probe region sequences from eight of the nine species used in the probe design were used as references for the initial mapping, while sequences from the Botryllus schlosseri genome (the 9th species) served as the primary reference. Consensus sequences were constructed from assembly clusters containing greater than an average of 250 reads. Ambiguity codes were employed for sites in which base frequencies could not be explained by a 1% sequencing error.

We would like to thank the following people who generously collected, shared, or facilitated sampling for this study: Carrie Craig, Derrick Cruz, Tony De Tomaso, Tom Frankovich, Erin Grey, Vanessa Guerra, Kristen Larson, Gretchen Lambert, Beth Moore, Rosana Rocha, Marc Ruis, Greg Ruiz, Shinsuke Saito, Yasunori Saito, Joelle Tirindelli, Linda Walters, Terry Gosliner, Christina Piotrowski, Rich Mooi, and Gary Williams. We would also like to thank Stefano Tiozzo and Megan Wilson for sending us transcriptomes, Christina Piotrowski, California Academy of Sciences, for help with sample curation and sampling, Carmela Gissi and Federica Montesanto for improving the manuscript, and Todd Newberry and Rick Grosberg for inspiring us to study botryllids. We thank the following students for contributing to sample processing, characterization, and tracking: Verena Wang, Ritchelle Quiambao, Karen Alroy, Elizabeth Sheets, Ryan Fergusson, KeChaunte Johnson, Tim Fuller, and Kelly Donahoe. We acknowledge The Philippines’ Bureau of Fisheries and Aquatic Resources and the National Fisheries Research and Development Institute for their role in the 2014–2016 Verde Island Passage (VIP) Expedition. All specimens from the Philippines were collected under Gratuitous Permit # GP-0077-14 from the municipalities of Puerto Galera and Batangas. The following samples were collected during the NSF Biotic Survey VIP Expedition: CASIZ 215739, CASIZ 215761, CASIZ 209556, CASIZ 209557, CASIZ 201262, CASIZ 201265, CASIZ 201274, CASIZ 217839, CASIZ 204513. This work was funded by an NSF EPSCoR Research Enhancement Grant. Sample acquisition and some analysis also benefited from projects funded by the National Science Foundation including NSF FSML Grant 0435033 and NSF DBI 1257630, California Sea Grant, the Smithsonian, and the Indian River Lagoon National Estuary Program. This study is dedicated to Yasunori Saito.