Molecular organization of the 5S rDNA gene type II in elasmobranchs

ABSTRACT

The 5S rDNA gene is a non-coding RNA that can be found in 2 copies (type I and type II) in bony and cartilaginous fish. Previous studies have pointed out that type II gene is a paralog derived from type I. We analyzed the molecular organization of 5S rDNA type II in elasmobranchs. Although the structure of the 5S rDNA is supposed to be highly conserved, our results show that the secondary structure in this group possesses some variability and is different than the consensus secondary structure. One of these differences in Selachii is an internal loop at nucleotides 7 and 112. These mutations observed in the transcribed region suggest an independent origin of the gene among Batoids and Selachii. All promoters were highly conserved with the exception of BoxA, possibly due to its affinity to polymerase III. This latter enzyme recognizes a dT₄ sequence as stop signal, however in Rajiformes this signal was doubled in length to dT_8. This could be an adaptation toward a higher efficiency in the termination process. Our results suggest that there is no TATA box in elasmobranchs in the NTS region. We also provide some evidence suggesting that the complexity of the microsatellites present in the NTS region play an important role in the 5S rRNA gene since it is significantly correlated with the length of the NTS.

KEYWORDS:

• Internal Control Region
• Microsatellites
• Non-transcribed Spacer
• RNA secondary structure
• TATA-box

Abbreviations

TFIIIA	=	Transcription factor IIIA
ICR	=	Internal control region
NTS	=	Non-transcribed Spacer

Introduction

The organization of ribosomal genes includes a transcribed followed by a non-transcribed region arranged in tandem within the nucleolar organizing region (NOR)¹ whose numbers vary among organisms.^2,3 In most eukaryotes, these genes are classified in 2 different groups: 45S and 5S rDNA. Following a typical organization, the 5S rDNA contains a 120 bp transcribed region that is highly conserved due to its essential function in stabilizing the ribosomal structure and enhancing the peptidyl transferase activity.⁴ This transcribed region is adjacent to a Non-transcribed Spacer (NTS) that varies in length among species.^1,5 This NTS region has a higher mutation rate than the transcribed region.¹ Insertions and deletions, rather than nucleotide substitutions, is believed to be the main cause of the observed variation.⁶

The 5S rDNA codes for a small structural RNA and is classified as a class III gene. This implies that it is transcribed by RNA polymerase III. RNA Polymerase III promoters are divided into 3 types according to their organization. The promoter for 5S rDNA is of type 1 and consists of BoxA, an Internal Element (IE), and BoxC; that is found within the transcribed region.⁷ These three elements are grouped in the Internal Control Region (ICR)⁸ and in elasmobranchs (shark, rays, and skates) have been reported to be respectively at positions 50 – 64, 67 – 72, and 80 – 97 bp of the transcript.^9,10 The transcribed region also contains a termination signal consisting of a poly-T relic, that is still recognized by RNA polymerases III during transcription. Another feature of 5S rDNA is that a TATA box is not strictly required for transcription. However, in mammals, it has been shown that there is a TATA sequence in the NTS region that may play an important role in regulating 5S rDNA expression.^5,12 In some fish species and in elasmobranchs, 5S rDNA has a TATA-like sequence, suggesting that this sequence influences gene expression.^{9,10,13,14,15,16} This TATA-like element in fish has not been experimentally investigated yet. In Xenopus borealis a deletion of 55 nucleotides in the upstream 5′ end of the 5S rRNA only expressed a part of the gene, showing that this region influences in the transcription process.¹⁷

Like many ribosomal genes, 5S rDNA has been shown to follow the model of concerted evolution.¹⁸ This model proposes that paralogous genes in a species are more similar among themselves than with orthologous genes across species as a result of homogenization.^19,20,21 This homogenization has been attributed to the replicative transposition, gene amplification, and gene conversion.^22,23 It has also been proposed that the unequal crossing over among the repetition units play an important role in this process.²⁴ Walsh,²⁵ using computer simulations and taking into account more parameters, showed that this process deleted tandem arrays. Later, this was supported by Stephan²⁶ who evaluated the structural patterns of tandem repetitions with the simulations of different recombination rates. Their results demonstrated that at high recombination rates, the repetitive structures are distributed in a regular pattern.

