Aligning evidence: concerns regarding multiple sequence alignments in estimating the phylogeny of the Nudibranchia suborder Doridina.

Molecular estimates of phylogenetic relationships rely heavily on multiple sequence alignment construction. There has been little consensus, however, on how to properly address issues pertaining to the alignment of variable regions. Here, we construct alignments from four commonly sequenced molecular markers (16S, 18S, 28S and cytochrome c oxidase subunit I) for the Nudibranchia using three different methodologies: (i) strict mathematical algorithm; (ii) exclusion of variable or divergent regions and (iii) manually curated, and examine how different alignment construction methods can affect phylogenetic signal and phylogenetic estimates for the suborder Doridina. Phylogenetic informativeness (PI) profiles suggest that the molecular markers tested lack the power to resolve relationships at the base of the Doridina, while being more robust at family-level classifications. This supports the lack of consistent resolution between the 19 families within the Doridina across all three alignments. Most of the 19 families were recovered as monophyletic, and instances of non-monophyletic families were consistently recovered between analyses. We conclude that the alignment of variable regions has some effect on phylogenetic estimates of the Doridina, but these effects can vary depending on the size and scope of the phylogenetic query and PI of molecular markers.

Introduction

The debate regarding multiple sequence alignment (MSA) construction [1-8] and molecular marker selection [9] for use in phylogenetic estimates has been well established in the literature. There is little consensus, however, on the most accurate and replicable approach when considering MSA assemblies [10]. The crux of the matter is how best to align and represent homology in highly variable regions of sequence data. Most approaches employ a combination of a mathematical algorithm and manual curation of variable regions. Unfortunately, in most instances manual curation is loosely defined—which makes subsequent studies difficult to replicate—and even though mathematical algorithms allow for studies to be easily replicated, they do not take into account evolutionary history and consider homology exclusively as similarity [3]. In addition, they often fail to accurately align variable regions [11].

The estimation of evolutionary relationships relies heavily on the ability to determine homology between any given set of characters. The wrongful designation of characters as homologous can confound interpretations of evolutionary relationships by creating homoplasies, which results in the loss of phylogenetic signal [12]. Thus, the determination of homology in either morphological character matrices or MSAs (computational or manually curated) has a direct impact on the accuracy and ability to replicate phylogenetic estimates. This has led to differing positions concerning variable regions and their importance in increasing [13,14] or decreasing [15-17] phylogenetic signal. These highly variable loop regions become increasingly more difficult to align as sequences become more divergent. This confounds the ability to determine accuracy in MSAs and may result in contradictory estimates of deeper evolutionary relationships. It has been proposed that the large sections created by highly variable loop regions in MSAs should be excluded because of the uncertainty in determining reliable estimates of positional homology [18]. Conversely, it has also been suggested that effects from MSAs should be explored in determining phylogenetic estimates due to their likelihood of containing phylogenetic signal [19,20].

Highly divergent sequences and a general lack of diverse molecular markers have made MSAs and resulting phylogenetic estimates for the Nudibranchia problematic [21-26]. For these reasons, we examine how MSA construction of common nDNA and mtDNA markers affect phylogenetic estimates of the diverse Nudibranchia clade Doridina. Dorid nudibranchs have been a robust model for investigating biochemical diversity [27-33], morphological evolution [34-42], colour evolution [43-49], population structure [50-52] and development [53-57]. Even though there have been previous phylogenetic estimates that included dorid nudibranchs (e.g. [35,36]), they did not include representatives of all 19 families within the Doridina, and thus, a comprehensive phylogenetic context for the evolution of these traits is deficient.

The Doridina suborder currently consists of 19 families and more than 2000 described species. The classification has been divided into five superfamilies (Bathydoridoidea, Doridoidea, Onchidoridoidea, Phyllidioidea and Polyceroidea) defined by morphological variation in gill and feeding structures [21,36,58,59]. Species that possess the ability to retract their gills into fully formed gill pockets, Eudoridoidea (=Cryptobranchia), are divided into two clades based on the presence or absence of radula, Labiostomata (=Doridoidea) and Porostomata (=Phyllidiodea), respectively [36]. Dorid nudibranchs that lack a fully formed gill pocket represent the Anthobranchia (=Phanerobranchia), which is subdivided into the Suctoria (=Onchidoridoidea) and the Non-Suctoria (=Polyceroidea) based on the presence or absence of a buccal pump. Even though morphological [36,59] and molecular [22,23,60,61] analyses have shown the Doridina to be monophyletic, no study has focused on a complete sampling of all currently recognized families. Thus, relationships between families are poorly understood due to conflicting phylogenies and inadequate sampling.

In this study, we use three common MSA methodologies: (i) strict mathematical algorithm; (ii) exclusion of variable or divergent regions and (iii) manually curated. We then discuss issues pertaining to the lack of resolution in phylogenetic estimates and limitations in resolving the base of the Doridina. We also address the possible effects alignment construction may have on phylogenetic informativeness (PI), and how phylogenetic signal may shift given alignment methodology. Lastly, we re-evaluate taxonomic classifications that were consistent between all three analyses.

Material and methods

Taxon sampling

We sampled (121 taxa) representatives from the five superfamilies and all 19 families that currently comprise the Doridina, and mined GenBank for specimens that had at least two sequences. To limit the amount of missing data, we obtained extractions or tissue samples from published specimens to sequence molecular markers that were previously unattained. Pleurobranchaea meckeli and Berthella martensi were chosen as outgroups because the Pleurobranchomorpha has been suggested to be sister to the Nudibranchia [23]. Additional members from the Aeolidina and Dendronotina were also used in testing the monophyly of the Doridina. Specimens, voucher numbers, collection sites and GenBank accession numbers are listed in table 1. Voucher specimens are located in the collections at California Academy of Sciences (CASIZ) and the Museum of National Scientific Center of Marine Biology, Vladivostok (MIMB).

