A comprehensive comparison of four species of Onchidiidae provides insights on the morphological and molecular adaptations of invertebrates from shallow seas to wetlands

The Onchidiidae family is ideal for studying the evolution of marine invertebrate species from sea to wetland environments. However, comparative studies of Onchidiidae species are rare. A total of 40 samples were collected from four species (10 specimens per onchidiid), and their histological and molecular differences were systematically evaluated to elucidate the morphological foundations underlying the adaptations of these species. A histological analysis was performed to compare the structures of respiratory organs (gill, lung sac, dorsal skin) among onchidiids, and transcriptome sequencing of four representative onchidiids was performed to investigate the molecular mechanisms associated with their respective habitats. Twenty-six SNP markers of Onchidium reevesii revealed some DNA polymorphisms determining visible traits. Non-muscle myosin heavy chain II (NMHC II) and myosin heavy chain (MyHC), which play essential roles in amphibian developmental processes, were found to be differentially expressed in different onchidiids and tissues. The species with higher terrestrial ability and increased integrated expression of Os-MHC (NMHC II gene) and the MyHC gene, illustrating that the expression levels of these genes were associated with the evolutionary degree. This study provides a comprehensive analysis of the adaptions of a diverse and widespread group of invertebrates, the Onchidiidae. Some onchidiids can breathe well through gills and skin when under seawater, and some can breathe well through lung sacs and skin when in wetlands. A histological comparison of respiratory organs and the relative expression levels of two genes provided insights into the adaptions of onchidiids that allowed their transition from shallow seas to wetlands. This work provides a valuable reference and might encourage further study.

Introduction specific protein as well as the primary protein in muscle [14]. In fact, the myosin heavy chain protein is related to the contraction of muscle and is thus of interest when analyzing muscle adaptations. Furthermore, since the discovery of myosins in non-muscle cells, it has been suggested that these proteins drive morphogenesis for successful development [15,16]. Non-muscle myosin II (NM II) is a suitable candidate for analyzing adaptions because it is present in all tissues and is related to morphological development [16]. Non-muscle myosin heavy chain II (NMHC II) is produced from non-muscle myosin II, which is composed of a pair of heavy chains and two pairs of light chains [17]. Therefore, these two genes are suitable candidates for comparing adaptations among the four species of interest. Interestingly, onchidiids express trait-associated genes in various tissues according to their specific habitats. The comparative analyses performed in this study provide insights into the seawater-to-land transition that has occurred in Onchidiidae.

Ethics statement
This study was carried out in strict accordance with the Guidelines on the Care and Use of Laboratory Animals issued by the Chinese Council on Animals Research and Guidelines of animal Care. The study was approved by the ethical committee of Shanghai Ocean University.

Sample collection
The adult individuals used in this study were collected between May and November. All the animals used in the transcriptome sequencing and qRT-PCR experiments were collected in August and were maintained at 27 ± 1˚C for one week prior to the study. Onchidium reevesii individuals were collected from Shanghai (31˚33 0 N, 121˚48 0 E); Paraoncidium reevesii and Platevindex mortoni individuals were collected from Xiamen, Fujian Province (24˚27 0 N, 1180 4 0 E); and Peronia verruculata individuals were collected from Zhanjiang, Guangdong Province (21˚11 0 N, 110˚24 0 E). All the individuals of the four species (10 specimens per onchidiid) were fed corn flour and reared at room temperature.

Stereomicroscopy, light microscopy and scanning electron microscopy
Three fresh adult specimens of each species were anaesthetized by ether, and their external morphologies were observed under a Olympus SZX16 stereomicroscope. Dorsal and ventral skin samples from the four species of Onchidiidae were dissected into small pieces, fixed in Bouin solution and embedded in paraffin wax [18]. Sections (5~6 μm) were cut on a Leica RM2035 microtome, stained with hematoxylin-eosin and observed under a Nikon Eclipse Ni light microscope.
Ten sections were selected randomly from each specimen for measurements of the skin, epidermis, dermis, stratum compactum and stratum spongiosum at six sites. The data were analyzed using the software package JMP Version 10.0 (SAS Inc., NC, USA).
For scanning electron microscopy (SEM) analysis of the Onchidiidae species, tissues were fixed in a mixture of methanol and glutaraldehyde for one week and then preserved in 75% alcohol. After this procedure, the materials were washed three times in phosphate buffer (pH 7.0) for 15 min each time, cleaned in an ultrasonic water bath for 2~3 min and then dehydrated in a series of increasingly concentrated ethanol solutions (30%, 50%, 70%, 80%, 90%, and 100% ethanol), with 15 min per solution. Finally, the samples were prepared for SEM using critical-point drying, sputter coated with gold using DMX-220 ion-plating equipment and then examined by SEM.