In elasmobranchs (as well as in several bony fishes²⁷) 5S rDNA occurs frequently in 2 forms: type I and II. Both types are differentiated by the length and the sequence of its NTS region. The longest of the 2 is designated as type II. The evolutionary history of type II gene suggests that it might have originated in the elasmobranch lineage by genome duplication. This further implies that it is paralogous to 5S rRNA type II in bony fish.²⁸ The aim of this study is to compile and analyze the molecular organization and the secondary structure of 5S rDNA previously as type II in elasmobranchs.

Results and discussion

Transcribed region: mutations in the ICR components and termination signal of the

RNA polymerase III

PCR amplification of the 5S rDNA type II in Urotrygon with the primers 5SRajid amplified a unique band of 1,5 Kbp (supplementary Fig. 1) and its partial sequencing revealed that this band is composed of the 5S rRNA gene as well as its associated NTS. The 5S rRNA segment contains all the elements of the ICR and the terminal signal of the polymerase III. Since the transcribed region has a functional role within the organism, it is expected to be strictly conserved. However, there are point mutations in the transcribed region covering both the inside and outside of the ICR. The BoxA has 7 positions with mutations (positions 51, 52, 56, 59, 60, 61, and 64). From these positions the nucleotide 60 is the one with highest inferred information (∼1 bit). The BoxC only shows 2 mutations (positions 84 and 92), while the IE has a single mutation with high information (∼1.8 bit) corresponding to the nucleotide 68. Outside the ICR 15 mutations were found, with the nucleotides 18, 24, 35, and 112 showing low information (≤ 1 bit) (Fig. 1).

Figure 1. Promoter and terminator elements description and entropy index for the transcribed region of the 5S rDNA. A) The promoter of the 5S rDNA contains the elements BoxA, the Internal Element, and the BoxC, altogether called Internal Control Region (ICR). The RNA polymerase recognizes a poly dT as stop signal. B) The entropy index for the 5S rRNA using all the 91 sequences retrieved from the GenBank. However, due the fact that the first 43 nucleotides in Potamotrygon motoro, Potamotrygon falkneri and Paratrygon aiereba transcript region sequences were missing, these were assumed to be equal to those in Urotrygon.

Unlike BoxC, the IE, and BoxA show a higher number of mutations and ratio between the number of mutations over length (BoxA = 7/14, IE = 1/5, and BoxC = 2/17). This might be due to the effects of the IE and BoxA in the affinity to the Transcription Factor IIIA (TFIIIA) compared to the BoxC.^29,30,31 It has been shown that BoxC interacts with the first 3 amino-terminal fingers of the TFIIIA, and that these are required for affinity to the 5S rDNA.^32,33 This may be the reason why BoxC is highly conserved among individuals. In Rajiformes the termination signal of the RNA polymerase III (poly-T tail) is longer (up to 4 dTs longer) than in other clades (up to 4 nucleotides as the expected range for eukaryotes³⁴). Some studies have shown gene transcripts (like the tRNA^Lys in rabbits and X. laevis, and the VA RNA-II in adenovirus) that contain extensive dT tail that are more efficient at terminating transcription.³⁵ Similar cases have been reported in Saccharomyces pombe and S. cerevisiae, where effective transcription required 5 and 6 thymidine fragments, respectively.³⁶ Apparently, a higher content of dT nucleotides in the termination signal make the transcription more robust limiting the influence of surrounding sequences. This allows that adjacent dC and dG nucleotides improve the effectiveness of the termination while the dA nucleotide decreases it.³⁷ However, there are some exceptions like the tRNA^Lys of X. laevis where an extension of residues dA downstream to the termination signal enhances its effectiveness.³⁸ In Rajiformes, both upstream and downstream sequences in the termination signal are similar to that of other species. This makes it difficult to conclude that the long tail of poly-T in Rajiformes is due to surrounding sequences that interfere with the effectiveness in the termination signal of the RNA polymerase III. Also, the poly-T tail has several adjacent nucleotides dC and dG that increase the effectiveness of the transcription termination.³⁷ Considering that RNA polymerase III often ignores the original termination signal and recognizes the posterior poly-dT cluster,³⁹ it is possible that the significant excess of dT residues (p < 0.01; Table 2) in the original termination signal in Rajiformes, with respect to other clades, is an adaptation to a better efficiency of transcription termination.

Table 1. List of the 5S rDNA type II gene sequences obtained from 14 elasmobranchs and 3 Petromyzontiformes species retrieve from the NCBI database.