Table 1.

Specimens successfully sequenced and used for molecular analyses. Voucher numbers, GenBank accession numbers and collection localities.

		GenBank accession number
specimen	voucher	16S	18S	28S	COI	locality
PLEUROBRANCHOIDEA
Pleurobranchaea meckeli (Blainville, 1825)		FJ917439	FJ917449	FJ917481	FJ917499	Mediterranean Sea, Spain
Pleurobranchidae
Berthella martensi (Pilsbry, 1896)	MZUCR 6982	HM162592	MF958319	MF958363	HM162683	Las Secas, Islas sin nombre, Panama
AEOLIDINA
Flabellinidae
Flabellina pedata (Montagu, 1815)		AF249247	AJ224788	—	AF249817	North Sea, Helgoland
DENDRONOTINA
Arminidae
Armina loveni (Bergh, 1866)		AF249243	AF249196	—	AF249781	North Sea, Kattegat
DORIDINA
Bathydoridoidea
Bathydorididae
Prodoris clavigera (Thiele, 1912)	CASIZ 167553	JX274067	MF958320	MF958364	JX274106	Elephant Island, Antarctica
Doridoidea
Actinocyclidae
Actinocyclus verrucosus Ehrenberg, 1831	CASIZ 189448	MF958311	MF958352	MF958397	MF958438	Kauai, Hawaii
Hallaxa indecora (Bergh, 1905)	CASIZ 179600	MF958302	MF958340	MF958386	—	South Loi Island, Marshall Islands
Hallaxa translucens Gosliner & S. Johnson, 1994	CASIZ 173447	EU982814	MF958341	MF958387	EU982760	Iles Radama, Madagascar
Hallaxa iju Gosliner & S. Johnson, 1994	CASIZ 175559	EU982813	—	—	EU982759	Maui, Hawaii
Cadlinidae
Aldisa sanguinea (J. G. Cooper, 1863)	CASIZ 182031	MF958309	MF958350	MF958394	MF958435	Marin Co., California
Aldisa sp.	CASIZ 175733	EU982818	MF958351	MF958395	MF958436	Tiger Reef, Malaysia
Cadlina cf. luteomarginata (MacFarland, 1905)	CASIZ 188599A	KJ653679	KP340317	KP340350	KM219678	Parksville, Vancouver Island, British Columbia
Cadlina laevis (Linnaeus, 1767)	CASIZ 175444	—	MF958359	MF958406	—	Scotland
Cadlina luarna (Marcus & Marcus, 1967)	CASIZ 175437	EU982768	—	—	EU982718	Guanacaste, Costa Rica
Cadlina modesta MacFarland, 1966	CASIZ 223286	MF958310	—	—	MF958437	Vista del Mar, San Luis Obispo Co., California
Cadlina pellucida (Risso, 1826)	CASIZ 175448	EU982774	—	MF958396	EU982724	Ilha de Pessegueiro, Baixo Alentejo Prov., Portugal
Cadlina sparsa (Odhner, 1921)	CASIZ 182932	EU982776	—	—	EU982726	La Jolla, San Diego Co., California
Chromodorididae
Cadlinella ornatissima (Risbec, 1928)	CASIZ 177420	MF958284	MF958325	MF958371	MF958415	Maricaban Island, Luzon, Philippines
Cadlinella ornatissima (Risbec, 1928)	CASIZ 175452	EU982779	—	—	EU982728	Mooloolaba, Queensland, Australia
Chromodoris alternata (Burn, 1957)	SAM D19281	AY458800	EF534031	—	EF535120	Port Phillip Bay, Victoria, Australia
Chromodoris ambigua Rudman, 1987	SAM D19260	AY458801	EF534038	—	EF535119	Port Phillip Bay, Victoria, Australia
Chromodoris quadricolor (Rüppell & Leuckart, 1830)		AF249241	AJ224773	—	AF249802	Red Sea, Egypt
Doriprismatica atromarginata (Cuvier, 1804)		—	AF249211	—	AF249789	Great Barrier Reef, Australia
'Felimare' elegans (Cantraine, 1835)		AF249238	AJ224779	—	AF249787	NE Atlantic, Spain
'Felimare' midatlantica (Gosliner, 1990)	CASIZ 175443	JQ727789	—	—	JQ727898	Islotes do Martinhal, Algarve, Portugal
'Felimare' picta verdensis Ortea, Valds & Garca-Gmez, 1996	CASIZ 179384	HM162594	MF958346	MF958389	HM162685	Gulf of Guinea, Ilha do Principe, Sao Tome & Principe
'Felimare' villafranca (Riaao, 1818)		AF249237	AJ224780	—	—	NE Atlantic, Spain
'Felimida' edmundsi Cervera, Garcia-Gmez & Ortea, 1989	CASIZ 179385	HM162595	MF958347	MF958390	HM162686	Gulf of Guinea, Ilha do Principe, Sao Tome & Principe
'Felimida' krohni (Verany, 1846)		AF249239	AJ224774	AY427445	AF249805	NE Atlantic, Spain
Goniobranchus geometricus (Risbec, 1928)	CASIZ 175549	JQ727717	—	—	JQ727842	Nosy Valiha, Iles Radama, Madagascar
Hypselodoris imperialis (Pease, 1860)	CASIZ 142952	EU982807	—	—	EU982754	Mala Wharf, Maui, Hawaii
Hypselodoris infucata (Rüppell & Leuckart, 1830)		FJ917427	FJ917442	FJ917467	FJ917485	NSW, Australia
Miamira magnifica Eliot, 1904	CASIZ 169951	EU982781	—	—	EU982731	Old Woman Island, Queensland, Australia
Thorunna daniellae (Kay and Young, 1969)	CASIZ 170055	EU982809	—	—	EU982756	Maalaea Marina, Maui, Hawaii
Tyrinna evelinae (Marcus, 1958)	CASIZ 175440	EU982811	—	MF958391	EU982757	Playa Ventana, Guanacaste, Costa Rica
Tyrinna nobilis Bergh, 1898	ZSM M20050508	EF534054	EF534035	—	EF535127
Dorididae
Aphelodoris sp. 1	CASIZ 176920	MF958293	MF958332	MF958379	MF958424	Oudekraal, Cape Prov., South Africa
Aphelodoris sp. 1		—	GQ326866	—	GQ292033	Auckland, New Zealand
Aphelodoris luctuosa (Cheeseman, 1882)		—	GQ326867	—	GQ292042	Auckland, New Zealand
Doris montereyensis (Cooper, 1862)	CASIZ 174493	MF958294	MF958333	—	MF958425	Battery Point, Crescent City, Del Norte Co., California
Doris odhneri MacFarland, 1966	CASIZ 188014	MF958295	MF958334	MF958380	—	Duxbury Reef, Marin Co., California
Doris sp. 8	CASIZ 192348	MF958306	MF958345	—	—	Red Sea, Saudi Arabia
Discodorididae
Asteronotus cespitosus (van Hasselt, 1824)	CASIZ 191096	MF958288	MF958328	MF958375	MF958419	Kranket Island, Papua New Guinea
Atagema cf osseosa (Kelaart, 1859)	CASIZ 185142	MF958296	MF958335	—	MF958426	Maui, Hawaii
Discodoris coerulescens Bergh, 1888	CASIZ 182850	MF958290	MF958330	MF958377	MF958421	Maricaban Island, Luzon, Philippines
Halgerda dalanghita Fahey & Gosliner, 1999	CASIZ 181264	MF958289	MF958329	MF958376	MF958420	Maricaban Island, Luzon, Philippines
Discodorididae
Peltodoris nobilis (MacFarland, 1905)	CASIZ 182223	EU982816	—	—	EU982761	Pillar Point, San Mateo Co., California
Platydoris sanguinea Bergh, 1905	CASIZ 177762	MF958285	MF958326	MF958372	MF958416	Maricaban Island, Batangaas Prov., Philippines
Rostanga calumus Rudman & Avern, 1989		FJ917427	—	FJ917468	FJ917485	NSW, Australia
Sclerodoris tuberculata Eliot, 1904	CASIZ 190788	MF958286	MF958327	MF958373	MF958417	Madang Prov., Papua New Guinea
Thordisa albomacula Chan & Gosliner, 2006	CASIZ 179590	MF958287	—	MF958374	MF958418	Kwajalein Atoll, Marshall Islands
ONCHIDORIDOIDEA
Akiodorididae
Armodoris anudeorum Valdés, Moran & Woods, 2011	LACM 3118	KP340290	GQ326879	KP340355	KP340387	McMurdo Sound, Ross Sea, Antarctica
Calycidorididae
Calycidoris guentheri Abraham, 1876	CASIZ 190966A	KP340301	KP340338	KP340371	KP340397	Chukchi Sea, Alaska
Diaphorodoris lirulatocauda Millen, 1985	CASIZ 184341	KP340307	KP340344	KP340377	KP340403	Duxbury Reef, Marin Co., California
Diaphorodoris luteocincta (Sars, 1870)	LACM 8.7A	KP340308	—	KP340378	KP340404	Bahia de Algeciras, Cadiz Prov., Spain
Diaphorodoris cf mitsuii	CASIZ 185986	KP340310	KP340345	KP340379	KP340406	Sepok Point, Philippines
Diaphorodoris papillata Portmann & Sandmeier, 1960	LACM 8.6A	KP340311	—	—	KP340407	Bahia de Algeciras, Cadiz Prov., Spain Peninsula
Corambidae
Corambe obscura (Verrill, 1870)	CASIZ 183942	KP340303	KP340340	KP340373	KP340399	New Castle Portsmouth Bay, New Hampshire
Corambe pacifica MacFarland & O'Donoghue, 1929	LACM 2007-2.6C	KP340305	KP340342	KP340375	KP340401	Long Beach Marina, Los Angeles Co., California
Corambe steinbergae (Lance, 1962)	CASIZ 190508	KP340306	—	—	KP340402	Pillar Point, San Mateo Co., California
Goniodorididae
Ancula gibbsoa (Risso, 1818)	CASIZ 182028	KP340291	KP340322	KP340356	KP340388	Cumberland Co., Maine
Ancula gibbsoa (Risso, 1818)	CASIZ 181211	MF958291	—	—	MF958422	Duxbury Reef, Marin Co., California
Goniodoris nodosa (Montagu, 1808)		AF249226	AJ224783	—	AF249788	NE Atlantic, Spain
Okenia kendi Gosliner, 2004	CASIZ 186125A	MF958297	—	MF958381	MF958427	Luzon, Philippines
Trapania reticulata Rudman, 1987	CASIZ 191431	MF958303	MF958342	—	MF958432	Tab Island, Papua New Guinea
Onchidorididae
Acanthodoris atrogriseata O'Donoghue, 1927	CASIZ 186000	KJ653646	KP340323	KP340357	KM219646	Puget Sound, Kitsap Co., Washington
Acanthodoris hudsoni MacFarland, 1905	CASIZ 179480	KJ653652	KP340324	KP340359	KM219650	Asilomar, Monterey Co., California
Acanthodoris nanaimoensis O'Donoghue, 1921	CASIZ 181569A	KJ653656	KP340325	KP340360	KM219657	Pillar Point, San Mateo Co., California
Acanthodoris pilosa (Abildgaard, 1789)	CASIZ 183941A	KJ653659	KP340326	KP340361	—	Cobscook Bay, Washington Co., Maine
Acanthodoris planca Fahey & Valdés, 2005	CASIZ 176116	KJ653669	KP340327	KP340362	KM219671	Table Bay, Western Cape Prov., South Africa
Acanthodoris rhodoceras Cockerell & Eliot, 1905	CASIZ 181572	KJ653671	KP340328	KP340363	KM219673	Pillar Point, San Mateo Co., California
Knoutsodonta brasiliensis (Alvim, Padula & Pimenta, 2011)	BNHS-Opistho-336	KC255225	—	—	KC255226
Knoutsodonta depressa (Alder & Hancock, 1842)	CASIZ 186769A	KP340315	KP340347	—	KP340409	Huelva, Spain
Knoutsodonta jannae (Millen, 1987)	CASIZ 175578	KP340296	KP340331	KP340366	KP340392	Pillar Point, San Mateo Co., California
Knoutsodonta oblonga (Alder & Hancock, 1845)	MN 3010A	—	KP340349	KP340385	KP340410	Mewstone, Skomer, United Kingdom
Onchidoris bilamellata (Linnaeus, 1767)	CASIZ 188593	KP340314	KP340346	KP340382	—	Puget Sound, Kitsap Co., Washington
Onchidoris muricata Muller, 1776	CASIZ 181312	KJ653676	KP340348	KP340383	KM219680	Asilomar, Monterey Co., California
Onchidoris evincta (Millen, 2006)	CASIZ 187758A	KP340293	KP340329	KP340364	KP340390	Puget Sound, Kitsap Co., Washington
Onchidoris macropompa Martynov, Korshunova, Sanamyan & Sanamyan, 2009	MIMB 34210	MF958292	MF958331	MF958378	MF958423	Avacha Bay, Kamchatka
Onchidoris proxima (Alder & Hancock, 1854)	CASIZ 183921A	KJ653673	KP340336	KP340369	KM219676	Passamaquody Bay Eastport, Washington Co., Maine
Onchidoris slavi Martynov, Korshunova, Sanamyan & Sanamyan, 2009	MIMB 34211	—	—	MF958409	MF958446	Avacha Bay, Kamchatka
Onchimira cavifer Martynov, Korshunova, Sanamyan & Sanamyan, 2009	MIMB 34209	MF958298	MF958336	MF958382	MF958428	Avacha Bay, Kamchatka
Onchimira cavifer Martynov, Korshunova, Sanamyan & Sanamyan, 2009		—	MF958360	MF958407	MF958445	Avacha Bay, Kamchatka
PHYLLIDIOIDEA
Dendrodorididae
Dendrodoris arborescens (Collingwood, 1881)	CMNH-ZM08965	—	AB917459	—	AB917436	Oka, Tateyama
Dendrodoris atromaculata (Alder & Hancock, 1864)	CASIZ 181231	MF958307	MF958348	MF958392	MF958434	Janao Bay, Luzon, Philippines
Dendrodoris denisoni (Angas, 1864)	CASIZ 177702	MF958308	MF958349	MF958393	—	Janao Bay, Luzon, Philippines
Dendrodoris fumata (Rüppell & Leuckart, 1830)	CASIZ 192304	—	MF958358	MF958405	MF958444	Red Sea, Saudi Arabia
Dendrodoris fumata (Rüppell & Leuckart, 1830)		—	AF249216	FJ917470	AF249799	Great Barrier Reef, Australia
Dendrodoris guttata (Odhner, 1917)	CMNH-ZM 08967	—	AB917461	—	AB917446	Okinoshima, Tateyama
Dendrodoris nigra (Stimpson, 1855)	CASIZ 182821	MF958318	MF958357	MF958404	MF958443	Maricaban Island, Luzon, Philippines
Dendrodoris nigra (Stimpson, 1855)		AF249242	AF249215	—	AF249795	Great Barrier Reef, Australia
Doriopsilla albopunctata (J. G. Cooper, 1863)	CPIC 00909	KR002428	—	—	KR002480	Long Beach, California
Doriopsilla bertschi Hoover, Lindsay, Goddard & Valdés, 2015	CPIC 01058	KR002462	—	—	KR002517	Bahía de los Ángeles, Baja California, Mexico
Doriopsilla davebehrensi Hoover, Lindsay, Goddard & Valdés, 2015	LACM 3419	KR002476	—	—	—	Bahía de los Ángeles, Baja California, Mexico
Doriopsilla fulva (MacFarland, 1905)	CPIC 00933	KR002444	—	—	KR002498	Malibu, California
Doriopsilla gemela Gosliner, Schaefer & Millen, 1999	CPIC 00938	KR002453	—	—	KR002506	Malibu, California
Doriopsilla janaina Marcus &. Marcus, 1967	CASIZ 173618	MF958312	MF958353	MF958398	—	Galápagos Islands, Ecuador
Doriopsilla miniata (Alder & Hancock, 1864)	CMNH ZM08970	—	AB917464	—	AB917457	Yoshio, Katsuura
Doriopsilla spaldingi Valdés & Behrens, 1998	CPIC 00908	KR002427	—	—	KR002479	San Pedro, Los Angele Co., California
Mandeliidae
Mandelia mirocornata Valdés & Gosliner, 1999	CASIZ 176266	MF958278	MF958321	MF958365	MF958411	Oudekraal, Cape Prov., South Africa
Phyllidiidae
Ceratophyllidia sp.	CASIZ 181247	MF958281	MF958323	MF958368	MF958413	Beatrice, Philippines
Phyllidia coelestis Bergh, 1905	CASIZ 190982	MF958279	—	MF958366	MF958412	Kranket Island, Madang Prov., Papua New Guinea
Phyllidiella nigra (van Hasselt, 1824)	CASIZ 186196A	MF958280	MF958322	MF958367	—	Maricaban Strait, Batangas Prov., Luzon, Philippines
Phyllidiella pustulosa (Cuvier, 1804)		AF249232	AF249208	—	—	Great Barrier Reef, Australia
Phyllidiopsis annae Brunckhorst, 1993	CASIZ 186138	MF958283	MF958324	MF958370	—	Philippines
Reticulidia halgerda Brunckhorst & Burn in Brunckhorst, 1990	CASIZ 186491	MF958282	—	MF958369	MF958414	Maricaban Island, Luzon, Philippines
POLYCEROIDEA
Aegiridae
Aegires albopunctatus MacFarland, 1905	CASIZ 182213	MF958313	MF958354	MF958399	MF958439	Marin Co., California
Aegires citrinus Pruvot-Fol, 1930	CASIZ 144027	MF958314	MF958355	MF958400	MF958440	Mooloolaba, Queensland, Australia
Aegires flores Fahey & Gosliner, 2004	CASIZ 191244	MF958316	—	MF958402	MF958442	Papua New Guinea
Aegires serenae (Gosliner and Behrens, 1997)	CASIZ 191285	MF958315	—	MF958401	MF958441	Papua New Guinea
Aegires villosus Farran, 1905	CASIZ 177563	MF958317	MF958356	MF958403	—	Maricaban Island, Luzon, Philippines
Gymnodorididae
Gymnodoris sp.	CASIZ 176781	—	MF958361	—	—	Pulau Penang, Malaysia
Hexabranchiidae
Hexabranchus sanguineus (Ruppell & Leuckart, 1828)	CASIZ 142942	MF958304	MF958343	—	—	Kapalua Bay, Maui, Hawaii
Hexabranchus sanguineus (Ruppell & Leuckart, 1828)	CAZIS 193381	MF958305	MF958344	MF958388	MF958433	Papua New Guinea
Okadaiidae
Vayssierea sp.	CASIZ 190731	—	MF958362	MF958408	—	Sunshine Coast, Kings Beach, Australia
Polyceridae
Kaloplocamus sp. 1	CASIZ 194412	MF958299	MF958337	MF958383	MF958429	South Madagascar, Madagascar
Limacia sp. 1	CASIZ 176312	HM162602	KP340320	KP340353	HM162692	False Bay, Western Cape Prov., South Africa
Limacia sp. 2	CASIZ 176276	HM162603	—	—	HM162693	Oudekraal, Cape Prov., South Africa
Nembrotha cristata Bergh, 1877	CASIZ 191428	MF958301	MF958339	MF958385	MF958431	Madang Prov., Papua New Guinea
Plocamopherus pecoso Valls and Gosliner, 2006	CASIZ 191587	MF958300	MF958338	MF958384	MF958430	Madang Prov., Papua New Guinea
Polycera quadrilineata (Müller, 1776)		AF249229	AJ224777	—	—	North Sea, Kattegat
Roboastra ricei Pola, Cervera & Gosliner, 2008	CASIZ 173900	HM162598	—	—	HM162688	Florida, 5 mi offshore of Loran Tower
Tambja marbellensis Schick & Cervera, 1998	CASIZ 180379	HM162599	—	—	HM162689	Setubal District, Outao, Portugal
Triopha catalinae (Cooper, 1863)	CASIZ 170648	HM162600	KP340321	KP340354	HM162690	Yacht Harbor, San Francisco, California
Triopha maculata MacFarland, 1905	CASIZ 181556	HM162601	—	—	HM181556	Duxbury Reef, Marin Co., California