Transcriptome sequencing and sequence analysis
Five active individuals of similar size from each species were selected randomly, and their mantles were removed and individually placed in labeled, RNase-free, 2-ml EP tubes with RNAlater (Qiagen: 76104). All the samples were sent to Genergy Biotechnology (Shanghai) Co., Ltd., for RNA extraction and sequencing. The cDNA library was sequenced by Genergy Biotechnology Company (Shanghai, China) using an Illumina HiSeq TM 2000 (Illumina, Inc., USA). The raw data were processed to remove nonsense sequences and then assembled using the short-reads assembling program Trinity [19,20].
Functional annotation of the transcriptome was performed using Blast2GO software [21][22][23]. For annotation, BLASTX alignment (e value<1e-5) between unigenes and protein databases, such as UniProt (www.uniprot.org) and NCBI NR (NCBI non-redundant nucleotide database, (http://www.ncbi.nlm.nih.gov/), was performed, and the best aligning results were used to annotate the protein functions. Unigene annotation provided functional annotations of unigenes and included protein sequence similarity, GO (Gene ontology, http://www. geneontology.org/GO.slims.shtml) [24] functional classification, and KEGG (Kyoto Encyclopedia of Genes and Genomes, http://www.genome.jp/kegg/) pathway analysis [25]. Representative sequences of 18S were obtained from the GenBank database for phylogenetic analysis. A phylogenetic tree was constructed using the Bayesian method with MrBayes version 3.2.4. The bootstrap test was employed based on 10,000 pseudo-replications to assess the reliability of the phylogenetic tree.

SNP marker development in Onchidium reevesii
We investigated unigenes from the transcriptome of Onchidium reevesii using SAMtools and calling SNPs [26]. Then, potential SNP loci of O. reevesii that differed from those of Peronia verruculata, Paraoncidium reevesii and Platevindex mortoni related to vascularization, muscle development, cuticularization and epidermis formation were selected. Primer pairs were designed using Primer Premier 5.0 (http://www.premierbiosoft.com) and synthetized by Map Biotech (Shanghai China). The primer pairs were then tested in 10 individuals as a preliminary screen. The primers that produced clearly defined bands were further tested for the analysis of polymorphisms in 60 individuals.
Genomic DNA was extracted from the mantle of living O. reevesii adults (collected from Chongming Island, Shanghai, China) using a phenol-chloroform extraction method [27]. Multiplex PCR conditions were standardized to a 20-μl volume containing 4 μl of Primer Mix, 1.6 μl of Mg 2+ , 0.4 μl of dNTP Mix, 10 μl of Ex Taq, 2 μl of DNA, and 2 μl of deionized water. The inactivated multiplex PCR product mix was used for SNaPshot multiple single-base extension reactions. The multiple single-base extension reactions were standardized to 25 μl comprising 5 μl of SNaPshot Multiplex Kit (ABI), 2 μl of multiplex PCR product mix, 1 μl of the extension primer mix, and 2 μl of deionized water. The extension reactions were performed under the following conditions: an initial denaturation at 95˚C for 10 s followed by 25 cycles of denaturation at 95˚C for 10 s, 40 s at 50˚C, and 30 s at 60˚C and a final extension step at 30˚C for 30 s. The products were then tested for polymorphisms using an ABI 3730XL sequencer.
The primary analysis of the ABI 3730XL sequencing data was performed with GeneMapper 4.0 (Applied Biosystems Co., Ltd., USA). The number of alleles (Na), the observed heterozygosity (H O ), the expected heterozygosity (H E ) and the deviations from Hardy-Weinberg equilibrium (HWE) for each locus were calculated with Popgene32 (Version 1.32). Bonferroni correction was used to correct the results. The polymorphism information content was calculated with Cervus 3.0 (http://www.fieldgenetics.com/pages/home.jsp).

Cloning and quantitative analysis of the Os-NMHC and MyHC genes
Samples of the dorsal skin, ventral skin, foot skin, lung sac, ganglion and ventricle were collected from the four species. The samples (three specimens per tissue) were immediately flash frozen in liquid N 2 and maintained at -80˚C until use. Total RNA was extracted from the tissues with RNAiso Plus (TaKaRa, Japan) according to the manufacturer's recommended protocol. Briefly, total RNA was obtained from a mixed extraction of tissues. The A260/280 and A260/230 ratios of the RNA were measured using a NanoDrop 2000 spectrophotometer (NanoDrop Technologies, USA) and equaled 1.90-2.10 and 2.00-2.50, respectively. cDNA was synthesized from the dorsal skin mRNA using an RT reagent kit with gDNA Eraser (TaKaRa, Japan), and the 3' and 5' ends of the cDNA were obtained using the RACE technique (TaKaRa, Japan).
Partial fragments of the Os-NMHC and MyHC genes were obtained from the de novo transcriptomic library. To confirm the fragment sequences, we used specific primers to amplify the partial fragments and re-sequenced the PCR products. The specific primers used for cloning the full-length cDNA of Os-NMHC and MyHC are provided in Table 1. The PCR cycling conditions were as follows: 94˚C for 5 min followed by 30 cycles of 94˚C for 30 s, 58˚C for 30 s, and 72˚C for 1 min. The smart 5'-RACE (5' Full RACE Kit, TaKaRa, Japan) and 3'-RACE kits (3'-Full RACE Core Set Ver. 2.0, Takara, Japan) were used per the manufacturers' instructions. The RACE-PCR product was ligated into pGEM-T Easy vector (Promega, USA) and transformed into competent Escherichia coli DH5-α cells. Using blue-white selection and PCR identification, positive clones were selected and sequenced. Concurrently, cDNAs of other tissues were synthesized for qRT-PCR analysis of Os-NMHC gene expression. In addition, a constitutively expressed gene, 18S, was used as an internal control to verify the fluorescent real-time RT-PCR reactions.
The expression of Os-NMHC and MyHC transcripts in different tissues was studied by fluorescent real-time RT-PCR. Quantitative RT-PCR was performed using the Light Cycler1 480 II instrument (Roche, Swiss) with a QuantiFast1 SYBR1 Green PCR kit (Qiagen, Germany). The reaction conditions were as follows: 94˚C for 5 min followed by 30 cycles of 94˚C for 30 s, 51˚C for 30 s, and 72˚C for 1 min and a final step at 72˚C for 5 min. Data were collected from each qRT-PCR experiment performed in triplicate and expressed as the means ± SE. All the primers used in this process are listed in Table 1.

Comparisons of morphological characteristics among the four species
The stereomicroscopy analysis revealed that the nodular papillae in the dorsal skin were most pronounced in Onchidium reevesii among the four species of Onchidiidae. The most prominent difference between Peronia verruculata and the other three species was the presence of dendritic gills located at the posterior end of the body. P. verruculata has these dendritic gills (Fig 1), which allow it to breathe well when submerged. Moreover, the skin over the gills is thin, and thus, the gills have greater permeability than the other parts of the back skin, although they also have thicker cuticular membranes than the gills of the other three species ( Fig  2). The surfaces of Onchidiidae species are covered with a layer of cuticular membrane, which becomes purple after staining. The epidermis of Onchidiidae species is composed of two to three cell layers, and the epidermis of P. verruculata is highly keratinized. The cells of each layer are abundant and closely arranged in O. reevesii and Platevindex mortoni, whereas in P. verruculata and Paraoncidium reevesii, they are sparsely arranged. We measured the dorsal skin thickness of the four species ( Table 2) and found that P. verruculata had the thickest dorsal skin and that P. reevesii had the thinnest dorsal skin. Granular and mucous glands, which are multicellular glands, were present in the skin of the Onchidiidae species (Fig 2). Onchidium reevesii had the most developed lung sacs, closely followed by Platevindex mortoni and Peronia verruculata, whereas Paraoncidium reevesii had the least-developed lung sacs (Fig 3). The structural differences among the lung sacs of the four species in the Onchidiidae family are described in Table 3. Specifically, Onchidium reevesii has developed reticular septa, secondary septa and third septa, whereas Paraoncidium reevesii only possesses undeveloped reticular septa ( Table 3). The degree of development of the lung sacs in the four species of Onchidiidae decreases in the order Onchidium reevesii, Peronia verruculata, Platevindex mortoni and Paraoncidium reevesii.
In addition to the GO analysis, a KEGG pathway mapping analysis based on the enzyme commission (EC) numbers, which is an alternative approach for categorizing gene functions with an emphasis on biochemical pathways, was performed using the assembled sequences. This analysis revealed that the detected unigenes participated in 129, 136, 138 and 134 pathways in Platevindex mortoni, Paraoncidium reevesii, Onchidium reevesii and Peronia verruculata, respectively. To determine the phylogenetic relationships among the Onchidiidae sequences and their orthologs in other mollusks and amphibians, we constructed a phylogenetic tree using MrBayes version 3.2.4 (Fig 4).

Development of SNP markers of Onchidium reevesii
SNPs are important molecular markers, and the SNPs of Onchidiidae developed in this study were valuable for understanding the species' respiratory traits and amphibious features. The proposed sites in the transcriptome sequences of Onchidium reevesii were searched using SAMtools, and 152,212 SNPs were detected after analysis. Fifty-seven alternative SNP loci of O. reevesii that differed from the SNPs of Peronia verruculata, Paraoncidium reevesii and Platevindex mortoni related to vascularization, muscle development, cuticularization and epidermis formation were selected for further study. Fifteen sequences were selected to design 57 pairs of primers for 57 loci. Forty-two pairs of primers were successfully amplified among these 57 loci, and 26 of these 42 pairs were selected to test for polymorphisms (Table 6). In  total, the observed and expected heterozygosities ranged from 0.2553 to 1.0000 and from 0.0000 to 0.7447, respectively ( Table 6). No genetic linkage was observed among these loci. Fifteen loci with ' Ã ' significantly departed from HWE after Bonferroni correction (P<0.05). Among the 26 SNP loci, three loci (S_Unigene508_c0_seq1_142, S_Unigene685_c0_seq1_3534, and S_Unige-ne508_c0_seq1_283) were related to epidermis formation, one locus (S_Unigene3026_c0_seq1_ 3726) was related to epidermis formation and muscle formation, three loci (S_Unigene512_c0_ seq1_971, S_Unigene512_c0_seq1_5524, and S_Unigene512_c0_seq1_5912) were related to both vascularization and muscle formation, one locus (S_Unigene11849_c0_seq1_804) was related to the formation of blood vessels and skin, and the remaining loci were related to vascularization. Finally, we selected the myosin heavy chain gene, which was related to S_Unigene512_c0_seq1_ 971 among the 26 SNPs, for the detection of muscle development in the four species.

mRNA expression of non-muscle myosin heavy chain II in different tissues and different species
In mammals, NMHC II has three isoforms, which are referred to as NMHC IIA, NMHC IIB and NMHC IIC [28,29]. However, Xenopus has two isoforms, II-A and II-B, and does not appear to have an II-C isoform [30]. In invertebrates, Drosophila has only a single isoform of NH II [31,32]. According to the taxonomic placement of Onchidiidae, we speculated that Onchidium contained at least one isoform. The total length of the Os-NMHC sequence is 6057 bp. Furthermore, the specific expression of a gene in tissues is typically related to its function in those tissues, and different tissues reveal the different adaptions of a species. To investigate tissue-specific expression, the levels of mRNA expression in dorsal skin, ventral skin, lung sac, ganglion and ventricle samples from the four species were quantified by qRT-PCR. The SNPs reflected genetic differences among the four species of Onchidiidae. However, the expression differences of genes related to phenotype remained unknown. According to our histological study of Onchidiidae and the transcriptome data, we determined that Os-NMHC reveals the adaption from seas to wetlands in the four species of Onchidiidae. Onchidium reevesii shows apomorphic characters, as is evident from its ability to live in more complex environments [8]. We obtained the full-length cDNA, submitted it to the GenBank database and obtained an accession number (KU663401). The expression of Os-NMHC in various tissues of adults of the four species of Onchidiidae was determined (Fig 5). The O. reevesii has developed lung-sac observed among the four species; accordingly, this species presented high expression of the Os-NMHC gene in the lung-sac. In Platevindex mortoni, Os-NMHC exhibited the highest expression level in ganglion tissue, whereas Os-NMHC expression was not detected in most tissues of Peronia verruculata. The expression levels in the same tissue differed among the species (P<0.05).

mRNA expression of myosin heavy chain in different tissues and different species
We cloned the MyHC gene (GenBank accession number: KU550708) of the four species and compared the expression levels in four types of tissues among the four species. The results showed that the full length of MyHC is 7566 bp. The expression level of MyHC was highest in Onchidium reevesii and lowest in Paraoncidium reevesii, and this difference was significant (P<0.05). The highest expression in the ventral skin and foot was observed in Platevindex mortoni. In the lung sacs, P. mortoni had the highest expression of MyHC, closely followed by O. reevesii, and P. reevesii showed the least expression.  We then analyzed the relative expression of the MyHC gene in three different tissues from the four Onchidiidae species to examine associations with living habitat (Fig 6). O. reevesii and P. mortoni are mainly terrestrial and burrow in mud or climb rocks to avoid the tide. However, P. reevesii and P. verruculata are mainly aquatic, and their movement requirements are lower. O. reevesii had high expression levels of MyHC in the dorsal skin, foot and lung sac, which are suited to their terrestrial habitat. P. mortoni expressed high levels in the foot, and this species Comparison of four species of Onchidiidae frequently climbs trees. The epidermis of P. verruculata showed the highest level of keratinization observed among the four species; accordingly, this species presented high expression of the MyHC gene in the dorsal skin. We speculated that the expression level of this gene is related to the species' respiration ability, moisture retention ability and defense capacity.

Discussion
Skin is an important respiratory organ for onchidiids, and skin with lower keratinization has higher permeability, which facilitates breathing in Onchidiidae species. In contrast, skin with higher keratinization helps retain moisture and protect against predators.
The skin of Onchidium reevesii, which is mainly terrestrial, has relatively weak respiratory function. The epidermis of this species is thick and functions in retaining moisture and protecting against predators, in accordance with its terrestrial characteristics [33,34]. The dorsal skin of the mainly aquatic Paraoncidium reevesii is thin and highly permeable; thus, its respiratory function is strong. Another aquatic species, Peronia verruculata, has a higher level of keratinization of the epidermis, but its gill skin is thin and suitable for breathing when submerged ( Table 2). For Platevindex mortoni, skin thickness and the number of blood sinuses are intermediate among the four species. This species lives mostly in the supratidal zone and mudbank, can remain in the sea for long periods, and can climb trees.
The secretions from the mucous glands are slimy and smooth, which can reduce the friction between skin and water and facilitates gas exchange and ion transport [35]. Moreover, the dense distribution of blood sinuses is a hallmark trait of aquatic species [36]. There is only a small amount of blood sinuses in the stratum spongiosum of Onchidium reevesii, whereas the blood sinuses in the remaining species are abundant.
The sequence of the diameters of the sac room and the small room are also showed differences among the four species. Dayrat called the respiratory organ of Onchidium vaigiense a lung sac, which is similar to the breathing bag of limacines [37]. The developed lung sacs of amphibians are more suitable for terrestrial life [38]. The efficiency of lung sac respiration depends on the sac's superficial area. Thus, species with larger superficial areas have stronger respiratory capacities. Because of its well-developed reticular diaphragm, thin connective tissue, rich blood capillaries and largest superficial area of the lung sac among the four species, O. reevesii has the strongest respiratory capacity for wetland living among the four species. In conclusion, the degree of lung sac development in the four species of Onchidiidae decreased in the order Onchidium reevesii, Peronia verruculata, Platevindex mortoni and Paraoncidium reevesii.
The evolution of Onchidiidae species and their amphibious features is reflected in their morphological characteristics. Transcriptome sequencing data allow the study of genes related to the morphological differences and amphibious features of these typical amphibious mollusks. Our data provide the best transcriptomic resource currently available for these four species. The transcriptome data were obtained by Illumina HiSeq TM 2000 sequencing, and the sequences were assembled and functionally annotated. Based on the annotated unigenes, GO and KEGG assignments were determined. This study established an excellent resource for future genetic or genomic studies of Onchidiidae variation and a platform for functional genomics and comparative genomic studies of mollusks.
In this study, SNP loci for O. reevesii were developed based on comparisons of its transcriptome sequences with those of the other three species. The different species adopted different strategies during their evolution. These loci of O. reevesii differed from those of Peronia verruculata, Paraoncidium reevesii and Platevindex mortoni and were related to vascularization and the formation of muscle and cuticle. These loci mutated readily and might reflect strong directional selection, which is important for understanding differences in respiratory traits and amphibious features among species. These loci might indicate reference genes and should be verified in other species. The SNPs might have undergone directional selection leading to adaptive evolution in Onchidiidae [39]. We selected the myosin heavy chain gene from among the SNP loci to investigate the different adaptations of the four species.
The dorsal skin plays significant roles in defending against predators and protecting against moisture loss. The developed ventral skin and foot skin are important for locomotion, and NMHC II can influence axon growth [40]. In some species, such as Drosophila [41], NMII is related to dorsal skin and is suggested to be involved in epidermal barrier functions [42]. Therefore, relevant to studies of epidermis adaption, NMHC II plays an important role in skin development. Our study revealed that the species that are mainly terrestrial (i.e., Onchidium reevesii and Platevindex mortoni) have a developed epidermis that can retain water and provide defense against predators. Their developed ventral skin and foot are suitable for rough terrain in the terrestrial environment. We found pronounced NMHC II expression in the skin (dorsal skin, ventral skin and foot skin) in the species that are mainly terrestrial, and we also found that P. mortoni showed high expression in the ventricle and ganglion. Previous studies found that the mutation of NMHC II affects the development of the heart [43,44]. Because the heart and ganglion are closely associated with feeding habits, this evidence explains how both O. reevesii and P. mortoni had the ability to adapt to terrestrial environments [45,46]. The observed higher-level expression of Os-NMHC in the respiratory and locomotive organs (dorsal skin, ventral skin, foot skin and lung sac) suggests that species that are mainly terrestrial (e.g., O. reevesii and P. mortoni) can easily adapt to complex land environments.
The myosin heavy chain protein is associated with muscle and is a suitable protein for analyzing muscle adaptation. Thus, its level of expression is associated with species habitat. The myosin heavy chain gene was thus selected from among the SNP loci for further examination. O. reevesii and P. mortoni mainly live in wetlands and must burrow in mud or climb rocks to avoid the tide. P. reevesii and P. verruculata, which are mainly aquatic, do not require strong terrestrial locomotor abilities. Moreover, the environment in which P. reevesii lives is more complex than the environments of the remaining species, and P. reevesii can climb trees. Therefore, P. reevesii has the strongest requirement for a developed foot, and its expression level of MyHC in its foot is much higher than the expression levels of this gene in the foots of P. reevesii and P. verruculata.

Conclusion
Onchidiidae constitute an interesting group of invertebrates that occupy habitats ranging from seas to wetlands, showing a gradual distribution. Onchidiidae is an intermediary form that connects the invertebrates on land to those in the sea. In addition, this amphibious mollusk family offers a model for understanding histological and genetic changes associated with the water-toland transition of invertebrates. Members of this group successfully made the transition from aquatic to terrestrial living and breathe with lung sacs, skin and gills to adapt to an amphibious lifestyle. Our histological analysis of respiratory organs from four onchidiids provides insights into their different adaptation. The 26 SNP loci that were selected can be used in further studies of respiratory traits and amphibious features. Moreover, we found that the relative expression levels of two genes are associated with the locomotor traits of the four species. We hope that this study provides a valuable reference point and a source of inspiration for future studies.