Species	Number of sequences	NCBI Accession code	Order
Entosphenus japonicus	2	X04308-X04309	Petromyzontiformes
Lampetra reissneri	1	X13038	Petromyzontiformes
Lethenteron camtshaticum	1	D00076	Petromyzontiformes
Rhizoprionodon lalandii	19	FJ517239-FJ517257	Carcharhiniformes
Rhizoprionodon porosus	16	FJ517223-FJ517238	Carcharhiniformes
Galeocerdo cuvier	1	FJ539131	Carcharhiniformes
Scyliorhinus caniculus	1	M24954	Carcharhiniformes
Alopias superciliosus	1	FJ539130	Lamniformes
Potamotrygon motoro	6	JF792331-JF792336	Myliobatiformes
Potamotrygon falkneri	3	JF92328-JF92330	Myliobatiformes
Paratrygon aiereba	5	JF92323-JF92327	Myliobatiformes
Raja asterias	5	DQ020522-DQ020526	Rajiformes
Raja clavata	6	DQ020527-DQ020532	Rajiformes
Dipturus oxyrinchus	6	DQ020541-DQ020546	Rajiformes
Raja miraletus	4	DQ020537-DQ020540	Rajiformes
Raja polystigma	4	DQ020533-DQ020536	Rajiformes
Raja montagui	1	AY278250	Rajiformes

Table 2. Contingency table of Chi-square to test the equality of the stop signal in all the elasmobranchii individuals studied here.

	Observed		Expected		Chi Squared
	Present	Absent	Present	Absent	Present	Absent	Total	p-Value
dT7	26	0	7,9	19,1	39,6	17,2	56,7
dT4	0	60	18,1	41,9	17,2	7,4	24,6
Total	26	60	26,0	60,0	56,7	24,6	81,3	<0,01

Changes in the transcribed region and secondary structure through elasmobranch

Evolution

The low information inferred for the positions 18, 24, 35, 60, and 112 is mainly due to variability in some taxonomic units. Mutations C18T, G60A, and T112M (where M = A or C) are characteristic of a division at a higher taxonomic level: Batoidea (Rajiformes and Myliobatiformes, in this study) with respect to Selachii (Carchariniformes and Lamniformes, in this study). However, the synapomorphic mutation of the position 112 is absent in the genus Galeocerdo. This might be due to the under-representation of this genus in the sampling (1 individual). Mutations A24C and T35C are synapomorphic in Myliobatiformes with respect to Rajiformes (Fig. 2). The synapomorphic mutations C18T, G60A, and T112M agree with the hypothesis of an independent evolution of Batoideans and Selachians based on the fossil record⁴¹ and also congruent with several molecular phylogenetic studies.^{42,43,44,45,46} However, Compagno's hypothesis,⁴⁷ that batoids are derived from Pristiophoriformes, cannot be tested with our data since the sharks of the superorder Squalimorphi are not represented in our sampling.

Figure 2. Evolutionary interpretation of the synapomorphic mutation events. (*)The synapomorphic mutation of the nucleotide 112 is absent in Galeocerdo cuvier. This might be due to the underrepresentation of this genus or a back mutation in that position.

The mutations at positions 18, 60 and 7, 112 form pairs in the secondary structure proposed by Barciszewska et al.⁴⁸ (Fig. 3A), and is similar to our secondary structure consensus (Fig. 3B) with the following exceptions: a) mismatch in the basepair (7,112) in sharks except in G. cuvier, b) within the loop C the basepairs (33,42) and (34,41) are joint, c) the bulge located at position 84 in Barciszewska et al. is now located at position 83, and d) basepair (80,96) mismatches in all the structures obtained. The aforementioned 7–112 mismatch (feature a, Fig. 3B) is located in helix I. Therefore, a large effect is expected on the transcription level since the TFIIIA binds to this site.^50,51 In addition, this feature is highly conserved in Selachii. An experimental validation of this claim is needed to assert this effect.

Figure 3. Secondary structure consensus of the 5S rRNA type II predicted by RNAalifold. The structure A) is proposed by Barciszewska et al.⁴⁸ which corresponds to the general model of 5S rRNA for eukaryotes. The structure B) is the consensus structure predicted for the 5S rRNA type II by RNAalifold. The main differences are those marked with red arrows: a) positions 7 and 112 with no pair in all sharks but G. cuvier; b) all the sequences fit well with the basepair 34C:41G and 33U:42A which reduce the loop C size; c) the bulge at the position 84 (with base dU in this study) proposed by Barciszewska et al.⁴⁸ is now located at the position 83 (dA). According to these authors the bulge located in the helix IV must be a purine; d) positions 80 and 96 are not joined in all the sequences.

