Comparative analysis of the ribosomal DNA repeat unit (rDNA) of Perna viridis (Linnaeus, 1758) and Perna canaliculus (Gmelin, 1791)

Mussel material and DNA extraction

Perna viridis samples were collected offshore of Leizhou Bay (20.968614°N, 110.325743°E) in China, and P. canaliculus samples were collected from Pegasus Bay (43.312595°S, 172.854629°E) in New Zealand. Three individuals from each species were used in this study. Total genomic DNA was extracted from muscle tissue using the TIANamp Marine Animals DNA Kit (TianGen Biotech Co., Ltd, Beijing, China) according to the manufacturer’s instructions.

PCR amplification, cloning and sequencing

The rDNA fragments were amplified using the primers listed in Table 1. The primers used to amplify the 18S and 28S rRNA genes were designed based on partial sequences corresponding to the 18S rRNA genes of P. viridis (EF613234) and 28S rRNA genes of Geukensia demissa (AY145405), respectively. The primers used to amplify the 18S–ITS–28S region were designed from the 3′end of P. viridis (EF613234) and 5′end of G. demissa (AY145405), while the IGS region were designed based on the 3′end of G. demissa (AY145405) and 5′end of P. viridis (EF613234). Each PCR had a total volume of 25 μl and contained the following: one μL of DNA template, 0.5 μL of each primer (10 μmol/L), 12.5 μL of 2 × TransTaq PCR SuperMix and 10.5 μL of H₂O. The reactions were initially denatured at 95 °C for 3 min and subjected to 35 cycles of denaturing at 95 °C for 45 s; annealing at 52 °C for 18S and 28S rRNA, 50 °C for 18S–ITS–28S and 56.5 °C for 28S–IGS–18S for 30 s, and extension at 72 °C for 1–3.5 min (1 kb/min). Finally, all reactions were subjected to a final extension at 72 °C for 7 min. The PCR products were electrophoresed on 0.8% agarose gels, extracted and purified. The purified PCR products were cloned using the pUCm-T Vector Cloning Kit and SanPrep Column Plasmid MiniPreps Kit (Sangon Biotech Co., Ltd., Shanghai, China) according to the manufacturer’s instructions. Three positive recombinant plasmid colonies corresponding to each amplified region were picked, cultivated and sequenced by Beijing Ruiboxingke Biotechnology Co., Ltd. (Beijing, China).

Sequence analysis

The complete rDNA unit was assembled using DNAMAN software (Kaukonen et al., 2000), and both sequences have been deposited in the GenBank database (accession numbers: MK419104–MK419109). The boundaries of each region were confirmed using the NCBI database and BLAST software (Johnson et al., 2008). Multiple sequence alignment was performed using T-Coffee software (Notredame, Higgins & Heringa, 2000) with manual adjustments. The general molecular features of P. viridis and P. canaliculus rDNA were calculated using MEGA 7.0 (Kumar, Stecher & Tamura, 2016). Sub-repeat (SR) fragments and inverted repeats of the IGS region were analysed with the Tandem Repeats Finder (Benson, 1999) and Unipro UGENE (Okonechnikov, Golosova & Fursov, 2012), respectively. Predictions of putative transcription initiation sites (TIS) and transcription termination sites (TTSs) in the IGS regions were based on a comparative analysis of the sequences with data for various marine animal species from the literature (Ki, Kim & Lee, 2009; Luchetti, Scanabissi & Mantovani, 2006; Chae et al., 2018; Guo et al., 2018). Cytosine-guanine (CpG) islands were identified using CpGPlot (Polanco & De La Vega, 1994). The sequence identity was confirmed using DNAMAN software. The genetic distance between P. viridis and P. canaliculus was calculated with MEGA 7.0 software according to a Kimura two-parameter model. Nucleotide diversity (π) was calculated using the software DnaSP 5.10 (Librado & Rozas, 2009).

For phylogenetic reconstruction, the ITS sequences of P. viridis and P. canaliculus were first determined. Next, nine complete or near-complete sequences from the following Mytilidae species were obtained from the GenBank database: Mytilus edulis (AY695798), M. unguiculatus (former M. coruscus) (MK157201), M. galloprovincialis (JX081670), M. trossulus (JX081669), Aulacomya maoriana (former Aulacomya atra maoriana) (DQ924557), P. perna southern Indian population (former P. indica) (JQ622200), P. perna Omani population (KC692037), P. perna Mauritanian population (former P. picta) (DQ924548) and Modiolus rectus (EF035114). Maximum likelihood (ML) and neighbour-joining (NJ) phylogenetic trees based on ITS sequences were constructed with 1,000 bootstrap replicates using MEGA 7.0 software and the GTR+G nucleotide substitution and Kimura two-parameter models, respectively. Crassostrea ariakensis (EU252081) was selected as the outgroup.

Complete rDNA sequences of P. viridis and P. canaliculus

The complete rDNA units of P. viridis and P. canaliculus were determined after sequencing, comparison and assembly. Figure 1 presents a scheme of the structural organisation of the rDNA repeat unit, which includes the following elements in order: 18S–ITS1–5.8S–ITS2–28S–IGS. The complete rDNA units of P. viridis and P. canaliculus varied from 8,432 to 8,616 bp and from 7,597 to 7,609 bp, respectively.

