Protocol used in Evidence supporting cryptic species complexes within two sessile microinvertebrates, Limnias melicerta and L. ceratophylli (Rotifera, Gnesiotrocha) v1

Sample collection and culture Aquatic plant samples were collected from habitats across the USA and a sediment sample from Australia (Tables S1.1, S1.2). Limnias melicerta and L. ceratophylli were identified and isolated from rehydrated sediments or by removing a piece of vegetation to which they were attached. Species identification was based on tube structure, shape of corona, antennae length, and the number of dorsal nodules [1]. Clonal lineages initiated from single females were cultured in modified MBL media [2] and fed a mixture of the algae Chlorella vulgaris Berijerinck, 1890 (The UTEX Culture Collection of Algae at the University of Texas at Austin [UTEX] strain 30) and Chlamydomonasreinhardtii Dangeard, 1888 (UTEX strain 90). Rotifers in this genus produce tubes of hardened secretions [3]; we added powdered carmine (Alfa Aesar, UK) to lab cultures to provide a supplementary matrix to aid tube construction and to increase their visibility in culture. Voucher specimens were deposited in the UTEP Biodiversity Collections at The University of Texas at El Paso (L. melicerta: UTEP:Zoo:43, 105-134; L. ceratophylli: UTEP:Zoo:32-42, UTEP:Zoo:52-61).Deposited specimens included approximately 10 individuals from each population preserved in 95% ethanol and 10 clonal individuals preserved in 4% buffered formalin for molecular analyses and identification, respectively. DNA Extraction and gene amplification DNA was extracted from one individual of each clonal lineage by adding 13 µl Chelex®-100 (Bio-Rad Laboratories, CA, USA) and incubating at 100°C for 10 min. DNA templates were stored at -80°C until used for amplification. Number of clonal lineages examined from each population is given in Tables S1.1 and S1.2. An approximate 630 bp portion of the cytochrome c oxidase subunit I (COI) gene was amplified using the primers LCO1490: 5' -GGTCAACAAATCATAAAGATATTGG-3' and HCO2198: 5'-TAAACTTCAGGGTGACCAAAAAATCA-3' [4]. The entire nuclear internal transcribed spacer region (ITS) was amplified using the primers ITS4: 5'-TCCTCCGCTTATTGATATGC-3' and ITS5: 5'-GGAAGTAAAAGTCGTAACAAGG-3' [5], and 865 bp of the 18S rRNA gene was amplified using primers 3F: 5’-GTTCGATTCCGGAGAGGG-3’ as modified by Giribet et al. [6] and primer 18Sbi: 5’-CTAGAGTCTCGTTCGTTATCGG-3’ as modified by Whiting et al. [7]. PCR reactions contained 10 µl of genomic DNA, 1 µl of each primer (500 ng/ µl), 22 µl HPLC grade sterile water, 1 µl GoTaq® G2 DNA Polymerase (Promega), 10 µl 5X PCR buffer B (10 mM MgCl2, pH 8.5) or 5X PCR buffer A (7.5 mM MgCl2, pH 8.5), followed by adding 5 µl dNTP mix (2.5 mM each of dATP, dCTP, dGTP, dTTP) at 80ºC. PCR cycles were run on a thermocycler (TECHNE TC-412) and consisted of an initial denaturation at 94ºC for 1 min, followed by denaturation at 94ºC for 1 min, annealing at 48ºC for 2 min and extension at 72ºC for 3 min for 35 cycles, and a final extension step at 72ºC for 7 min. To verify the size of amplification products we used electrophoresis, and we purified them using GENECLEAN® kits (MP Biomedicals, LLC) before sequencing. Sequencing was done at UTEP’s BBRC Genomic Analysis Core Facility on an Applied Biosystems 3130xl Genetic Analyzer using BigDye Terminator v3.1 Cycle Sequencing Kits (Applied Biosystems). GenBank accession numbers for all sequences obtained are given in Tables S1.1 (L. melicerta) and S1.2 (L. ceratophylli). The COI gene sequences of L. melicerta (accession number, KT870155.1) and L. ceratophylli(KT870157.1) from GenBank are not included in our analyses for two reasons. 1) The COI sequence of L. melicerta KT870155.1 is 330 bp, which was much shorter than COI sequences obtained in this study (623 bp). 2) The COI sequence of L. ceratophylli KT870157.1 grouped with cryptic species M of L. melicerta in phylogenetic analyses. Additional 18S rRNA sequences (L. melicerta: KM873599.1, L. ceratophylli: KM873598.1) and a COI sequence from L. melicerta (KT870154.1) were included from GenBank. Sinantherina socialis(Linneaus, 1758) and Ptygura pilula (Cubitt, 1872) were included as outgroups in phylogenetic analyses for the COI gene. Floscularia conifera (Hudson, 1886) andPtygura brachiata (Hudson, 1886) were used as outgroup taxa in phylogenetic analyses based on ITS region. For analysis of 18S rRNA sequences, we used Collotheca campanulata as the outgroup (Table S1.1). Genetic Diversity FinchTV v 1.4.0 [8] was used to check sequences manually, especially for potential double peaks in the ITS region sequences. The ITS region alignment was uploaded to the SeqPhase online tool (http://seqphase.mpg.de/seqphase/) to phase the sequences as described by Flot [9]. Contigs for all sequences were made using CAP 3 [10] and were aligned using MAFFT v 7 [11]. Mesquite v 3.2 [12] was used to manually check the alignments and to translate COI gene sequences to proteins. To measure substitution saturation, we used DAMBE v 6 [13]. Number of polymorphic sites, number of parsimony informative sites, number of haplotypes, haplotype diversity (h), and nucleotide diversity (π) were calculated using DnaSp v 5.10.01 [14], and uncorrected pairwise sequences distances ("p") were calculated in Mega v 7.0 [15]. A haplotype network was constructed using the median joining method in Network v 5.0.3 [16]. Species delimitation Models for sequence evolution were TPM2uF+I+G for the COI gene, TPM1uF+I for the ITS region, and JC for the 18S rRNA gene as determined using Jmodeltest2 [17,18] available at the Cipres Science Gateway 3.3 [19]. To construct the phylogenetic trees, Bayesian analysis was run for 107 generations with two parallel runs and a 25% burn in period using MrBayes v 3.2.6 on XSEDE high-throughput computing resources available at Cipres Science Gateway [19]. Phylogenetic analyses were implemented in BEAST and *BEAST [20] using GTR+I+G model of sequence evolution for the COI gene and GTR+I model for the ITS region. TPM2uF and TPM1uF models are not available in BEAST. However, both of these models are classified under the GTR model. Thus GTR was used in both instances. To determine the number of evolutionary entities (putative cryptic species), we used Generalized Mixed Yule Coalescent (GMYC, [21]), Poisson Tree Process (PTP, [22]), Automatic Barcoding Gap Discovery (ABGD, [23]), *BEAST v 1.8.3 [20], and Bayesian Species Delimitation (BSD) implemented in Bayesian Phylogenetics and Phylogeography software (BPP v 3.1, [24–26]). We used BEAST v 1.8.3 [20] to construct ultrametric trees and *BEAST v 1.8.3 [20] for species delimitation. Both analyses were run for the COI gene and ITS region sequences separately for 107 generations, with sampling every 1,000 generations. Tracer v 1.6.0 [27] was used to check the effective sample size (ESS>200) and to verify convergence. Consensus trees were obtained using TreeAnnotator v 1.8.3 with a 25% burn in. Ultrametric trees were used for species delimitation in single threshold and multiple threshold GMYC [28] (http://species.h-its.org/gmyc/, accessed June 12, 2018), Bayesian GMYC (bGMYC) [29] and PTP methods. bGMYC was run using the R package bGMYC v 1.0.2 for 100,000 iterations with sampling every 1,000 iterations. We ran PTP by uploading the ultrametric trees to the online tool available at http://species.h-its.org/ptp/(accessed June 12, 2018) andused default settings. ABGD delimitation was done by uploading the sequence alignment to the online tool available at wwwabi.snv.jussieu.fr/public/abgd/ under the default settings (accessed June 12, 2018). *BEAST v 1.8.3 [20] was run under assumptions regarding the number of species for both the COI gene and ITS region sequences. Lineages having posterior probabilities > 0.90 were retained in *BEAST analyses. To run BSD, the phylogenetic tree based on Bayesian inference was used as the guide tree, and we used the joint species delimitation and tree estimation method (unguided species delimitation) that does not rely on the topology of the guide tree. Isolation by Distance Geographic distance matrices were constructed using Geographic Distance Matrix Generator v 1.2.3 [30]. To test the correlation between genetic variation and geographic distances (log transformed; km) among populations, Mantel tests with 10,000 permutations were run using the R package ecodist v 1.2.9 [31]. Trophi morphology Trophi were prepared for scanning electron microscopy (SEM) by dissolving rotifer tissue in ~5% sodium hypochlorite, rinsing with deionized water 10-15 times, and air-drying on circular cover slips at room temperature [32]. Trophi were coated with gold/palladium using a Gatan 682 PECS sputter coater. SEM images were obtained at 20 kV using a Hitachi S-4800 system. Trophi were prepared for individuals from one clonal lineage from each examined population. We used a geometric morphometric approach to study variation in the shape and size of 92 trophi for L. melicerta and 60 trophi for L. ceratophylli among putative cryptic species. This method uses Cartesian coordinates for a set of anatomical landmarks [33]. SEM images were obtained from caudal and frontal views of the trophi. Using TPS series software [34], nine landmarks on the caudal view and 10 on the frontal view of the trophi were digitized (Fig 1). Configuration of landmarks were analyzed using Generalized Procrustes Analysis [33]. Trophi size was calculated as Centroid Size (CS): i.e., the square root of the sum of squared distances between landmarks and their centroid [35]. Variation in the shape of trophi based on landmarks was examined using Discriminant Analysis in SPSS v 24 [36]. Because trophi centroid size was not normally distributed, variation in size among putative cryptic species was tested using a non-parametric Kruskal-Wallis test, and between the two morphospecies using a non-parametric Man-Whitney U test implemented in SPSS v 24 [36]. References 1. Wallace RL, Kordbacheh A, Walsh EJ. Key to the currently recognized species of Limnias Schrank, 1803 (Rotifera, Monogononta, Gnesiotrocha, Flosculariidae). Zootaxa. 2018;4442: 307–318. doi.org/10.11646/zootaxa.4442.2.7. 2. Stemberger RS. A general approach to the culture of planktonic rotifers. Can J Fish Aquat Sci. 1981;38: 721–724. doi:10.1139/f81-095. 3. Wallace RL, Nogrady T, Snell TW, Ricci C. Rotifera 1: Biology, ecology and systematics. 2nd ed. Nogrady T, editor. Kenobi Productions, 2006, 293 pp. 4. Folmer O, Black M, Hoeh W, Lutz R, Vrijenhoek R. DNA primers for amplification of mitochondrial cytochrome coxidase subunit I from diverse metazoan invertebrates. Mol Mar Biol Biotechnol. 1994;3: 294–299. doi:10.1371/journal.pone.0013102. 5. White T, Bruns T, Lee S, Taylor J. Amplification and direct sequencing of fungal ribosomal RNA genes for phylogenetics. In: Innis M, Gelfand D, Shinsky J, White T, editors. PCR Protocols: A guide to methods and applications. Academic Press, 1990, pp. 315–322. citeulike-article-id:671166. 6. Giribet G, Carranza S, Baguna J, Riutort M, Ribera C. First molecular evidence for the existence of a Tardigrada + Arthropoda clade. Mol Biol Evol. 1996;13: 76–84. 7. Whiting MF, Carpenter JC, Wheeler QD, Wheeler WC. The Strepsiptera problem: phylogeny of the holometabolous insect orders inferred from 18S and 28S ribosomal DNA sequences and morphology. Syst Biol. 1997;46: 1–68. 8. Geospiza.com. FinchTV. Seattle, WA; 2014: http://www.geospiza.com/Products/finchtv.shtml. 9. Flot JF. seqphase: a web tool for interconverting phase input/output files and fasta sequence alignments. Mol Ecol Resour. 2010;10: 162–166. doi:10.1111/j.1755-0998.2009.02732.x. 10. Huang X, Madan A. CAP3: A DNA sequence assembly program. Genome Res. 1999;9: 868–877. 11. Katoh K, Standley DM. MAFFT. Multiple sequence alignment software version 7: improvements in performance and usability. Mol Biol Evol. 2013;30: 772–780. doi:10.1093/molbev/mst010. 12. Maddison WP, Maddidon DR. Mesquite: a modular system for evolutionary analysis. 2017: Version 3.2. http://mesquiteproject.org. 13. Xia X. DAMBE6: New tools for microbial genomics, phylogenetics and molecular evolution. J Hered. 2017;108: 431–437. 14. Librado P, Rozas J. DnaSP v5: a software for comprehensive analysis of DNA polymorphism data. Bioinformatics. 2009;25: 1451–1452. doi:10.1093/bioinformatics/btp187. 15. Kumar S, Stecher G, Tamura K. MEGA7: Molecular Evolutionary Genetics Analysis version 7.0 for bigger datasets. Mol Biol Evol. 2015;33: 1870-4. doi:10.1093/molbev/mst197. 16. Bandelt HJ, Forster P, Rohl A. Median-joining networks for inferring intraspecific phylogenies. Mol Biol Evol. 1999;16: 37–48. 17. Darriba D, Taboada GL, Doallo R, Posada D. jModelTest 2: more models, new heuristics and parallel computing. Nat Methods. 2012;9: 772. http://dx.doi.org/10.1038/nmeth.2109. 18. Guindon S, Gascuel O. A simple, fast, and accurate algorithm to estimate large phylogenies by maximum likelihood. Syst Biol. 2003;52: 696–704. doi:10.1080/10635150390235520. 19. Miller MA, Pfeiffer W, Schwartz T. Creating the CIPRES science gateway for inference of large phylogenetic trees. 2010. Gateway Computing Environments Workshop, 2010, pp. 1–8. doi:10.1109/GCE.2010.5676129. 20. Drummond AJ, Suchard MA, Xie D, Rambaut A. Bayesian phylogenetics with BEAUti and the BEAST 1.7. Mol Biol Evol. 2012;29: 1969–1973. doi:10.1093/molbev/mss075. 21. Pons J, Barraclough TG, Gomez-Zurita J, Cardoso A, Duran DP, Hazell S, et al. Sequence-based species delimitation for the DNA taxonomy of undescribed insects. Syst Biol. 2006;55: 595–609. doi:10.1080/10635150600852011. 22. Zhang J, Kapli P, Pavlidis P, Stamatakis A. A general species delimitation method with applications to phylogenetic placements. Bioinformatics. 2013;29: 2869–76. doi:10.1093/bioinformatics/btt499. 23. Puillandre N, Lambert A, Brouillet S, Achaz G. ABGD, Automatic Barcode Gap Discovery for primary species delimitation. Mol Ecol. 2012;21: 1864–1877. doi:10.1111/j.1365-294X.2011.05239.x. 24. Rannala B, Yang Z. Bayes estimation of species divergence times and ancestral population sizes using DNA sequences from multiple loci. Genetics. 2003;164: 1645–56. doi:10.1073/pnas.0913022107. 25. Yang Z, Rannala B. Bayesian species delimitation using multilocus sequence data. Proc Natl Acad Sci. 2010;107: 9264–9269. doi:10.1073/pnas.0913022107. 26. Yang Z, Rannala B. Unguided species delimitation using DNA sequence data from multiple loci. Mol Biol Evol. 2014;31: 3125–3135. doi:10.1093/molbev/msu279. 27. Rambaut A, Suchard MA, Xie D, Drummond AJ. Tracer v1.6. http://beast.bio.ed.ac.uk/Tracer. citeulike-article-id:13604975. 28. Fujisawa T, Barraclough TG. Delimiting species using single-locus data and the generalized mixed yule coalescent approach: a revised method and evaluation on simulated data sets. Syst Biol. 2013; 62: 707–724. doi:10.1093/sysbio/syt033. 29. Reid NM, Carstens BC. Phylogenetic estimation error can decrease the accuracy of species delimitation: a Bayesian implementation of the general mixed yule coalescent model. BMC Evol Biol. 2012;12: 196. doi:10.1186/1471-2148-12-196. 30. Ersts PJ. Geographic distance matrix generator version 1.2.3. Nat Hist. 2006; http://biodiversityinformatics.amnh.org/open_source/gdmg/documentation.php. 31. Goslee SC, Urban DL. The ecodist package for dissimilarity-based analysis of ecological data. J Stat Software. 2007;22: 1–19. https://www.jstatsoft.org/v022/i07. 32. Segers H. Rotifera of some lakes in the floodplain of the River Niger (Imo State, Nigeria). New species and other taxonomic considerations. Hydrobiologia. 1993;250: 39–61. doi:10.1007/BF00007494. 33. Adams DC, Rohlf FJ, Slice DE. Geometric morphometrics: Ten years of progress following the revolution. Ital J Zool. 2004;71: 5–16. doi:10.1080/11250000409356545. 34. Rohlf FJ. The tps series of software. Hystrix. 2015;26: 1–4. doi:10.4404/hystrix-26.1-11264. 35. Cavalcanti MJ, Monteiro LR, Lopes PRD. Landmark-based morphometric analysis in selected species of serranid fishes (Perciformes: Teleostei). Zool Stud. 1999;38: 287–294. 36. IBM Corp. IBM SPSS statistics for windows, version 24.0. IBM Corp. 2016.