The results obtained by mutual information show 5 pairs of positions with scores above the critical value (M_ijc ≈0.212), these positions are (60,18), (112,18), (92,60), (112,60), (60,24) with M_ij equal to 0.99, 0.51, 0.24, 0.45, and 0.23, respectively. Moreover, the covariance-like measure shows 3 pairs of positions with values above its critical value (C_ijc ≈0.172), these positions are (18,51), (18,60), and (35,60) with C_ij equal to 0.31, 0.49, and 0.44, respectively. Although there is some evidences of the formation of base triples in RNA structures,⁸² it is well known that the 5S rRNA structure form only base pairs.⁸³ Given that one nucleotide can join only to another one nucleotide, we opted to choose only the pair position (18,60), since this is the only pair common to the 2 analyses and has the highest score. We used this information to constrain the secondary structure prediction for each haplotype. This resulted in 25 different secondary structures which can be classified into 3 groups (Fig. 4). Group α represents 84% of all the sequences studied and have a structure similar to the consensus structure, with 3 arms. The structure of the group β have 4 arms and represents 5% of all sequences. The group γ comprises 11% of all the sequences and has a structure with 4 arms arranged differently from those of group β. Despite the amount of different predicted secondary structures, the RNA sequence fit well to the consensus secondary structure (Fig. 3B) with exceptions in some positions due to SNP's. This lower accuracy in the 5S rRNA secondary structure prediction with respect to other RNA molecules has been reported previously⁵⁶ (and references therein). This could be due to possible interactions that are needed for the proper folding of the 5S rRNA. For example, within the group γ (Fig. 4), the sequence FJ517193.1 has a structure with a free energy of -37.41 kcal/mol. This is notably different from the others due to the mutation G116A which seems to destabilize helix I. When this mutation is reverted, the resulting secondary structure is similar to those of the group α and its free energy is reduced to -43.60 kcal/mol. However, experimental validation of this hypothesis is needed. In addition, we also marked those sequences whose promoter elements are different to the consensus sequence of metazoan 5S rRNA in the tree of structures (supplementary Fig. 2). We found not only that these sequences are interspersed in the tree (among the 3 groups), but also those sequences with promoter elements that fit the consensus sequences. These results suggest that SNP's both inside and outside the ICR can change the pattern of folding into another conformation with lower free energy.

Figure 4. Secondary structural groups of the 5S rRNA type II defined by neighbor joining. The results show 3 different structural groups: α, β and γ. The group α is composed of structures that are similar to our consensus secondary structure of the 5S rRNA type II, while the groups β and γ are composed of structures different to our consensus. Red numbers represent the bootstrap values.

The compensatory mutation between nucleotides 18 and 60 indicates that helix II in the secondary structure is important for transcription. This may be because it is the region where TFIIIA's 7 zinc fingers attach.^49,50 In fact, it has been proven that a minimum fragment of 5S RNA that includes the helices I, II, and V, and the loops A and E, is sufficient to bind to TFIIIA.^50,51 This is an essential component of the initiation complex in the transcription of RNA polymerase III.^52,53,54 Nucleotides 24, 35 and 40 are not tightly bound to the secondary structure and can mutate more freely. Also the loops B and C containing these mutations don't have an important role for the union of the TFIIIA's zinc fingers to the 5S rRNA molecule.⁵⁵

NTS region