The length, GC content, pairwise identity and variable sites of each region are shown in Table 2. An analysis of the P. viridis and P. canaliculus sequences revealed low levels of genetic variation in the 18S, 5.8S and 28S rRNA regions, which had pairwise identities of 99.56%, 100% and 98.94%, respectively. The variable (V) sites in 18S and 28S rRNA made respective contributions of only 0.83% and 1.22%, indicating that the coding regions (18S, 5.8S and 28S rRNA) were even more highly conserved than the complete rDNA unit (6.25%). Differences between P. viridis and P. canaliculus were observed mainly in non-coding regions (ITS1, ITS2 and IGS). In intra-species comparisons of three individuals per species, the sequence similarities of ITS1 and ITS2 in P. viridis were 99.38% and 99.87%, respectively, and these regions differed by only five and one base, respectively. Among P. canaliculus individuals, both regions had sequence similarities of 100%. The respective sizes of ITS1 and ITS2 were 269 and 263 bp in P. viridis and 304 and 258 bp in P. canaliculus. Differences in the lengths and variable sites of the ITS1 and ITS2 regions resulted in relatively low sequence identities (66.99% in ITS1, 67.03% in ITS2).

Comparisons of the IGS sequences at various levels revealed that the highest (91.84% for P. viridis, 98.31% for P. canaliculus) and lowest similarities (33.32%) were found at the individual and inter-species levels, respectively, and the variable sites between the IGS sequences of the two species accounted for 15.65% of all sites in this region (Tables 2 and 3). The sizes of the IGS sequences ranged from 1,399 to 1,411 bp and 2,267 to 2,451 bp in P. canaliculus and P. viridis, respectively. This result suggested that the IGS sequences of them had large differences in length and a relatively high level of sequence heterogeneity. Due to heterogeneity in the length of IGSs at the individual and inter-species levels of the genus Perna, the longest IGS sequence of each species was chosen as a representative to outline the canonical structural organisation of the IGSs from P. canaliculus and P. viridis. The detailed analysis of the molecular structure of IGS sequences of them contains six distinct regions: TTS, TIS, NTS, ETS, SR and CpG islands.

Beginning and ending sequences of transcripts, ETS and NTS

A poly (T) tract (5′-TTTTCGTTTGCCTTTTTCGTTCCTTTTTTTTT-3′) was identified in the P. viridis IGS region from 171 to 202 bp. The sequence 5′-TTACTTGT-3′ was detected in the P. canaliculus IGS region from 120 to 127 bp. Both sequences were considered putative TTSs. Accordingly, the positions of the 3′ETS sequences in the P. viridis and P. canaliculus 5′IGS regions were determined to correspond to bp 1–202 and bp 1–127, respectively. The sequence ‘TTATTATGTGGAGTGGG’ was considered a putative TIS for RNA polymerase I, and the +1 position of the TIS was the C/T nucleotide in P. viridis and P. canaliculus. The 5′ETSs were localised between the putative TIS and the beginning of the 18S rRNA gene and had lengths of 581 bp and 500–503 bp in P. viridis and P. canaliculus, respectively. Non-transcribed spacers comprised the remainder of IGS. The length, nucleotide diversity, intra- and inter-species sequence identity of 3′ETS, NTS and 5′ETS in P. viridis and P. canaliculus were shown in Table 4, comparision with 3′ETS and 5′ETS, NTS had lower intra- and inter-species sequence identity and higher nucleotide diversity. The high level of individual, intra- and inter-species divergence within the NTS sequences was attributed to lower functional constraints. In contrast, the nucleotide sequences were more conserved in the ETS region than in the NTSs, despite differences in the length of the former between P. viridis and P. canaliculus.

Sub-repeat regions and CpG islands in the IGS sequence

Next, the DNA repeats in the IGS sequences were analysed to explore the IGS genetic structure in greater detail. The identification of the numbers and lengths of the repetitive sequence motifs occurring in P. viridis and P. canaliculus enabled the definition of three different SRs with lengths of 40–131 bp (Figs. 1 and 2). A comparison of the IGS regions indicated an extremely high level of sequence diversity within the SR region. The P. canaliculus sequence included only one type of SR (length: 131 bp) which was repeated 2.9 times, and three poly (A) tracts in the NTS area. The P. viridis sequence included two repeat motifs: SR1 (length: 93 bp), which was repeated 9.3–11.5 times within the NTS, and SR2 (length: 40 bp), which was repeated 2.5 times within the 5′ETS. Differences between individuals of the same species were due to SR copy number variation. No inverted repeats were detected in either Perna species.

The GC content of the IGS region was slightly higher in P. viridis than in P. canaliculus (Table 2). Similarly, the number and lengths of the CpG islands were greater in P. viridis than in P. canaliculus. In P. viridis, three CpG sites were found to cover nearly the entire IGS sequence (lengths: 440, 1,078 and 537 bp). By contrast, the P. canaliculus IGS region included two CpG islands located at the NTS and 5′ETS (lengths: 361 and 484 bp).

Phylogenetic tree analysis based on ITS sequences

Maximum likelihood and NJ phylogenetic trees were constructed based on the ITS sequences of Mytilidae species (Fig. 3). Clustering pattern analysis revealed that the NJ and ML trees were identical, both including three major clusters, and species of the same genus clearly clustered together with relatively high supporting values. For both P. viridis and P. canaliculus, the individuals were classified into single clusters with high sequence identity. All Mytilidae species included in the phylogenetic analysis belong to the subfamily Mytilinae, except M. rectus which belongs to the subfamily Modiolinae. During tree formation, M. rectus and A. maoriana initially clustered together, and then clustered with Perna and Mytilus species.