Extraction, PCR protocols and sequencing

Genomic DNA was extracted using the Qiagen DNeasy Blood & Tissue extraction kit (Valencia, CA). Amplification and sequencing for targeted gene fragments 16S, 18S, 28S and cytochrome c oxidase subunit I (COI) followed protocols used by Hallas & Gosliner [62]. Amplified fragments were sequenced at the Center for Comparative Genomics located at the California Academy of Sciences and National Scientific Center of Marine Biology.

Multiple sequence alignments

We assembled and edited sequences in Geneious Pro v. 9.1.7 [63] and BioEdit [64]. We aligned rDNA fragments (16S, 18S and 28S) using three different methodologies to examine conflicts regarding estimated phylogenies. For our first method, we used a computer algorithm E-INS-i in MAFFT (MA) [65]. The E-INS-i algorithm is designed to handle sequence data with several conserved regions embedded in long unalignable gapped regions. Our second method excluded all variable regions from the initial MAFFT alignment of rDNA using Gblocks (GA) [11]. We determined blocks using a less stringent selection by allowing for smaller final blocks, gap positions within the final blocks and less strict flank positions. For our third method, we used the initial MAFFT alignment from method one, but manually curated variable regions (CA). This was done by hand to correct for possible inappropriate alignment of sequence regions. For each alignment method, we concatenated all four targeted molecular markers into single MSAs that resulted in three separate concatenated datasets.