TATA Box - The NTS region of the 5S rDNA type II is more conserved within species than expected considering the overall substitution rate over the wider phylogenetic range.^57,58 This can be explained either by concerted evolution or because of strong selective pressures due to the regulation of the gene.^5,59 In bony fishes, a TATA-like box was found between position -20 and -30 upstream from the transcription start site.^13,14,60 This suggests a possible influence on the level of transcription. However, despite the presence of a TATA-like box in Rhizoprionodon (with an AATT motif starting at nucleotide -27)¹⁰ and in Rajiformes (as TACA at nucleotide -25),⁹ in both Urotrygon species there is not a TATA-like element within the positions −30 and −20 (5′-TGTGAAGAGGC-3') similar to those reported previously in Rhizoprionodon and Rajiformes. This discrepancy in both sequence and location with respect to a TATA-like box among the 3 orders may indicate that these TATA-like boxes are not really essential factors. Furthermore, since the promoter of 5S rDNA is of type 1, all of its elements are internal and therefore not necessary for transcription.¹¹ However, more studies are needed to determinate whether these TATA boxes have some level of influence on the transcription process. There are hybrid promoters like the U6 snRNA gene in S. cerevisiae that contain both internal promoters and a TATA box^61,62 which contributes to better transcription.^63,64,65

Microsatellites - The different microsatellites arrangements between genera at the NTS region are the main contributors to the variation in the NTS length. It has been shown that there is a correlation between the length of the NTS region and the number of motif repetitions within each region of microsatellites in Rajiformes.⁹ But in addition, it is possible that there exists a correlation among the number of microsatellite cluster regions and the length of the NTS region. This process is possible because the NTS shows several clusters of repetitions such as (dNdN)₂ (see supplementary Fig. 3). To confirm this correlation, a linear regression was performed. However, missing data in Urotrygon, Potamotrygon and Paratrygon might bias the analyses. To test whether the missing region could be safely ignored, we performed the Clark-Evans and the Diggle-Cressie-Loosmore-Ford (DCLF) test in Rhizoprionodon, Raja and Dipturus to check if the microsatellites were distributed randomly, regularly or aggregatively. If the distribution is random or regular, and assuming that the NTS of Urotrygon, Potamotrygon and Paratrygon behaves similarly, the regression won't be affected by the absence of the fragments. Our results show that indeed the distribution of the microsatellite is random in some species, regular in others, but never aggregated (Table 3). Therefore, the partial sequence of the NTS of Urotrygon, Potamotrygon and Paratrygon species can be used for the linear regression.

Table 3. Statistical results for the distribution of the microsallite region inside the NTS. Significant p-values for R. montagui and R. miraletus indicate that its distribution is regular. For other species the distribution is random.

Species	Clarc-Evans testtwo tail (p-value)	DCLF(p-value)
R. lalandii	0.67	0.87
R. porosus	0.22	0.60
R. polystigma	0.28	0.73
R. montagui	0.03	0.65
R. miraletus	0.03	0.87
R. clavata	0.55	0.30
D. oxyrinchus	0.29	0.80
R. asterias	0.26	0.18

The linear regression analysis shows a strong relationship between these 2 variables (r = 0.95, p − value < 0.01; Fig. 5), even if both species of Urotrygon are excluded from the analysis (r = 0.75, p − value < 0.01). With this, we speculate that the microsatellite clusters play a functional role in transcription where a mechanism of unwinding is providing and used during the transcription process.⁶⁶ However, according to Pinhal et al.¹⁰ this is not possible because the NTS region of the Rhizoprionodon genus does not present a high microsatellite arrangement complexity like in Rajidae. This low complexity could be due to the short length of the NTS region in Rhizoprionodon where the amount of microsatellites needed to maintain the arrangement and other functions is lower (directly proportional). Therefore it is possible to propose that fewer microsatellites are needed to maintain the arrangement and other functions, but this is a speculative claim.

Figure 5. Linear regression analysis for the length of the NTS and the complexity of the microsatellites present in the NTS region. The analysis shows a high coefficient determination, which was significant. Solid line represents the linear regression according to the equation Y = 10.78 + 0.03(X) + e . Dashed lines represent the 95% confidence interval for the median value of Y_i in any X_i . Thin lines are the prediction limits of 95% for new observations.

Species belonging to Rajiformes have a highly similar NTS which Pasolini et al.⁹ interpret as evidence of the existence of a similar NTS organization in the 5S rRNA genes of the common ancestor of the Rajiformes families from Upper Cretaceous (97 Mya). However, this high level of similarity may be the product of a more recent event of speciation. Valsecchi et al.⁶⁷ by molecular studies of the control region and the 16S rRNA estimated that the divergence of the Rajini tribe could be less than 7 Mya. This estimate is also supported by the paleogeography and paleoclimate events that gave rise to the Mediterranean fauna where the samples of the genus Raja and Dipturus were obtained.^9,68

Conclusions

In this study we show that mutations in the primary structure and the conformation of the secondary structure of the 5S rRNA type II can be tracked in their evolution among the elasmobranchii group. Despite the fact that the secondary structure is highly conserved, there are several mutations present. We advise that future studies take into account the tertiary structure of the 5S rRNA since it might be more informative with respect to the possible changes in the spatial conformation. Moreover, mutations that maintain the secondary structure and others characteristics (e.g. terminal signal for RNA polymerase III and the low variation in the promoters), indicates that this non-coding gene is under strong purifying selection. Also, our study shows that the complexity of the microsatellite present throughout the NTS region is positively correlated with its length, which supports the idea that microsatellite, and therefore the NTS regions, play an important role in the function and evolution of the gene.

Materials and methods

Sample collection and storage

Samples from Urotrygon rogersi and Urotrygon aspidura were collected by artisanal shrimp trawling in Juanchaco and La Bocana localities in the Colombian Pacific Ocean (3 o 52 N − 77 o 18 W ). From each individual approximately 1 cm³ of muscle was extracted and preserved in Longmire's buffer⁶⁹ as suggested by Hleap et al.⁷⁰ until DNA extraction.

DNA extraction, PCR amplification, and DNA sequencing

The DNA extraction was performed with proteinase-K digestion following a phenol/chloroform purification with the conditions reported by Hleap et al.⁷⁰. The PCR amplification of the 5S rDNA type II was performed using the primer set 5SRajid designed to amplify specifically the 5S rDNA type II in the genera Raja and Dipturus.⁹ The primer 5SRajid-F (5′-GGAATACCAGGTGCCGTAGG-3') binds between the positions 55 - 74 of the positive strand, whereas the primer 5SRajid-R (5′-GGTATTCCCAGGCGGTCT-3') binds between the positions 88 – 105 in the negative strand. The PCR reaction contained 1X Taq Buffer, 1 mM Mg, 200 pmol dNTP's, 10 μM of each primer and 1 U of Taq polymerase. Thermal conditions were the following: an initial denaturation step at 95 ° C for 3 min, 30 cycles of 95 ° C for 30 s, 56.5 ° C for 30 s and 72 o C for 1 min. Finally an extension at 72 o C for 7 min was made. A negative control (all above except DNA, the volume of which was replaced by water) was used to test for any contamination. The PCR products were assessed by electrophoresis in 1.0% agarose gels and visualized by ethidium bromide staining under ultraviolet illumination. The PCR purification was done with the polyethylene/glycol (PEG) protocol.⁷¹ Nine sequences were gathered (5 samples of U. aspidura and 4 of U. rogersi) through the standard Sanger sequencing service of Macrogen Inc.. in Korea and were submitted to GenBank (www.ncbi.nlh.nih.gov) under the accession numbers JQ670255 – JQ670263.

In silico data collection, sequence alignment and RNA secondary structure prediction

A total of 82 sequences are reported for the 5S rDNA type II among 13 species of elasmobranchs and 3 Petromyzontiformes were obtained from GenBank (Table 1). Sequences from Galeocerdo cuvier, Scyliorhinus canicula and Alopias superciliosus were assigned to belong to the type II 5S rRNA by BLAST (nucleotide blast).⁷² The BLAST search was optimized for highly similar sequences (megablast) showing higher similarity with the 5S rRNA type II among all elasmobranchii to the exclusion of type I. Similarly, the 5S rDNA sequences of Raja montagui was assigned to belong to type II since it showed higher similarity to the type II sequence of Raja and Dipturus genus instead of type I. The authors⁶⁸ also reported that this sequence corresponded to the heaviest band of 2 bands obtained. Therefore, the study included 91 specimens (82 retrieve from NCBI database, plus 5 of U. aspidura and 4 of U. rogersi described in this study). All sequences were aligned with MUSCLE⁷³ implemented in MEGA V5.0⁷⁴ with default options. In order to ensure the alignment of homologous positions, low complexity regions were manually removed. However, the microsatellite often present in such regions were taken into account to discuss the NTS sequence.

The level of conservation of each nucleotide position of 5S rRNA was assessed with WebLogo⁷⁵ which uses information gained at each position.

The consensus secondary structure of the 5S rRNA type II was obtained through RNAalifold⁷⁶ using the matrix of sequences. However, since a consensus structure can overshadow the structure variability of the molecule, we also made the secondary structure prediction for each haplotype sequence. The haplotype sequences were retrieved with a python script that only keeps unique sequences. Then we use the covariation among positions with the mutual information (M_ij) using:(Eq. 1) ${\displaystyle M_{ij}={\displaystyle {\sum }_{k\in B}{f}_{ij}(XY)lo{g}_2\frac{f_{ij}(XY)}{f_i(X)f_j(Y)}}.}$(Eq. 1) where f_i(X) is the frequency of base X at aligned position i and f_ij(XY) is the frequency of finding X in i and Y in j and (XY) is the set of allowed basepair B = {GC,CG,AU,UA,GU,UG} so that (XY) ∈ B. Mutual information is stringent respect to the information contained in consistent mutations such as GC↔GU. When one position is conserved the result of M_ij is equal to zero. Therefore, we use also the covariation-like measure proposed by Hofacker et al.⁷⁷ divided by 2 since C_ij make use of a symmetric matrix:(Eq. 2) ${\displaystyle C_{ij}=\frac{{\sum }_{XY;X\text{'}Y\text{'}}{f}_{ij}(XY)D(XY;X\text{'}Y\text{'})f_{ij}(X\text{'}Y\text{'})}{2}.}$(Eq. 2) where D(XY;X'Y') is the Hamming distance of (XY) and (X'Y') and both (XY) and (X'Y') ∈ B⁷⁷. Null models were fitted to test for the statistical significance of the results of both M_ij and C_ij . The null models were computed by 1) making 2 vectors with a length equal to the number of sequences studied, where each vector has a composition of 25% for each possible base, 2) estimating M_ij or C_ij between these 2 vectors and storing the results, 3) repeating steps 1 and 2 for 1.0 × 10⁵ iterations, 4) making a cumulative distribution function with the results and finding the value X_ijc such as the result of the integral from X_ijc to X_ij max is equal to 0.01, where X_ij represent both M_ij or C_ij and X_ij max represents the maximum possible value of M_ij or C_ij (2 and 0.75, respectively), and finally 5) using the pair of columns from the real dataset whose values of M_ij or C_ij were higher than its respective X_ijc to constrain the secondary structure prediction of each haplotype in RNAalifold. To define the groups of structures obtained, we built a matrix using the dot-bracket notation of each structure and then we applied the distance analysis of Neighbor Joining with 1000 replicates in SeaView v4.4.⁷⁸

