Placozoans are a phylum of nonbilaterian marine animals currently represented by a single described species, Trichoplax adhaerens, Schulze 1883. Placozoans arguably show the simplest animal morphology, which is identical among isolates collected worldwide, despite an apparently sizeable genetic diversity within the phylum. Here, we use a comparative genomics approach for a deeper appreciation of the structure and causes of the deeply diverging lineages in the Placozoa. We generated a high-quality draft genome of the genetic lineage H13 isolated from Hong Kong and compared it to the distantly related T. adhaerens. We uncovered substantial structural differences between the two genomes that point to a deep genomic separation and provide support that adaptation by gene duplication is likely a crucial mechanism in placozoan speciation. We further provide genetic evidence for reproductively isolated species and suggest a genus-level difference of H13 to T. adhaerens, justifying the designation of H13 as a new species, Hoilungia hongkongensis nov. gen. , nov. spec. , now the second described placozoan species and the first in a new genus. Our multilevel comparative genomics approach is, therefore, likely to prove valuable for species distinctions in other cryptic microscopic animal groups that lack diagnostic morphological characters, such as some nematodes, copepods, rotifers, or mites.
Placozoans Grell, 1971, are small, benthic marine animals found worldwide in various habitats [1–6]. To date, only a single species has been described, Trichoplax adhaerens Schulze 1883. Animals are flat and have a typically disc-like morphology but have the capacity to change shape [7–9]. The lack of symmetry axes, neurons, and defined muscle cells, and the presence of only six morphologically distinguishable somatic cell types ([9,10]; Fig 1, S1 Fig), makes the Placozoa morphologically the most simply organized animals. The prominent placozoan modes of reproduction are asexual, i. e. , binary fission and budding [8,9, 11–13] that produce genetically identical clones. Sexual reproduction has rarely been observed under laboratory condition [14–19], but both oocytes and sperm cells have been reported [14,17,19], and fertilization, likely coupled with genetic exchange, was confirmed based on structural similarities of the placozoan eggshell with the fertilization membrane of other animal groups [16]. No sexually reproducing individual has ever been reported from the wild. However, the occurrence and success of sexual reproduction in the field have been demonstrated by DNA sequence analyses, since nuclear-encoded marker genes have revealed the occurrence of allele sharing and mixing of heterozygous alleles in a natural placozoan population isolated from a Caribbean habitat [20]. These molecular signatures for genetic exchange prove that sexual reproduction does occur and that the life cycle is completed in the natural environment. However, all efforts to follow the placozoan embryonic development in the laboratory have failed to date. All embryos died at an early stage during development, never reaching beyond the 128-cell stage [19]. The fragmentation of the nucleus in the zygote [21] was previously suggested as the reason for the termination of development, although this has been questioned [19]. This ambiguity and scarcity of information has, therefore, left us with a large knowledge gap regarding the life history of the Placozoa and has resulted in speculations of the existence of a missing life stage (compare [22]). The genome of the diploid T. adhaerens was sequenced previously [22], revealing that this morphologically very simple animal harbors a rich repertoire of gene families [22]. These families are known from bilaterian animals and are typically associated with a considerable cell type diversity, a complex body plan, developmental processes, and behavioral responses to external stimuli [10,23–31]. Additionally, single-gene molecular phylogenetics have identified a sizeable cryptic diversity within placozoans collected worldwide; but while their gross morphology is highly plastic, morphologically, all isolates fit the description of T. adhaerens [32] (Fig 1). The high intraspecific shape variability, coupled with an ultraconserved internal structure (Fig 1, S1 Fig), does not allow the establishment of reliable diagnostic morphological characters in the Placozoa, hindering attempts to characterize their diversity. While these single-marker studies provided clear indications that additional species may be uncovered in the Placozoa, two fundamental questions remain: how different are placozoans at the nuclear genome level, and what can we learn from comparative genomics about the evolution and diversity of placozoans? To address these questions, we generated a high-quality draft genome of a placozoan lineage that is genetically distantly related to T. adhaerens [3,5] and performed a multilevel comparison, including genome synteny, gene clustering, gene ontology (GO) term enrichment, allele sharing, and cross-phylum comparative distance analyses. This approach, together with the morphological characterization of the lineage H13, allowed us to assign a taxonomic status to morphologically cryptic taxa and led to the establishment of the second placozoan species in a new placozoan genus.
Based on mitochondrial 16S ribosomal DNA (rDNA) analyses, the genetic lineage H13 is among the most distantly related haplotype to T. adhaerens (lineage H1) [5], whose nuclear genome has been sequenced previously [22]. We hypothesized that the substantial 16S rDNA divergence might also be reflected on the whole-genome scale and, therefore, targeted H13 for nuclear genome sequencing. To assemble the genome of H13—a new species described here, called H. hongkongensis nov. gen. , nov. spec. (Fig 1, S1 Fig; see species description in Material and methods; Tables 1 and 2) —we generated 24 Gb of paired-end reads and 320 Mb of Moleculo (Illumina Artificial Long Synthetic) reads. Our final, highly complete 87-megabase assembly contained 669 high-quality and contamination-filtered contigs with an N50 of 407 kb (S1 Table; S2–S4 Figs), 7 megabases smaller than the T. adhaerens contig assembly. The overall calculated genome heterozygosity (based on single-nucleotide polymorphism [SNP] counts, see S2 Table) was 1. 6%, which is moderate for a marine animal but about average when compared to arthropods and high in comparison to terrestrial chordates [33]. This value cannot be compared to T. adhaerens because of the low genome coverage of the latter, which does not allow haplotype phasing. We annotated the genome with a combination of 15. 3 Gb of RNA-Seq and ab initio methods to yield 12,010 genes (S1 Table, S1 & S2 Data). A high percentage of raw reads mapped back to the genome (S3 Table), and between 90. 8%–95. 3% of the 978 genes in the BUSCO v3 Metazoa dataset were identified in the transcriptome and the ab initio gene models, respectively (S4 Table). Together, this suggests an almost complete assembly and annotation, in which 96. 5% of the genes in the H. hongkongensis genome were expressed in what are commonly considered adult animals. In our gene set, H. hongkongensis had 490 more genes than the 11,520 genes reported in the original T. adhaerens annotation from 2008 [22]. We reannotated T. adhaerens with AUGUSTUS and found an additional 1,001 proteins and also managed to complete formerly partial proteins (for T. adhaerens Blast2GO protein annotations see S3 Data). This approach added 4. 4 Mb of exons to the T. adhaerens annotation, an increase of 28% of exonic base pairs to the original annotation. The new T. adhaerens annotation now has 511 more genes than H. hongkongensis, which accounts for some portion of the size difference between the two genomes. Moleculo reads also enabled us to assemble very large reference contigs, the longest being over 2 Mb. We compared the organization of genes in H. hongkongensis to the 10 longest scaffolds in the T. adhaerens genome (size range 2. 4–13. 2 Mb; accounting for 66% of the T. adhaerens assembly). We found 144 contigs >100 kb from H. hongkongensis that aligned to these 10 scaffolds, accounting for 69% of the H. hongkongensis assembly (Fig 2A). Mean gene collinearity (i. e. , the same genes in the same direction) in this reduced genome representation was in the range of 69. 5% to 78. 8% (mean 73. 6% ± 5. 5%; see S5 Table). The mean number of genes per syntenic block was 33. 8 (±25. 2) in the reduced set and 33. 9 (±24. 7) when comparing full genomes (S5 Fig), which indicates that the reduced set is representative for both complete genomes. Although much of the gene order is conserved between the two species, we counted 2,101 genes (out of the 8,260 genes in the 10 scaffolds) that were inverted or translocated within the same scaffold relative to the order in the T. adhaerens scaffolds. These numbers seem low when compared to the fast-evolving bilaterian genus Drosophila [34,35] or Caenorhabditis [36], but they are in the range of rearrangements found between mouse and human [37]. Comparison to Bilateria, however, might be misleading (see also results on genetic distances below), and genome rearrangement events might be more favored in some bilaterian taxa because of inherent genomic traits such as transposon-induced rearrangement hotspots [38]. Nonetheless, the high percentage of rearrangements between T. adhaerens and H. hongkongensis is clear evidence for a deep genetic separation of both lineages. To estimate how divergent the two placozoan genomes are at the sequence level, we calculated genetic distances for 6,554 one-to-one orthologs. Between H. hongkongensis and T. adhaerens, genetic distances ranged from 0. 9% to 80. 1% (mean 28. 3% ± 12. 9%) for proteins and 7. 4% to 80. 7% (mean 28. 5% ± 9. 9%) for coding sequences (CDSs), respectively (Fig 2B). To assess if certain genes are under positive (diversifying) selection, indicative of functional evolution, we calculated the ratio of nonsynonymous to synonymous nucleotide substitutions (dN/dS ratio [39]) for each H. hongkongensis and T. adhaerens one-to-one ortholog pair. Results show that most orthologs (97%) are under strong purifying selection (dN/dS < 0. 5). One might hypothesize that strong purifying selection pressure is the reason for the phenotypic stasis we see in modern placozoans. However, more placozoan genomes in the phylum are clearly needed to test this hypothesis. Despite this strong tendency toward purifying selection, a high proportion of orthologs (46%) showed larger protein distance than CDS distance and, therefore, an accumulation of double or triple mutations in already mutated codons, which led to amino acid substitutions (S6 Fig). Only 3 of the 6,554 one-to-one orthologs had dN/dS ratios slightly >1, indicating positive selection (S7 Data; see S6 Fig for an estimate of mutation saturation in codons). One of these seems placozoan specific, since it could not be annotated because of missing UniProt BLAST hits and InterPro domains, respectively. For the second, GO annotation and InterPro IDs indicate a role in telomere maintenance. The third positively selected gene (CYP11A1) is putatively a cholesterol side-chain cleavage enzyme acting in the mitochondrion. The roughly 4x coverage of the genome with long Moleculo reads (N50 of 5. 4 kb) allowed the assembly of large haplocontigs (i. e. , contigs representing both haplotypes of the genome). This phasing information for large parts of the genome facilitated the isolation of 2,870 one-to-one orthologs with both full-length alleles after a highly stringent filtering procedure. Only by using the phasing information we were able to show that many orthologs with high allelic variation in H. hongkongensis were also profoundly different between the species (S7 Fig). This indicates that genetic sequence adaptation already takes place at the population level and is further magnified between species in the same genes. The Markov cluster (MCL) analysis identified 6,644 true one-to-one orthologs (for an overview of ortholog categories, see Material and methods and [40]) for both placozoan species (55% of all proteins in H. hongkongensis and 53% in T. adhaerens, respectively) (S8 Fig). A fraction of 465 (3. 8%) H. hongkongensis and 1,036 (8. 3%) T. adhaerens proteins, respectively, did not have reciprocal BLAST hits. The difference in the non-BLAST hits almost perfectly matches the differences in total gene numbers, which is probably an indication that genes without a homolog in H. hongkongensis account at least partially for the slightly higher gene number in T. adhaerens. A high proportion of proteins had BLAST hits to the UniProt database, and only 15. 4% (1,859) and 19. 0% (2,384) of H. hongkongensis and T. adhaerens proteins, respectively, did not have BLAST hits to metazoans included in UniProt. Placozoan-specific duplications constitute a significant proportion of both proteomes, with 3,943 (32. 8%) co-orthologs in H. hongkongensis and 3,484 (27. 8%) in T. adhaerens. The enrichment analyses for the proteins in each non-BLAST-hit bin identified unique GO terms in all three GO categories among the first five most significantly enriched GO terms (S4 & S5 Data). The same applies to one-to-many and many-to-one co-orthologs in both species. The enrichment analyses further indicate that both placozoan species have multiple co-orthologs associated with G-protein-coupled receptor (GPCR) signaling. A rich repertoire of GPCRs has been identified in T. adhaerens [22], but here, we were able to identify independent GPCR duplications in H. hongkongensis and T. adhaerens, respectively (S6 Data). Furthermore, we identified multiple enriched GO terms related to synaptic activity in all co-ortholog categories (S5 Data) and both placozoan species. This points to a plethora of independent duplication events in gene families related to sensory capacities. Despite lacking neurons (based on traditional morphological classifications), T. adhaerens has previously been shown to stain positive for FMRFamide [10,41] and recently even to change behavior when exposed to physiologically relevant levels of neuropeptides [31]. Based on the identification of vast and independent gene family expansions in both placozoans, we propose that adaptation in the Placozoa, ultimately leading to speciation, is coupled with independent gene duplications as suggested, for example, for bacteria, yeast, plants, and other animals (compare [42–45]). H. hongkongensis was isolated from a stream running through a mangrove with rapid drops in salinity and temperature, especially during heavy rainfall in the summer. We hypothesize that the presence of multiple divergent copies of genes involved in various processes, such as behavior and metabolism (compare [42,43]), in addition to a situation-dependent expressional fine-tuning of these copies was necessary for adaptation to this habitat and would facilitate speciation. We furthermore propose that the presence of multiple copies of genes and their expression does not affect the phenotype but instead provides a genetic toolkit for gradual physiological responses to (changes in) the environment. All internal Linnaean ranks within the Placozoa are, as yet, undefined [5]. Despite efforts to identify them, reliable diagnostic morphological characters, commonly used for defining animal species, are lacking in the Placozoa [46]. Thus, all present taxonomic definitions in the phylum must solely rely on diagnostic molecular characters. In other taxonomic groups (e. g. , bacteria and archaea [47], protists [48,49], and fungi [50]), purely sequence-based approaches and working models for the distinction of taxa have been proposed and are generally well established and widely accepted [51]. In animals, such methods (which may be based on distances, on trees, or on allele sharing; [52]) are currently under development and have been used in rare cases to identify and describe cryptic species [53]. In a first step to converting the identified genomic differences into a taxonomically meaningful system, we studied reproductive isolation by addressing allele sharing within placozoan isolates from different localities. To identify reproductive isolation, a conspecificity matrix (CM) was generated [54]. The CM was based on three nuclear genes encoding ribosomal proteins and clearly identified reproductive isolation between placozoan clades (Fig 3). This approach extends a previous study that has uncovered sexual reproduction only within one placozoan haplotype (H8) [20] and provides clear evidence that the previously established placozoan clades (based on 16S genotyping) are reproductively isolated biological species. We have shown that biological species exist in the Placozoa. Previous studies have furthermore provided first indications for the existence of deeper differences between placozoan lineages [1,3], with as-yet-unknown correspondence to, for example, the Linnaean ranks of genus, family, order, and class. However, these observed deeper divergences were based on single marker genes only, and no diagnostic morphological traits could be identified to establish a firm, higher-level, systematic framework in the Placozoa. To further estimate the level of taxonomic relatedness between T. adhaerens and the new placozoan species H. hongkongensis (strain H13), and in an attempt to initiate a higher-level taxonomic system for the Placozoa, we performed cross-phylum multimarker sequence divergence analyses. To do so, we compared the variation between the two placozoans to variation within the other three nonbilaterian phyla, Cnidaria, Ctenophora, and Porifera (compare [1]), as well as the bilaterian phylum Chordata. Marker sets included a nuclear protein set of 212 concatenated proteins (dataset 1, a taxon-extended matrix from [56]; S7–S9 Tables; see Fig 4) as well as 5 selected genes with different substitution rates (S9–S14 Figs), all commonly used for DNA barcoding and molecular systematics. Across individual markers, it appears that the phylogenetic ranks are most robust in the Cnidaria, in which the partitioning of molecular variation matches the established taxonomy, in that Linnaean ranks consistently correspond to the greater distance between groups (Fig 4; S9–S14 Figs). The same is true for the Chordata, which was included in our distance calculations for the 212 nuclear protein set as an example of a bilaterian phylum with a high taxonomic coverage (many genomes are available for this group). However, distances in chordates are, in general, much lower when compared to the overall more similar nonbilaterian phyla. This indicates that (i) genetic distances and corresponding Linnaean rank assignments in Chordata cannot be compared to nonbilaterian lineages and (ii) that comparisons among nonbilaterians are better suited to guide taxonomic ranking of the two placozoan species. We consequently used genetic distances in the Cnidaria as an approximation and comparative guideline for the higher systematic categorization of the new placozoan species. Genetic distances between H. hongkongensis and T. adhaerens were higher than those for the Cnidaria in five of the six marker sets at the generic level but lower at the family level for all markers (S14 Fig, S10 Table), which, cautiously interpreted, supports genus-level genetic differences between the two placozoans. A clear split of the Placozoa in the molecular groups “A” and “B” was previously shown by the rearrangement pattern of mitochondrial genomes [61] and compensatory base changes in the internal transcribed spacer 2 (ITS2) [55]. The conspecificity analysis, the high amount of genomic rearrangement, and the large-scale independent gene duplication history, as well as the genetic distances in six independent datasets, strongly support this split (Fig 3). Since clades were identified as the primary taxonomic units—i. e. , biological species—these two previously identified higher-level placozoan “groups” consequently represent at least the genus level in the Linnaean hierarchical system. We therefore establish the new genus Hoilungia for the former group “A” (clades III–VII), which is, so far, the single sister genus to Trichoplax (former group “B”; clades I and II). Future research efforts focusing on genome sequencing of additional placozoan clades/species will likely help to establish a broader and more detailed systematic framework for the Placozoa and provide further insights into the mechanisms and driving forces of speciation in this enigmatic marine phylum. Recent discussions about the phylogenetic position of placozoans have largely been based on the T. adhaerens genome. A better sampling of placozoan genomic diversity is, however, needed [62] to address their placement in the metazoan tree of life. In this context, it is important to first assess if adding another placozoan genus would break up the long placozoan branch. The inclusion of a single representative of a clade with a very long terminal branch, or fast-evolving taxa that can have random amino acid sequence similarities, may result in erroneous groupings in a phylogeny (so-called “long-branch attraction artefacts”) [63,64]. To address these questions, we generated a highly (taxa) condensed version of the full protein matrix from Cannon and colleagues [56] (termed dataset 2; with less than 11% missing characters and 194 genes). We additionally created a Dayhoff 6-state recoded matrix [65,66] of this second set to reduce amino acid compositional heterogeneity, which is also known to be a source of phylogenetic error [67,68]. Phylogenetic analyses were performed on these two matrices (protein and Dayhoff-6 recoded), using the site-heterogeneous CAT-GTR model in PhyloBayes-MPI [69] and using the site-homogenous GTR model both in Phylobayes-MPI and RAxML (RAxML, protein only) [70], as well as the LG model in RAxML (protein only). The resulting trees (S15–S20 Figs) of the highly dense gene matrix (S21 Fig) suggest a sister group relationship of the Placozoa to a Cnidaria + Bilateria clade with both CAT-GTR (Protein, Dayhoff-6 recoded, S15–S17 Figs) and GTR models (Protein, S18 Fig) in PhyloBayes, or these relationships are unresolved (RAxML, protein, both GTR, S19 Fig, and LG, S20 Fig). This is in agreement with some previous findings [56,64,71–74] and with recent studies using a large gene set and intense quality controls [64] as well as improved modeling of compositional heterogeneity [68]. In addition, the sister group relationship of the Placozoa to the Cnidaria + Bilateria clade is corroborated by independent data—namely, the analysis of metazoan genome gene content [73,75,76].
Phylum: Placozoa, Grell 1971 [77] Type Family: Trichoplacidae, Bütschli and Hatschek 1905 [78], synonymized with “Trichoplaciden” (in German original), Haeckel 1896 [79]. Diagnosis: We assign all currently known 19 placozoan genetic lineages (16S haplotypes H1-H19; [5]) to the Trichoplacidae. The description of T. adhaerens Schulze 1883 applies to all. Type Genus: Hoilungia, nov. gen. , Eitel, Schierwater, and Wörheide Hoilungia is the second genus of the family Trichoplacidae. Etymology: Hoilungia, pseudo-Latinized from “Hoi Lung, ” Cantonese, meaning “sea dragon, ” which is based on the shape-shifting dragon king in Chinese mythology. Diagnosis: Gross and fine morphology appear similar among all placozoans studied to date. We therefore use molecular diagnostics to define Linnaean ranks. Among all tested markers, the mitochondrial large ribosomal subunit 16S rDNA appears to be the most variable among placozoans and other nonbilaterian phyla, and the mean pairwise distance is closest to that calculated for the nuclear dataset in most cases (S14 Fig). This marker also best mirrored classical taxonomy in the Porifera and Cnidaria (S11 Fig; in Ctenophora, 16S rDNA is highly derived and hard to identify [80]). According to these data, molecular diagnostics based on differences in the 16S rDNA appear to be suitable for current and future designation of species in the Placozoa, which is in agreement with previous results [3]. Diagnostics are here, therefore, defined by nucleotide substitutions in the 16S rDNA. Full-length 16S rDNA sequences of T. adhaerens and H. hongkongensis (clonal strain “M2RS3-2”), as well as for the undescribed Placozoa sp. H4 and sp. H8, were aligned with MAFFT v7. 273 [81] using the GINSI option and otherwise default settings. Ambiguously aligned 5′ and 3′ sequence ends were removed. To this alignment, we added all currently available placozoan 16S haplotype sequences [5] using MAFFT [added option: —add]. The final alignment contained all 19 placozoan haplotypes and had a length of 2,551 nucleotides (including gaps). The region for identification of diagnostic nucleotides was restricted to a part of the 16S alignment that was previously shown to be suitable and sufficient for molecular haplotype discrimination [1,3, 5]. We furthermore restricted the identification of diagnostics to stem regions of this rDNA to omit uncertainties in future taxonomic assignment due to ambiguously aligned loop regions. To identify molecular diagnostics for the genus Hoilungia, we screened for molecular synapomorphies (nucleotide exchanges) within the placozoan 16S group “A” (clades III–VII; [5,61]) versus group “B” (clades I and II). Molecular diagnostics for Hoilungia and Trichoplax are summarized in Table 1. Type species: H. hongkongensis, nov. spec. , Eitel, Schierwater, and Wörheide. Diagnosis: To identify molecular species diagnostics, we determined unique substitutions (based on the alignment used for genus diagnostics before) for H. hongkongensis (clade V) in comparison to the other Hoilungia clades (III, IV, VI, and VII). Molecular diagnostics for H. hongkongensis are summarized in Table 2. Type locality: A single specimen of H. hongkongensis (clonal strain “M2RS3-2”) was isolated in the Ho Chung River close to a small mangrove at Heung Chung village, Hong Kong (22. 352728N 114. 251733E), on June 6,2012. Type specimen: One specimen of H. hongkongensis (clonal strain “M2RS3-2”) has been mounted and deposited at the Bayerische Staatssammlung für Paläontologie und Geologie in München, Germany, under voucher number SNSB-BSPG. GW30216. Clonal individuals have been stored in ethanol as paratypes under voucher number SNSB-BSPG. GW30217 in addition to a DNA extraction under voucher number SNSB-BSPG. GW30218. Etymology: hongkongensis, from “Hong Kong, ” and “-ensis, ” Latin, suffix referring to place of origin, as specimens are at present endemic to Hong Kong. The full name “Hoilungia hongkongensis” thus means “Hong Kong sea dragon. ” Two strains were used for this project: The “M2RS3-2” strain was used for the DNA sequencing (the “DNA strain”) and the “M153E-2” strain (the “RNA strain”) for the transcriptome. Both strains descend from a single placozoan individual each, which was isolated from mangroves/mangrove associates at two different sites in Hong Kong (SAR, China). The DNA strain was isolated from a dead mussel shell collected in the Ho Chung River close to a small mangrove at Heung Chung village (22. 352728N 114. 251733E) on June 6,2012. The habitat undergoes daily changes in salinity, and on the day of collection, the salinity was 20 psu. The RNA strain was isolated from collection traps (for details on slide sampling, see [82]) connected to mangrove associates (Hibiscus sp.) and high shore mangrove (Excoecaria sp.) trees at Tai Tam Tuk (22. 244708N 114. 221978E) on March 30,2012. Both clonal cultures were cultured in 14 cm glass Petri dishes as described [19], with a pure Pyrenomonas helgolandii algae culture (strain ID 28. 87, Culture Collection of Algae, Georg-August-Universität Göttingen). The two different strains were used for DNA and RNA sequencing, respectively, to identify polymorphisms in these strains living in the same habitat but at two hydrogeographically distinct sampling sites (northeast versus southeast Hong Kong). Animals were transferred in 20% BSA in artificial seawater, high-pressure frozen in a Wohlwend HPF Compact 02, and stored in liquid nitrogen. Samples were processed from −90 °C to room temperature for Epon embedding in a Leica AFS unit as follows: they were fixed and contrasted in 0. 1% tannic acid in acetone for 24 h and washed 4 times for 15 min in acetone; samples were then incubated in 2% Osmium tetroxide in acetone while the temperature was increased stepwise to −40 °C within the next 23 h; samples were then washed and progressively infiltrated in Epon: acetone mixes (1: 2,2: 1) and pure Epon while temperature was further raised from −40 °C to room temperature over 6 h. They were then polymerized in Epon. Seventy-nm ultrathin sections were cut on a Leica Ultracut and picked up on a copper slot grid 2 × 1 mm coated with a polystyrene film. Sections were poststained with uranyl acetate 2% in distilled water for 10 min, rinsed several times with distilled water followed by Reynolds lead citrate in distilled water for 10 min, and rinsed several times with distilled water. Micrographs were taken with a Transmission Electron Microscope Philips CM100 at an acceleration voltage of 80 kV with a TVIPS TemCam-F416 digital camera. A" lavalamp" kmer/GC plot was generated (S2 Fig) to yield a high-resolution plot of read counts per %GC and 31 bp kmer coverage using the Jellyfish kmer counter and a set of custom Python scripts (kmersorter. py and fastqdumps2histo. py; for details on the procedure, see https: //github. com/wrf/lavaLampPlot). In contrast to the conceptually similar approach Blobtools [108], we used raw reads instead of contigs to yield a high-resolution plot of read counts per %GC and 31 mer coverage. The plot identified two read clouds with high counts at a kmer coverage of 80–140x (heterozygous “read cloud”) and 160–260x (homozygous “read cloud”), respectively. Additional “read clouds” at 270–320x and 380–410x coverage mark repetitive sequence stretches. Another “read cloud” was found at a low coverage of 20–50x. Reads within this cloud and their pairs were extracted with kmersorter. py [added options: -s 0. 16 -b 50 -w 0. 40 -T -k 31] and fastqdumps2histo. py. Bowtie2 v2. 2. 5 [109] [added options: -q—no-sq] was used to map the 580,092 extracted reads to the 19 previously identified bacterial contigs (see section “Contamination screening”). More than 86% of these reads mapped to the bacterial contigs, confirming the bacterial origin of the reads within the low-coverage “read cloud. ” Read counts identified a relatively high abundance of bacterial cells, and the GC content was similar to the host genome. To estimate the per-base genome coverage, paired-end reads were mapped to the softmasked reference assembly with Bowtie2 v2. 2. 5 [added options: -q—no-unal—no-sq) and sorted with SAMtools v1. 3. 1 [110]. The bam file was used to create a bedgraph file in BEDtools v2. 25. 0 [111] by invoking the genomecov operation [added options: -ibam stdin -bga]. A custom Python script (bedgraph2histo. py) [added options: -m 2000] was used to create a coverage histogram table. Since 81. 4% of the genome falls within the second peak (165–332x coverage with a maximum at 248x), most of the genome was merged in the reference assembly (S3 Fig). To identify collinearity between the two placozoan species, all H. hongkongensis contigs >100 kb were aligned to the longest 10 T. adhaerens scaffolds (accounting for 70. 3 Mb or 66. 5% of the genome assembly; including 5. 7-Mb gaps) with default settings. For generating the alignments, LASTZ v1. 02. 00 [115] (implemented as a plugin in Geneious) was used. Of the 222 H. hongkongensis contigs >100 kb, a total of 144 (accounting for 60. 6 Mb or 69. 4% of the genome assembly) aligned to the 10 longest T. adhaerens scaffolds. Aligned H. hongkongensis contigs were extracted from the assembly, sorted, and occasionally reverse complemented to be oriented according to the T. adhaerens scaffolds. Gene annotations (GFF) of contigs as well as protein sequences were extracted for the target scaffolds/contigs sets of both species. A MCScanX run [116] [added option: -a] was performed for each target set, using the extracted T. adhaerens and H. hongkongensis GFFs together with the reciprocal best 5 BLASTP hits [added options: -evalue 1e-10 -max_target_seqs 5 -outfmt 6] between and among proteins of both placozoans. Dual synteny line plots of the resulting collinearity files were visualized in VGSC v1. 1 [117] [added options: -tp DualSynteny] and combined to Fig 2A. In addition, bar plots were generated for the 10 T. adhaerens scaffolds and the matching 144 H. hongkongensis contigs in VGSC [added option: -tp Bar]. Bar plots were mapped onto the DualSyntheny plots to show collinearity within each set and macrosynteny between both genomes. The percentage of collinearity between the T. adhaerens scaffolds and H. hongkongensis contigs was calculated in MCScanX, and results for the 10 scaffolds are given in S5 Table. The mean collinearity was calculated as the sum of the individual collinearities for the 10 T. adhaerens scaffolds multiplied by a size-correction faction for each scaffold (i. e. , percent coverage of the evaluated 70. 4 Mb of the T. adhaerens genome). Syntenic block sizes and the number of blocks were calculated using the custom Python script microsynteny. py (described in [118]) with skipping no more than 1 gene [added option: -s 1] and otherwise default options. To identify allele sharing or reproductive isolation, 3 genes encoding for ribosomal proteins were amplified via PCR, using degenerate primers designed based on the T. adhaerens genomic sequence, as well as a previously sequenced EST library of lineage H4 [19]. Primer sequences to amplify gDNA (including intronic sequence) for the ribosomal proteins L9 and L32, as well as ribosomal protein P1, were as follows: PCRs were run with an initial denaturation of 3 min at 94 °C; followed by 40 cycles of 30 s of denaturation at 94 °C, 30 s of annealing at 60 °C, and 1. 5 min of elongation at 72 °C; and finished with a final elongation for 3 min at 72 °C. The BIOTAQ system was used (Bioline, London). A list of samples used for amplification is provided as S6 Table. Sequencing was performed by Macrogen (South Korea). Alleles were identified as double peaks in standard sequencing in the case of heterozygous alleles. The phasing of SNPs was inferred from homozygous sequences as well as the sequence of allelic variants in closely related haplotypes, for which phasing information was available because of the long Sanger reads. To check for reproductive isolation and to identify conspecific isolates, haplowebs [123] were generated for each marker as well as a CM [54] for combined markers using the online tool HaplowebMaker (https: //eeg-ebe. github. io/HaplowebMaker/; Spöri & Flot, in prep.). The resulting conspecificity scores were plotted in R using the heatmap3 package [124], sorted according to a UPGMA tree (JC69 model) of the three concatenated genomic sequences (with indels removed). If present, both alleles of an isolate were merged, and the consensus sequence was used to generate the tree. dN/dS ratios—as well as fractions of unchanged codons, synonymous, and nonsynonymous sites—were calculated based on a custom Python script (alignmentdnds. py) using regapped CDS alignments and untrimmed protein alignments (S6 Fig). Codons with any ambiguous bases and gapped sites were ignored. Clustering into homologs and co-orthologs was performed with a custom python script (makehomologs. py) [added options: -s 1 -p 234 -H 200]. The script calls the MCL v12-068 algorithm [127], which uses the output of a local all-versus-all BLASTP search [added options: -evalue 1e-3 -outfmt 6] of all H. hongkongensis and T. adhaerens proteins. To identify enriched GO terms in non-BLAST hits as well as in four co-ortholog categories (one-to-many, many-to-one, many-to-many, and many-to-zero), an enrichment analysis was performed for the three main GO categories (Biological Process, Cellular Component, and Molecular Function) using topGO [128]. Only enriched GO terms with a p-value <0. 05 were kept, based on the classic Fisher test. Ortholog categories (see also [40]) are defined as (1) one-to-one: Only one ortholog is found in each species; (2) one-to-many: One ortholog in this species, but many co-orthologs in the other species. The gene was duplicated in the other species from the ancestral copy after speciation; (3) many-to-one: More than one co-ortholog in this species but only one in the other species. The gene was duplicated in this species from the ancestral copy after speciation; (4) many-to-many: More than one co-ortholog in this and the other species. At least two gene duplications could be found from an ancestral gene in the common ancestor of both species—one duplication in this species, and a second one in the other species; (5) many-to-zero: Many co-orthologs in this species but none in the other. In this case, the gene was duplicated from an ancestral copy in this species after speciation and likely lost in the other species. To estimate molecular differences between H. hongkongensis and T. adhaerens and to bring these into a taxonomic context, we measured genetic distance using an extended data matrix of 212 nuclear proteins set up by Cannon and colleagues [55]. This data matrix was chosen as it includes a comparable number of sites for a diverse taxonomic range and is, therefore, also suitable for phylogenetic analyses. In addition, genetic distances were measured for 5 standard barcoding (“selected”) markers—namely, nuclear ribosomal subunits 18S (S9 Fig) and 28S (S10 Fig), mitochondrial large ribosomal subunit 16S (S11 Fig), and the mitochondrial proteins cytochrome c oxidase subunit 1 (CO1) (S12 Fig) and NADH dehydrogenase subunit 1 (ND1) (S13 Fig). An overview of means for all distances of all six marker sets is provided as S14 Fig. The incorporation of datasets from four individual categories (nuclear protein versus nuclear rDNA versus mitochondrial protein versus mitochondrial rDNA) enabled the comparison among markers with different substitution rates. To assess the effect of adding a second placozoan species on the placement of the Placozoa in the animal tree of life and to estimate branch lengths to the two placozoan species, dataset 1 was further condensed to generate a highly complete protein matrix (dataset 2). This set had only 10. 8% missing characters in 58 taxa, including 32 nonbilaterians and 2 outgroups with an almost complete gene set (194 genes, see also gene density matrix in S21 Fig). It has been demonstrated that the CAT model (specifically CAT-GTR) implemented in PhyloBayes [132] fits phylogenomic amino acid supermatrices containing nonbilaterians best [73,133], and obviously, only best-fitting evolutionary models should be used in probabilistic phylogenetic analyses to reduce systematic errors [133]. However, the computational burden of reaching convergence of analyses using the CAT-GTR model can be prohibitive. It is also well known that phylogenomic datasets frequently suffer from compositional heterogeneity that might negatively influence phylogeny estimation [134–136]. Compositional heterogeneities can be reduced by the so-called Dayhoff recoding [65,137,138], which combines amino acids with similar physicochemical properties into one of six categories. Through this reduction of character space, lineage-specific compositional heterogeneities are lessened—at the cost, however, of losing phylogenetic signal [67]. However, another advantage of Dayhoff recoding is a significant reduction of computation time needed to reach convergence. The protein as well as the Dayhoff 6-state recoded dataset 2 were analyzed with PhyloBayes-MPI v1. 7 [69,132], employing the CAT-GTR model, on the Linux cluster of the Leibniz Rechenzentrum (http: //www. lrz. de) in Garching bei München, running 2 chains (each on 112 CPUs) each until reaching convergence, as estimated by using tracecomp and bpcomp programs of the PhyloBayes package (see PhyloBayes manual for details). Furthermore, to evaluate the effect of using less-fitting site-homogeneous evolutionary models on the phylogenetic relationships of the Placozoa, we conducted a PhyloBayes-MPI analysis as above but with the GTR model (see for example [73], [68]), and also two maximum-likelihood analyses in RAxML: one with the GTR model using RAxML-NG v0. 5. 1b [139] [added options: —model PROTGTR+G+I—bs-trees 100—data-type AA] and one with the LG model using RAxML v8. 2 [70] [added options: -f a -x 670 -m PROTGAMMAILG -p 220 -N 100]. The LG model was used as it was the best-fitting site-homogeneous model in 210 of the 212 gene partitions determined by ProtTest v3. 4 [140]. Phylogenetic trees are shown as S15–S20 Figs.
Placozoans are a phylum of tiny (approximately 1 mm) marine animals that are found worldwide in temperate and tropical waters. They are characterized by morphological simplicity, with only a handful of cell types, no neurons, no tissue organization, and even no axial polarity. Since the original description of Trichoplax adhaerens 135 years ago, no additional accepted species has been established, leaving the Placozoa as the only animal phylum with only a single formally described species. While classical morphological species identification has failed to reveal further species, single-gene DNA sequence analyses have identified a broad and deep genetic diversity within the Placozoa. To address the significance of this deep genetic diversity in this morphologically uniform phylum, and to better understand its consequences for speciation processes, general biology, and species delimitation in the Placozoa, we sequenced the genome of the placozoan isolate “H13, ” a lineage distantly genetically related to T. adhaerens. Our multilevel genomic comparisons with the T. adhaerens genome show considerable differences in the general structure of the genome and the makeup and history of various gene families of biological relevance to habitat adaptation. Based on comparative genomics, we here describe the second placozoan species and show that it belongs to a new genus.
taxonomy
invertebrates
animals
trichoplax adhaerens
placozoa
phylogenetics
data management
phylogenetic analysis
energy-producing organelles
mitochondria
bioenergetics
cellular structures and organelles
research and analysis methods
sequence analysis
computer and information sciences
sequence alignment
bioinformatics
comparative genomics
short reports
evolutionary systematics
biochemistry
eukaryota
cell biology
cnidaria
database and informatics methods
genetics
biology and life sciences
genomics
evolutionary biology
computational biology
organisms
Comparative genomics and the nature of placozoan species