Phylogenetic informativeness

We estimated PI for each molecular marker per alignment [9] through the PhyDesign web interface [66], and estimated each PI profile using the site-rates model HyPhy [67]. First, we generated ultrametric trees by converting our concatenated Bayesian inference (BI) phylogenetic estimates in Mesquite [68]. Owing to the lack of a fossil record for the Doridina, our ultrametric trees are not in known time units. However, it has been shown that PI profiles can be used effectively even if divergence time estimates are absent [69]. Then, we analysed alignments in relation to each of their resulting estimated phylogeny to gain a greater understanding of how PI can change based on the alignment approach.

Phylogenetic analyses

We determined evolutionary models using PartitionFinder v.1.1.1 [70], and partitioned our concatenated datasets by rDNA fragment as well as codon position for COI. We also analysed our nDNA and mtDNA fragments separately to investigate possible conflicting evolutionary histories. We did this by aligning both nDNA and mtDNA datasets using the algorithm E-INS-i in MAFFT, and used the same partitioning scheme as in our concatenated datasets.

We analysed all our datasets using BI and maximum-likelihood (ML). BI searchers were run using MrBayes v. 3.2.1 [71], and convergence was checked in TRACER v.1.5 [72]. The datasets were run for 5 × 107 generations with Markov chains sampled every 1000 generations. The standard 25% burn-in was calculated and remaining tree estimates were used to create a 50% majority rule consensus tree of calculated posterior probabilities. Posterior probabilities (pp) that exceeded 0.95 were considered strongly supported, and values 0.94 and below were interpreted as having low support. For our ML analyses, we calculated non-parametric bootstrapping (bs) values and the ML tree simultaneously in RAxML v. 7.2.6 [73]. We used the same partitioning scheme as in our BI search, but used the evolution model GTR+I− and executed fast bootstrapping runs for 5 × 104 iterations. Bootstrap values 70 or higher were considered strongly supported, while all other values were evaluated as weakly supported [74].