Analysis of the molecular structure of 5S rDNA type II

Analysis was performed by direct comparison of gene sequences between species, focused in the 5S rRNA, internal control region (ICR), microsatellites, TATA box sequences and termination signal of the RNA polymerase III. The microsatellite distribution along the NTS region was analyzed with the Clark-Evans test⁷⁹ and the Diggle-Cressie-Loosmore-Ford (DCLF) test⁸⁰ with 100000 iterations of a Monte Carlo simulation implemented in R v3.0.1 with the module Spatstat.⁸¹ The spatial reference for each microsatellite was its midpoint. In order to prove that a relationship between the complexity of microsatellite regions and the length of the NTS exist, we performed a linear regression analysis. For those samples of Urotrygon only the sequenced regions were taken into account. Also, a chi-square test with Yates's correction was performed to prove statistically that the stop signal in 5S rDNA is equal in all the individuals (Table 2).

Author's contributions

JSH and HCH coordinated the study. JSH collected the samples used in this study. SIC completed the laboratory work, the in silico analysis, the results interpretation and wrote the manuscript. JSH, HCH and CB helped in the analysis, the results interpretation and revised the manuscript. All authors read and approved the final manuscript.

Disclosure of potential conflicts of interest

No potential conflicts of interest were disclosed.

Acknowledgments

We thank Hernán Paredes and Florentino Cuero and sons for the support in fieldwork and Agustin Martínez, owner of the motorboat “Arraiján” to furnish his ship for some of the sampling. To Jitka M. Krejci, Conor Mehan, and Johanna Zaglauer for the review of the English version, as well as the critical review of this manuscript. The SQUALUS Foundation logistically supported this work through contributions of PADI Foundation and IEA scholarships for the project “Biodiversidad e historia de vida de rayas del Pacífico vallecaucano.” The Universidad del Valle supported this work through research grant CI #7809 for the project entitled “Efecto de la selección sobre la heredabilidad estimada en algunos caracteres morfométricos de historia de vida en dos especies del género Urotrygon del Pacífico vallecaucano” and the research grant CI #7922 for the project entitled “Determinación de la relación génica entre el gen 5S rDNA tipo I y tipo II entre peces óseos y cartilaginosos.” NSERC supported this work through the grant No. 120504858. The government of Canada partially sponsored this project through the DFATD scholarship.