Free Publication in 2019
Alberto Arias-Pérez 1, *, Daniel Ramírez-Torres 1, María E. Rodríguez 1, Silvia Portela-Bens 1, Emilio García-Suarez 1, Manuel A. Merlo 1, Aglaya García-Angulo 1, Ismael Cross 1, Thomas Liehr 2, Laureana Rebordinos 1
Received: July 31, 2018 | Accepted: September 27, 2018 | Published: October 24, 2018
OBM Genetics 2018, Volume 2, Issue 4 doi:10.21926/obm.genet.1804044
Academic Editors: Thomas Liehr
Special Issue: Applications of Fluorescence in Situ Hybridization
Recommended citation: Arias-Pérez A, Ramírez-Torres D, E. Rodríguez M, Portela-Bens S, García-Suarez E, Merlo MA, García-Angulo A, Cross I, Liehr T, Rebordinos L. In Silico Detection and FISH Analysis to Determine Location of miRNAs in Solea senegalensis Chromosomes Using BACs. OBM Genetics 2018;2(4):044; doi:10.21926/obm.genet.1804044.
© 2018 by the authors. This is an open access article distributed under the conditions of the Creative Commons by Attribution License, which permits unrestricted use, distribution, and reproduction in any medium or format, provided the original work is correctly cited.
MicroRNAs (miRNAs) are small, non-coding RNAs with lengths ranging from 18 to 25 nucleotides (nt). These RNAs play a very important role in gene expression by regulating mRNA cleavage and translation. The biogenesis of miRNA is a complex process and several pathways are known [1,2]. As a general rule, miRNA synthesis begins in the nucleus where pri-miRNAs (capped and polyadenylated primary nascent transcripts) are transcribed by RNA polymerase II. Then, these pri-miRNAs are processed to shorter precursor molecules (approximately 70 nucleotides) known as pre-miRNAs. This step is carried out by Drosha, an RNAse type III enzyme, along with several other polypeptides. Later the pre-miRNA is transported to the cytoplasm where Dicer and associated proteins generate a duplex that, on average, is 22 nt long. Usually one of the strands will form part of an active miRNA-induced silencing complex (miRISC) [1,2,3].
MicroRNAs mainly influence gene expression by post-transcriptional repression of mRNA in cytoplasm, but some of them are imported back into the nucleus where they regulate the expression of non-coding RNAs . The regulation of mRNA cleavage and translation is related to base complementarity between the miRNA and its target. A perfect pairing induces target degradation while imperfect pairing inhibits translation elongation . Nucleotides at positions 2-8 from the 5’ end of a mature miRNA play a fundamental role in miRNA-to-target mRNA pairing. This “seed” sequence is usually highly conserved among taxa, and is used to classify miRNAs into families .
These non-coding RNAs were initially discovered in Caenorhabditis elegans [7,8] and it took some years to uncover the importance that these non-coding RNAs have in gene expression. Currently, there is extensive ongoing research in both animals and plants in relation to the role of miRNAs in different biological processes including skeletal muscle development , cancer , and responses to stress . In teleosts, the number of studies is still quite limited (but see , , and ). Currently, the number of mature miRNA sequences available at miRBase (v. 22) is 3687, corresponding to only 16 species. It is clear that considerable effort is needed to comprehend how miRNAs generally affect gene expression in fish.
The Senegalese sole, Solea senegalensis (Kaup 1858), is a flatfish species belonging to the Pleuronectiformes order. It is distributed in the Atlantic Ocean from the Gulf of Biscay to the Northwest coast of Africa, and it is less frequently found in the Western Mediterranean . It is a predominantly littoral and benthonic marine species that shows great potential for marine aquaculture diversification due to its growth rates and flesh quality . Nevertheless, the existence of sexual dysfunction of males reared in captivity, high larval mortality, and diseases have hampered its production . In recent years, an increasing number of studies have assessed different genetic and genomic aspects of S. senegalensis such as transcriptomics , integration of genetic maps , and cytogenomics . There are also various studies assessing gene expression , but to date only one study has explored miRNA sequences in S. senegalensis . This study assessed the role of miRNAs in the thermal plasticity of growth, identifying 320 conserved miRNAs and observed a strong influence of temperature on expression of miRNAs.
The integration of sequence information with data on chromosomal physical location is useful for understanding genetic processes and is essential for comparative analysis, QTL mapping, and positional cloning of genes . The combination of fluorescence in situ hybridization (FISH) and bacterial artificial chromosome (BAC) libraries provides an efficient approach for anchoring genomic and linkage sequence data onto physical chromosomes .
In this study, we employed an in silico approach to identify potential miRNAs and their target genes in an S. senegalensis BAC library. This information was then combined with BAC-FISH data. The results will help improve the understanding of the regulatory roles of these miRNAs in Senegalese Sole and other fish, but will also be useful for understanding their evolution and distribution.
Mature miRNA sequences (35828 entries in 223 species) were downloaded from miRBase (Release 22, March 2018; http://www.mirbase.org). In addition, mature sequences from different fish species were also considered: Culter alburnus , Cynoglossus semilaevis , Cyprinus carpio , Hippoglossus hippoglossus , Paralichthys olivaceus [28,29], Scophthalmus maximus , and Solea senegalensis . To avoid redundant results, non-unique mature miRNA sequences were merged (FASTA Unique sequences, https://www.ncbi.nlm.nih.gov/CBBresearch/Spouge/html_ncbi/html/fasta/uniqueseq.cgi).
A nucleotide database was created containing 2981 sequences from 91 BAC clones of a Solea senegalensis BAC library with construction details described in , using the basic local alignment search tool (NCBI suite BLAST 2.6.0 ). This database was then interrogated using blastn command for potentially conserved miRNAs, employing the following parameters: task “blastn”, percent identity of 0.7, percent query coverage per high-scoring segment pair of 1, an expectation value threshold of 0.01, and low complexity sequence data was filtered using the DUST algorithm. All BLAST results with a maximum of 4 mismatches and no gaps were used for further analysis. Potential mature miRNAs and surrounding flanking regions (100 bp on each side, when possible) were extracted with the BLAST tool blastdbcmd. These sequences were used with MIReNA 2.0 , which filters secondary structures using five criteria to predict pre-miRNAs. The program RNAfold 2.4.6  from the ViennaRNA package was used to generate RNA secondary structure plots. Finally, pre-miRNAs were manually inspected to remove possible false positives.
BAC-FISH data was obtained from previous works [18,23,31] and by employing the FISH techniques described in . All applicable international, national, and/or institutional guidelines for the care and use of animals were followed.
To study the evolutionary relationship among the potential miRNAs detected for S. senegalensis and those of other teleosts, the stem-loop sequences available at miRBase were downloaded. When possible, stem-loop sequences from the references mentioned in the previous section were also used. The MEGA X software  was used to align the sequences (via the MUSCLE algorithm with default parameters) and to build a tree by the maximum likelihood method (by means of partial deletion with site coverage cut-off of 95% with 1000 bootstrap replications).
The reference transcriptome for S. senegalensis (v4.1) was downloaded from Solea DB  and used for target prediction using miRanda v3.3a . Only the last part of the sequences was used, given that the 3'-UTR region for most S. senegalensis genes has not yet been resolved as well as to speed up the computational analysis. All sequences were conservatively trimmed to a maximum of 1000 nt, given that previous data on teleosts showed a median 3'-UTR length of 795 nt . The parameters employed were: minimum score of 50, free energy ≤ 20 kcal/mol, and strict alignment in the seed region.
The advent of NGS technologies and the feasibility of large-scale genomic projects have made sequencing data available on an unprecedented scale. This makes it possible to search different types of sequences and even to analyze the chromosome structure of non-model organisms. In the present study, an in silico approach was used to screen miRNAs, which play a fundamental role in the regulation of gene expression, in BAC data of S. senegalensis. The BLAST analysis of 2981 contigs (with an average length of about 4800 bp) from 91 clones of a BAC library of S. senegalensis yielded 272 matches. These matches, together with some flanking sequences, were extracted for pre-miRNA and secondary structure prediction with MIReNA and RNAfold. This resulted in 124 potential miRNAs; however, given that some of them map to approximately the same position due to differences in the mature sequence length, 45 groups were obtained (Table 1 and Table S1). The largest group contained 16 mature sequences ranging in length from 20 to 23 bp. When considering only one mature sequence per group, that being the one with the lowest e-value, most of them were observed once but others were found several times: for example, miR-467g and miR-8915 appeared nine and five times, respectively. Other miRNAs appeared three times (miR-669l-3p, miR-467e-3p) and twice (let-7f, dme-miR-4968-3p, cse-PC-5722-5p, and let-7f-5p).
Table 1 Potential miRNAs identified in Solea senegalensis BAC sequences (full details in Table S1). Group numbers (column 1) with two integers separated by a dot indicate multiple mature miRNAs matching to the same or almost the same position. When a miRNA name (column 3) corresponds to a group name (column 2), this means that this sequence (column 4) is the 5’ or 3’ arm of the other mature sequence. Sequence names (column 4) contain the following four items of information separated by dashes: chromosome number (99 indicates hybridization to multiple chromosomes, na indicates no hybridization); BAC name (plate number is the first two characters and plate position is the last three characters); sequencing batch; and contig number. Other keys given at foot of table.
|Group||Name||miRNA homolog||Sequence||Positions||Side||NM||Strand||LP (nt)||A U||MFE||MFEI|
Key: NM: number of mismatches; LP: length of the precursor; A U: A U content; MFE: minimum folding free energy (kcal/mol); MFEI: minimal folding free energy index. Species abbreviations: Abu Astatotilapia burtoni, Cca Cyprinus carpio, Dre Danio rerio, Fru Fugu rubripes, Gmo Gadus morhua, Hhi Hippoglossus hippoglossus, Ipu Ictalurus punctatus, Mze Metriaclima zebra, Nbr Neolamprologus brichardi, Ola Oryzias latipes, Oni Oreochromis niloticus, Pny Pundamilia nyererei, Pol Paralichthys olivaceus, Ssa Salmo salar, Sse Solea senegalensis, Tni Tetraodon nigroviridis.
Among the 124 potential miRNAs, 52 have been not reported in fish (Table S1) and the remaining ones have been reported in 1 to 13 species (average 2.75). The species showing the greatest number of records were C. carpio (22), P. olivaceus (22), Danio rerio (20), and S. senegalensis (19); the sequences identified in the greatest number of species were aca-let-7a-5p (13), aca-let-7e-5p (11), cal-let-7-5p-13 (11), cal-miR-153-3p-1 (10), aca-let-7c-5p (10), and aca-miR-205a (10). To date, only one work has studied miRNAs in S. senegalensis , and the information available for the teleost species is limited (but see [3,12,13]). Therefore, the results obtained here will contribute to a better understanding of the evolution of miRNAs in teleosts.
To assign miRNAs to a particular chromosome, BAC-FISH data was employed [18,23,31]. The BAC library was screened for genes of interest and then positive clones were sequenced and annotated. Following this, FISH probes were prepared and hybridization was carried out (Figure 1). Most of the groups (26 of 45) could not be assigned to a particular chromosome, but 8, 3, and 2 groups were assigned to chromosomes 1, 4, and 6, respectively. For other chromosomes (12, 16, 18, 19, and 21) only one group was observed. Finally, one group of miRNAs corresponded to a BAC that hybridizes with multiple chromosomes. Although chromosome 1 seems to have a higher amount of miRNAs, the amount of data available does not allow significant differences to be inferred. This is the case in D. rerio, where miRNA genes have shown unequal distribution among chromosomes . The analysis of the distribution of miRNAs in the Senegalese sole chromosomes might help in studying the biological implications of the genomic context of miRNA. In addition, the integration of QTL and miRNA could be used to investigate biological aspects such as disease susceptibility .
Figure 1 FISH of BACs 46P22 (green), 03C15 (pink), 08A23 (red) and 09J04 (blue).
The length of mature miRNAs ranged between 18 and 25 nt with a mode of 22; a length of 21 and 22 nt accounted for 52% of the 124 matches. The average length of the predicted precursor was 80 nt with a minimum and a maximum length of 60 and 125 nt, respectively. The secondary structures computed from these sequences showed the typical stem-loop miRNA structure (Supplementary data 1), except for sse-PC-5722-2 and sse-PC-5722-2 (groups 8 and 22) that had a more cruciform shape. This putative miRNA was described in C. semilaevis  by the sequencing of small RNA libraries. The precursor found in this study and in that of C. semilaevis differed markedly in size (112 vs 61 nt), and the alignment between these sequences was poor. This could be due to two distinct genomic loci expressing identical mature sequences or to the precursor found here belonging to another RNA class.
The placement of the mature miRNA sequences on the 5' or 3' arm of the stem-loop structure was almost equally distributed among the groups of miRNAs (23 and 22 results for the 3' and 5' arms, respectively), and the members of all groups were placed on the same side. Although usually only one of the strands of the miRNA duplex plays a role in the processes of gene regulation, there is evidence that some miRNAs can provide two mature functional miRNAs . Our results revealed that there were several instances where significant matches for both strands of the precursor were observed: (a) 6 (cca-miR-460-5p) and 7 (sma-miR-460-3p*), (b) 10 (hhi-miR-153b) and 11 (dre-miR-153b-5p), (c) 12 (cse-PC-15046-3p) and 13 (cse-PC-2293-5p), (d) 25 (aca-let-7e-5p) and 26 (aca-let-7e-3p), (e) 27 (prd-let-7-5p) and 28 (pol-let-7a-3-p3), (f) 29 (efu-miR-205) and 30 (ola-miR-205), (g) 38 (cse-miR-454-3p) and 39 (tgu-miR-454-5p), (h) 40 (ola-miR-301a-3p) and 41 (cal-miR-301b-5p-1), and (i) 44 (cal-miR-130a-5p-2) and 45 (aca-miR-130b-3p).
Several miRNAs were located relatively close to one another: (a) sse-miR-8915-2 and sse-miR-8915-3; (b) sse-miR-8915-4 and sse-miR-8915-5; (c) sse-let-7e and sse-let-7a; and (d) sse-miR-454, sse-miR-301b, and sse-miR-130c. This suggests that they might belong to a miRNA gene cluster as has been observed for many miRNAs , which commonly form part of polycistronic transcripts. It is believed that these gene clusters could help coordinate the regulation of biological processes. The identification of the promotor regions of these miRNAs could help determine if these sequences form a true cluster.
In regard to the putative pre-miRNA, the minimum folding free energy for secondary structures ranged from -69.3 to -22.4 kcal/mol (with an average of -34.79), and the adenosine and uridine content (A U) ranged from 16% to 90% (with an average of 58%). These average values are similar to values observed in other fish [25,39,40]. Again, precursors showing similarity with the mature sequence cse-PC-5722-5p stand out because of their low A U percentage (0.16).
The minimum folding free energy value is an important parameter for evaluating the potential secondary structure of RNAs. miRNAs have lower values than some RNAs, but not others. Therefore, additional criteria are required, such as the minimum folding free energy index. This parameter ranged from 0.70 to 3.70 (with an average of 1.09) in our results. Studies in different fish have shown similar or lower values, such as Huang et al.'s study in eleven fish species (range: 0.58 to 0.91, average: 0.71) ; C. semilaevis: (range: 0.51 to 1.31, average: 0.90) ; and in Boleophthalmus pectinirostris (range: 0.51 to 1.25, average: 0.90) . According to the results obtained by Zhang et al. 's , a majority of miRNA precursors have values over 0.85, allowing miRNAs to be discriminated from other coding and non-coding RNA. In our case, 80% of the results had values higher than 0.85, indicating that the predicted miRNAs are true positives and not other types of RNA.
The precursors of several miRNAs (let-7, miR-130, miR-153, miR-205, miR-301, miR-454, and miR-460) were compared with available stem-loop sequences of the same family from other fish (Supplementary data 2 and Figure 2).
Figure 2 Maximum likelihood phylogenetic tree for precursor sequence miR-130 from Solea senegalensis (marked with a black dot) and other fish. Species abbreviations: Abu: Astatotilapia burtoni, Cca (sequences from ) and Ccr: Cyprinus carpio, Dre: Danio rerio, Eel: Electrophorus electricus, Fru: Fugu rubripes, Gmo: Gadus morhua, Hhi: Hippoglossus hippoglossus, Ipu: Ictalurus punctatus, Mze: Metriaclima zebra, Nbr: Neolamprologus brichardi, Oni: Oreochromis niloticus, Ola: Oryzias latipes, Pol: Paralichthys olivaceus, Pny: Pundamilia nyererei, Ssa: Salmo salar, Tni: Tetraodon nigroviridis.
The miRNA family with the largest number of sequences was let-7 with 215; this includes sequences for closely-related mature sequences (from let-7a to let-7j). Although most bootstrap values were low (below 70), sequences usually formed clusters composed of sequences of the only same type (e.g. let-7d, let-7f, let-7g, and let-7i). This points to a high degree of evolutionary conservation. Among the sequences showing similarity to let-7 in S. senegalensis, the sequence sse-let-7a formed a group with let-7a miRNAs, including one Pleuronectiform (C. semilaevis, cse-let-7a-5p). The sequence sse-let-7e was part of a cluster with good bootstrap support (88) containing all let-7e sequences, including cse-let-7e-5p from C. semilaevis, plus two sequences from C. alburnus (cal-let-7-5p-22 and cal-let-7-5p-11).
The tree for mir-130 (Fig. 2) did not show a clear clustering among related sequences (mir-130a, mir-130b, and mir-130c). Sequence sse-miR-130a formed a group displaying good bootstrap support (99) with a combination of mir-130a and mir-130b sequences. On the other hand, sse-miR-130c clustered with low bootstrap support (63) with several Cichliformes and one Cypriniform (mir-130 and mir-130a sequences).
Among the mir-153 sequences available for fish, there were three types of closely-related mature sequences: mir-153a, mir-153b, and mir-153c. Most mir-153a sequences, apart from ssa-mir-153a-1 and ssa-mir-153a-2, formed a cluster (with a bootstrap value of 81). For mir-153b and mir-153c there was not a clear separation, and sequence sse-miR-153b clustered with 153b and 153c miRNAs.
Two main clusters (with a bootstrap value of 70) were observed for mir-301, one containing mostly mir-301a sequences and the other mostly mir-301b and mir-301c sequences. The sequence sse-miR-301b was positioned close to mir-301b sequences, including cse-miR-301b-5p from C. semilaevis.
No clear clustering was observed for sse-miR-460 with miR-460 sequences, despite having a sequence for C. semilaevis, a Pleuronectiformes species. Nevertheless, well-supported clusters were observed for Cichliformes and Cypriniformes species (as well as a Siluriformes species). Two other miRNAs, miR-205 and miR-454, did not show a clear clustering.
In addition, the precursor sequences of two miRNAs, miR-467 and miR-8915, appeared multiple times and were aligned. Three of the miR-467 precursor sequences were identical (sse-miR-467g-1, sse-miR-467g-2, and sse-miR-467g-4), and the others ranged in length from 64 to 105 nt. When the maximum likelihood tree was obtained (data not shown), two groups were observed (with a bootstrap value of 100), one with the identical sequences and the other with the remaining ones. The precursor sequences for sse-miR-8915 were almost identical, with one sequence being 1 nt shorter (sse-miR-8915-1) and the other showing a substitution (sse-miR-8915-4, position 15, with an A instead of a G). This strongly suggests that miR-467 and miR-8915 genes might have suffered gene duplications during evolution, and that there are several loci producing the same mature sequence.
The observation that some miRNAs are considerably conserved among teleosts suggests that the information obtained in S. senegalensis might be transferable to other fish species.
Of the 59514 sequences present in the reference transcriptome for S. senegalensis, 12387 were putative targets according to miRanda, with some of them appearing up to seven times. Among the 34 unique sequences with similarity to mature miRNAs that showed potential targets, sse-miR-3584 and sse-miR-4698 produced the most results (with more than twice as many potential targets as the next miRNA), with 3234 and 2748 hits, respectively.
The top ten most frequent gene ontology terms (out of a total of 5857) are presented in Table 2. The predicted targets are implicated in a great variety of processes such as ion binding and nucleic acid binding, or are membrane-related. The sequences showing the highest score for each of the 34 mature sequences are shown in table S2. For example, sse-PC-15046 (solea_v4.1_unigene226900) targets a gene in D. rerio  that encodes an E3 ubiquitin ligase that targets a receptor in the Notch signalling pathway, which is involved in a wide range of developmental processes. The sse-miR-153b mature sequence targets a gene of the odorant receptor gene family; these genes exhibit several transmembrane domains and show great diversity in teleost fish . The number of predicted target genes observed implies that miRNAs in this species modulate a wide range of processes. This information is essential for understanding development and other physiological processes, such as larval and sex development as well as disease resistance, and it will contribute to improving the aquaculture process in this and other species. Nevertheless, these results need to be taken with caution. Several factors might have a considerable impact on the rate of false positive results: (a) the lack of information on gene structure without properly defined 3' UTR regions, as is the case for most S. senegalensis genes, and (b) the difficulty of predicting and investigating the precise way that a miRNA and its mRNA target pair interact . Experimental validation is required to confirm the predictions obtained here.
Table 2 Ten most frequent gene ontology terms found for the targets of 34 putative miRNAs in Solea senegalensis.
|zinc ion binding
|integral component of membrane
|nucleic acid binding
|DNA-binding transcription factor activity
|calcium ion binding
In this work, an in silico approach was used to identify miRNAs in S. senegalensis from sequences of a BAC library. Using the mature sequences available at miRBase and other sources, 45 groups of sequences showed similarity to mature miRNAs. Of these, 19 were assigned to Senegalese Sole chromosomes using BAC-FISH data. Phylogenetic trees showed different levels of support for clustering among types of miRNAs, and S. senegalensis sequences usually were close to other Pleuronectiformes. Thirty-four unique mature sequences targeted 12387 sequences from an S. senegalensis reference transcriptome. This suggests that miRNAs might modulate a wide range of processes in this species. These results might be helpful for understanding processes such as sexual differentiation, larval development, and disease susceptibility, which hamper the commercial utilization of this species in aquaculture. The phylogenetic conservation of some miRNA implies that the information obtained in the Senegalese sole might be transferable to other aquaculture fish species. Nevertheless, experimental validation by RT-PCR, northern blotting, or microarrays is required.
The following additional materials are uploaded at the page of this paper.
1. Table S1. Potential miRNAs identified in Solea senegalensis BAC sequences. Group numbers (col. 1) with two integers separated by a dot indicate multiple mature miRNAs matching to the same or almost the same position. When a miRNA name (col. 3) corresponds to a group name (col. 2) this means that this sequence (col. 4) is the 5’ or 3’ arm of the other mature sequence. Sequence names (col. 4) contain the following four items of information separated by dashes: chromosome number (99: hybridization to multiple chromosomes, na: no hybridization); BAC name (plate number: first two characters, and plate position: last three characters); sequencing batch; and contig number.
2. Table S2. Sequences from Solea senegalensis reference transcriptome (v4.1) showing the highest score as targets of 34 putative miRNAs.
3. Supplementary data 1. Secondary structures for the 124 miRNA precursors found for Solea senegalensis. For miRNA group identification numbers see table S1.
4. Supplementary data 2. Maximum likelihood phylogenetic trees in Newick format for precursor sequences let-7, miR-130, miR-153, miR-250, miR-301, miR-454 and miR-460, from Solea senegalensis and other fish. Species abbreviations: Abu: Astatotilapia burtoni, Cca (sequences from Huang et al. Gene Reports, 2017, 8: 30-36) and Ccr: Cyprinus carpio, Dre: Danio rerio, Eel: Electrophorus electricus, Fru: Fugu rubripes, Gmo: Gadus morhua, Hhi: Hippoglossus hippoglossus, Ipu: Ictalurus punctatus, Mze: Metriaclima zebra, Nbr: Neolamprologus brichardi, Oni: Oreochromis niloticus, Ola: Oryzias latipes, Pol: Paralichthys olivaceus, Pny: Pundamilia nyererei, Ssa: Salmo salar, Tni: Tetraodon nigroviridis.
5. Supplementary data 3. Target prediction results from Solea senegalensis reference transcriptome (v4.1) and 34 putative miRNAs using miRanda. Unigen_id: unigene identification code, Functional_description: functional description (see Solea-DB for details), GO_terms: gene ontology terms, mature_miRNA: name of the mature miRNA, Tot Score: total score, Tot Energy: total energy, Max Score: maximum score, Max Energy: maximum energy, Len1: mature miRNA sequence length, Len2: length of the reference sequence examined, Positions: potential target sites.
AAP, DRT and LR conceived and designed the analysis and AAP and DRT developed the pipeline. LR, IC and TL provided guidance and support. MER, SPB, EGS, MAM and AGA contributed to BAC-FISH data collection and analysis. AAP, DRT and LR wrote the manuscript. All the authors have read and contributed to the manuscript.
This study has been supported by the Spanish Ministerio de Ciencia e Innovación MICINN-FEDER (project AGL2014-51860-C2-1-P awarded to Laureana Rebordinos).
The authors have declared that no competing interests exist.