Results

Molecular data

Sequences obtained for phylogenetic analyses and PI profiles are labelled in table 1, and all alignments have been deposited in TreeBASE (http://purl.org/phylo/treebase/phylows/study/TB2:S20396). As expected, our concatenated alignments (MA, GA and CA) varied in length as well as number of parsimony informative characters (table 2). Surprisingly, there were more parsimony informative characters in the MA than the CA. The GA had far fewer parsimony informative characters than the other two alignments generated. Evolutionary model GTR+I+I− was selected by PartitionFinder v. 1.1.1 for each partition based on the Akaike information criterion. Interestingly, specific species of Dendrodoris (D. nigra, D. arborescens, D. guttata and D. fumata) had mtDNA regions that were highly divergent.

Table 2.

Summary of multiple sequence alignment variation between alignment methodologies for 16S, 18S, 28S and COI.

	length (base pairs)				parsimony informative characters
alignment	16S	18S	28S	COI	16S	18S	28S	COI
MAFFT	498	2578	960	661	315	485	388	370
Gblocks	390	1546	585	661	236	351	174	370
Curated	584	2722	1077	661	308	396	324	370

The net PI profiles depict the amount of signal through time of each molecular marker in relation to the respective estimated phylogeny (figure 1a–c). Regardless of the alignment, net PI profiles show most phylogenetic signal resides towards the tips of the trees. The amount of overall phylogenetic signal was lowest in our GA estimate, which suggests that information was lost when variable regions were excluded. Lastly, our MA appeared to have slightly more PI than our CA.

Figure 1.

Phylogenetic informativeness profiles estimated by PhyDesign. (a) MAFFT alignment. (b) Gblocks alignment. (c) Curated alignment.

Phylogeny

BI posterior probability estimates generally resulted in higher overall support of relationships than ML non-parametric bootstrap values. Branch lengths for basal nodes were extremely short (figure 2a–c), especially in regions of the tree that resulted in low resolution for relationships between superfamilies (figure 3a–c). There were, however, some general similarities between topologies of both phylogenetic searches and among alignments.

Figure 2.

Phylograms of Doridina phylogenetic estimates from all three different alignment constructions. Topologies represent Bayesian estimates. (a) MAFFT alignment. (b) Gblocks alignment. (c) Curated alignment.

Figure 3.

Cladogram of Doridina phylogenetic estimates. Topology represents Bayesian estimate. Branches are coloured based on superfamily designations, which are pictured adjacent to the phylogenies. Only four of the five families are depicted due to the inability to obtain an image of Bathydoridoidea. Circles represent posterior probabilities (top) and non-parametric bootstrap support values (bottom). Closed circles indicate high Bayesian and ML support (pp ≥ 0.95; bs ≥ 75). Red circles indicate moderate support values (pp: 0.95–0.90; bs: 75–70). Open circles indicate no support (pp < 0.90; bs < 70). Relationships that were not recovered by ML analysis are represented by dashes. (a) MAFFT alignment. (b) Gblocks alignment. (c) Curated alignment.

The Doridina suborder was recovered as monophyletic with high support in all analyses estimates (CA: pp = 1.00, bs = 91; GA: pp = 1.00, bs = 86; MA: pp = 1.00, bs = 86; figures 4–6). GA was the only analysis that recovered the Aegiridae sister to all other members of the Doridina, but there was no support. Both MA and CA recovered Prodoris sister to the rest of the Doridina.

Figure 4.

Phylogenetic estimate from MAFFT Alignment. Topology represents Bayesian estimate, with posterior probabilities (pp) and non-parametric bootstrap (bs) support values depicted above and below each branch, respectively. Relationships that were not recovered by ML analysis are represented by dashes. Branches are coloured based on family designations and represented on the exterior of the phylogeny.

Figure 6.

Phylogenetic estimate from Curated alignment. Topology represents Bayesian estimate, with posterior probabilities (pp) and non-parametric bootstrap (bs) support values depicted above and below each branch, respectively. Relationships that were not recovered by ML analysis are represented by dashes. Branches are coloured based on family designations and represented on the exterior of the phylogeny.

Phylogenetic estimate from Gblocks alignment. Topology represents Bayesian estimate, with posterior probabilities (pp) and non-parametric bootstrap (bs) support values depicted above and below each branch, respectively. Relationships that were not recovered by ML analysis are represented by dashes. Branches are coloured based on family designations and represented on the exterior of the phylogeny.

Onchidoridoidea and Phyllidioidea were the only two superfamilies recovered as monophyletic in any of the three analyses (figure 3a–c). Onchidoridoidea was recovered as monophyletic only in the CA (pp = 0.99, bs = 46). Phyllidioidea was monophyletic with little to no support in CA (pp = 0.94; bs = 47) and MA (pp = 0.75, bs = 34) analyses. Most families, however, were recovered as monophyletic with high pp and bs support across all three alignments analysed. The Phyllidiidae and Calycidorididae were the only two families that were not supported in all three analyses. In addition, Chromodorididae, Dendrodorididae, Dorididae and Polyceridae were consistently recovered as not monophyletic.

The Doriopsilla clade was recovered more closely related (GA: pp = 1.00, bs = 50; MA: pp = 1.00, bs = 65) to a monophyletic Phyllidiidae than to the Dendrodoris clade. This strongly suggests the Dendrodorididae is not monophyletic. Interestingly, Cadlinidae was recovered as sister to the Phyllidioidea in MA (pp = 0.67, bs = 13) and CA (pp = 0.64), but in our GA estimate Cadlinidae was sister (pp = 1.00, bs = 17) also to a clade formed by Phyllidiidae and Doriopsilla. Aphelodoris was recovered nested within Discodorididae and sister to Atagema cf. osseosa (CA: pp = 0.63, bs = 52; GA: pp = 0.98, bs = 73; MA: pp = 0.95, bs = 67). Furthermore, Vayssierea sp. and Gymnodoris sp. were recovered nested in the Polyceridae and sister to Polycera quadrilineata (CA: pp = 1.00, bs = 64; GA: pp = 0.98, bs = 77; MA: pp = 0.53, bs = 86).

The Polyceroidea, excluding Aegiridae and Hexabranchidae, was consistently recovered closely related to the Chromodorididae (CA: pp = 0.69; GA: pp = 0.99, bs = 36; MA: pp = 0.97, bs = 53). In addition, all three analyses recovered similar topologies that suggest Cadlinella orniatisma to be sister to Hexabranchus sanguineus. MA is the only phylogenetic estimate that recovered significant support for this relationship (pp = 0.96, bs = 82). Onchimira cavifera was recovered nested in the Onchidorididae and sister (CA: pp = 0.98, bs = 82; GA: pp = 0.77, bs = 51; MA: pp = 0.99, bs = 72) to a clade formed by Knoutsodonta and Onchidoris (CA: pp = 1.00, bs = 86; GA: pp = 1.00, bs = 87; MA: pp = 1.00, bs = 91).

There were only two instances where long branches were evident in the Doridina; Vayssierea sp. and a clade formed by Dendrodoris arborescens, D. fumata, D. guttata and D. nigra. Analysis of nDNA (electronic supplementary material, figure S1) and mtDNA data (electronic supplementary material, figure S2) suggests that the long branches may be an artefact of mtDNA data. The nDNA tree depicts no significantly long branches within the Doridina and recovers a monophyletic Dendrodoris clade, which contradicts our mtDNA tree that recovered a polyphyletic Dendrodoris.

Discussion

Resolving the Doridina

Disparity between phylogenetic estimates can confound interpretations and conclusions regarding processes and patterns of lineage diversification. Conflicting estimates are usually consequences of opposing methodologies, which have centred on taxonomic sampling [75-77], molecular markers [12,78], phylogenetic estimates [79-81] and alignment construction [2,3,14,79]. We examined molecular markers, which have been used in varying arrangements, that are most commonly used in nudibranch phylogenetics [21,23,24,49,60-62,82-92].

The large spikes observed in rDNA fragments for our MA and CA are probably a result of highly variable regions or ambiguous sequence calls, which ML is poor at estimating [66], thus overestimating PI towards the tips of their respective tree. Furthermore, PI profiles suggest that information was lost when these regions were excluded from our analyses, and resulted in surprisingly high loss of parsimony informative characters (25%–55%) for rDNA markers. The removal of these variable regions has been shown to negatively affect phylogenetic estimates [9] and may explain why some relationships were not consistently recovered. This supports the position that highly variable loop regions can be vital in resolving some phylogenetic relationships [8]. Unfortunately, we were unable to increase the resolution of the dorid tree by any of our three MSA construction methods. All three PI profiles were fairly consistent in depicting similar curves, and as we have been able to show, these markers are more appropriate for phylogenetic estimates at family- or higher-level classifications.

Another issue we encountered was the potential noise that was incorporated into our estimates by additional taxon sampling. Even though PI profiles suggest these markers were informative at family level, our increased taxonomic sampling may have hindered our ability to recover consistency across analyses. For example, we only recovered a monophyletic Onchidoridoidea in our CA phylogenetic estimate. By contrast, Hallas & Gosliner [62] recovered a mostly resolved monophyletic Onchidoridoidea with significant pp and bs support. Their taxonomic sampling, however, was much more focused and included histone 3 as an additional molecular marker. These contradictory estimates illustrate issues that can result from inappropriate taxonomic sampling [93] and noise incorporated into analyses with inclusion of highly divergent taxa.

In a few instances, however, our expansive taxonomic sampling has illuminated the relationships of some problematic groups, specifically Aphelodoris, Cadlinidae, Cadlinella, Hexabranchidae and Polyceridae, but in relation to the Onchidoridoidea our estimates contradicted previous highly supported hypotheses. Even though we were able to include Onchimira cavifera, it is relevant to state that we were unable to procure other morphologically unique species that may have affected our ability to resolve some family relationships in the Doridina (e.g. Colga, Goslineria, Hoplodoris, Kalinga, Murphydoris, Otinodoris). These findings illustrate that each phylogenetic query has its own set of challenges and optimal sampling strategy [94], and that the focus for each investigation should be carefully calculated.

Doridina relationships

We were not able to confidently investigate patterns of biogeographical, morphological or chemical evolution due to the lack of resolution at the base of our phylogenetic estimates. The present study, however, offers some consistent new insights into Doridina relationships, in part due to our increased taxonomic sampling.

This work reinforces the conclusion from previous studies that traditional phanerobranch and cryptobranch groupings are not monophyletic [36,59,62]. Even though there was only moderate support in our GA, Cadlinidae does appear to be closely related to at least some members of the Porostomata, despite the ambiguous position of some porostomes such as Dendrodoris and Mandelia. In addition, we consistently recovered Gymnodorididae and Okadaiidae nested within the Polyceridae, which together are closely related to the Chromodorididae and Hexabranchidae. Unexpectedly, Cadlinella was recovered sister to the Hexabranchidae in our phylogenetic estimates. Cadlinella was originally included into the Chromodorididae based on morphological similarities [95] and further supported by molecular studies [84,93]. Our broader taxon sampling, however, consistently recovered Cadlinella sister to Hexabranchus and that both of these taxa, together with the Polyceridae, are closely related to the Chromodorididae. This also is supported by the sperm ultrastructure of Cadlinella, which has been shown to be divergent from members of the Chromodorididae [96].

In addition, the yellowish northeastern Pacific species of Doriopsilla were thought to represent a species complex of closely related taxa, but this assumption was not tested by including any species from outside the complex, other than the outgroup taxon D. spauldingi [97]. In our analysis, we included D. miniata from Japan and D. janaina, another eastern Pacific species that has divergent colouration. In all three of our analyses, the ‘species complex’ suggested by Hoover et al. [97] includes members of two separate lineages, rather than a single radiation. This suggests that these species with yellowish colouration and white spots evolved similar colouration convergently rather than by means of radiation from a single common ancestor.

Lastly, the evolution of the gill pocket is further confounded by the recovery of Onchimira cavifera nested within the Onchidorididae. Onchimira cavifera was described as having both cryptobranch and phanerobranch characteristics, and hypothesized as a missing link in the current understanding of gill reduction [37,62,98,99]. Onchimira possess all the characteristics of a phanerobranch: buccal pump, rectangular rachidian tooth and hooked shaped first lateral tooth, but also possesses a fully formed gill pocket and retractable gill, which is typical of cryptobranch dorids. Surprisingly, Onchimira is not closely related to the only other two members of the Onchidorididae that possess similar gill structure to the Cryptobranchia, Calycidoris and Diaphorodoris [62], but instead nested within the Onchidorididae. Based on our estimates, it is unclear how or under what conditions the gill pocket might have evolved or was lost throughout the Doridina because of the lack of resolution at the base of the tree.

Molecular evolution

It is unclear why there are such large inconsistencies between mtDNA and nDNA phylogenies regarding Dendrodoris. To confirm if there was sequencing error, we examined additional specimens of D. fumata and D. nigra to compare to those on GenBank [61,100]. Surprisingly, all sequences collected were identical. Our inclusion of D. atromaculata and D. denisoni, however, suggested that there are possible highly divergent regions among mtDNA sequences. We were unable to compare other species of Dendrodoris from Hirose et al. [100] because they only analysed COI. A complete sampling of Dendrodoris is needed to fully comprehend the discrepancies between mtDNA and nDNA sequences. Furthermore, there appears to be no molecular distinction between Aegires citrinus and A. serenae. Both species are clearly defined by morphological characteristics [101], but both the mtDNA and the nDNA suggest they are in fact the same species. Much like in the Dendrodorididae, it is unclear what molecular mechanisms might have influenced our observations. Further investigations are needed, but are beyond the scope of our data.

Conclusion

We decided to take an approach that used three common methods used in MSA construction for Nudibranchia phylogenetics. As expected, our findings suggested that MSA methodology affected phylogenetic estimates of the Doridina, especially regarding how we decided to align highly variable rDNA regions. We were able to show that the most commonly sequenced molecular markers for the Nudibranchia lacked the robustness to resolve the base of the dorid tree, and manipulation of highly variable regions affected our ability to recover consistent phylogenetic estimates. This effect, however, is most probably dependent upon the size and scope of the phylogenetic query and amount of missing data. These markers are better suited for higher-level classifications as suggested by our PI profiles. Even though the base of the Doridina was unresolved, family-level classifications were mostly supported across our three analyses, and families that were recovered as non-monophyletic were consistent between alignments. Our analyses suggest that the exclusion of variable regions may have weakened our ability to resolve the base of the Doridina, but previous studies that used much larger datasets have benefited from removing these regions (e.g. [109,110]).

Even though the focus of the present study was to understand MSA construction, our estimates of the Doridina also give a frame of reference for allowing more intensive queries into specific family evolutions. For example, the evolution of caryophyllidia in the Discodorididae, molecular evolution in Aegiridae and Dendrodorididae, or the relationships pertaining to the Polyceridae, Okakaiidae, Gymnodorididae, Hexabranchidae and Chromodorididae, which are some of the most morphologically unique and chemically distinct families.

Nudibranch studies unquestionably suffer from a lack of abundant and diverse molecular markers. Studies have argued that increasing molecular markers could resolve problematic relationships [94,102], but an increase in molecular markers does not resolve issues regarding homology and variable region alignments. Automated filtering protocols allow for MSAs to be easily replicated and eliminate the uncertainty of manual curation of alignments; however, these methods are not without error. In addition, there has been little consensus on the soundest method of increasing signal to resolve phylogenies [6]. Genomic tools, which only have recently been used to investigate nudibranch [103] and larger opisthobranch phylogenetics [104,105], have potential of resolving dorid relationships. However, genomic applications also suffer from alignment and homology issues [106,107]. Phylogenetic resolution of the Doridina can greatly benefit from a genomic approach, but it is important to emphasize the critical role MSAs and homology have on phylogenetic studies. Owing to the varying size and scope of molecular and taxonomic sampling, we strongly recommend the exploration of multiple MSA construction methods that can aid in the selection of an approach that best suits the data.