Functional evolutionary convergence of long noncoding RNAs involved in embryonic development

1 This is the Author Accepted Manuscript file of the article: Olazagoitia-Garmendia A*, Senovilla-Ganzo R*, García-Moreno F, Castellanos-Rubio A. (2023) Functional evolutionary convergence of long noncoding RNAs involved in embryonic development. Commun Biol. 5;6(1):908. doi: 10.1038/s42003-023-05278-z. 2 Functional evolutionary convergence of long noncoding RNAs involved in embryonic development. Ane Olazagoitia-Garmendia1,2*, Rodrigo Senovilla-Ganzo3*, Fernando Garcia-Moreno1,3,4, Ainara Castellanos-Rubio1,2,4,5,& 1 University of the Basque Country, UPV-EHU, Leioa, Spain 2 Biobizkaia Health Research Institute, Barakaldo, Spain 3 Achucarro, Basque Center for Neuroscience, Leioa, Spain 4 Ikerbasque, Basque Foundation for Science, Bilbao, Spain 5 CIBERDEM, Madrid, Spain *These authors contributed equally to this work &Corresponding author: Ainara Castellanos Rubio, PhD, Barrio Sarriena, s/n, Maria Goyri Building, Lab 1.11, UPV/EHU, Leioa, 48940, Bizkaia, Spain, Phone: +34946007762; Fax: +34946013145, Email: [email protected]; ORCID: 0000-0003-3192-9544 Abstract Long noncoding RNAs have been identified in most vertebrates, but the functional characterization of these molecules is challenging, mainly due to the lack of linear sequence homology between species. In this work, we aimed to find functional evolutionary convergent lncRNAs involved in development by screening of k-mer content (non linear similarity) and secondary structure-based approaches combining in silico, in vitro and in vivo validation analysis. From the Madagascar gecko genes, we have found a non-orthologous lncRNA with a similar kmer content and structurally concordant with the human lncRNA EVX1AS. Analysis of functionrelated characteristics together with locus-specific targeting of human and gecko EVX1AS (i.e. CRISPR Display) in human neuroepithelial cells and chicken mesencephalon have confirmed that gecko EVX1AS-like lncRNA mimics human EVX1AS function and induces EVX1 expression independently of the target species. Our data shows functional convergence of non-homologous lncRNAs and presents a useful approach for the definition and manipulation of lncRNA function within different model organisms. 3 Introduction Long noncoding RNAs (lncRNAs) are RNA molecules longer than 200bp in length that do not have coding potential1. LncRNAs are gaining importance due to their involvement in a wide range of biological processes, and some of them have been described to be implicated in different aspects of embryonic development2. However, the study of lncRNA relevance through their evolutionary conservation has been challenging due to their lack of linear sequence homology among species3–5. Evolutionary conservation is widely used as an indicator of the functional significance of newly discovered genes, and the simple search for homology at the nucleotide level has proven to be valuable for protein-coding genes. However, lncRNAs with similar functions often lack linear sequence homology which implies lncRNA function cannot be readily assigned from their nucleotide sequence. K-mer based comparison methods have been demonstrated to be useful to find functionally related lncRNAs with different spatial arrangements of related sequence motifs, where a k-mer is defined as all possible combinations of a continuous sequence of nucleotides of a given length k. K-mer-based classification has been demonstrated to be a powerful approach to detect recurrent relationships between motif sequence and function in lncRNAs even in the lack of evolutionary conservation6,7. Additionally, lncRNAs with the same secondary or tertiary structure can exert identical molecular functions despite divergent nucleotide sequences8, thus analyzing structural equivalence could also help identify lncRNAs with evolutionary preserved mechanisms9. Moreover, the molecular role of lncRNAs is tied to other characteristics such as the subcellular localization, the abundance within the cell or the interactions with other molecules10,11. In this work, we have taken advantage of k-mer and structure similarity analyses to find evolutionary convergent lncRNAs involved in development. For this purpose, we used Madagascar ground gecko (Paroedura pictus) as the target species to find functionally conserved lncRNAs related to embryonic development. The Madagascar gecko is a useful species to investigate the evolutionary path of vertebrate features due to its phylogenetic position within the squamates order of reptiles12,13. Additionally, we have performed several in vitro and in vivo analyses to investigate the functional cross-species conservation of the candidate lncRNAs. For the in vivo analyses, we have used chicken (Gallus gallus) embryos as an extra-phyletic species to both human and gecko species. The chicken is the most suitable animal model for experimental embryology, and it is a perfect playground for genetic engineering during development14–16. Using these approaches, we found that human EVX1AS and gecko EVX1AS-like regulate the coding gene EVX1 independently of the recipient species molecular machinery. This evolutionary 4 functional convergence of the two non-syntenic lncRNAs emerged independently in the two species, evolving in parallel to play equivalent functions from non-orthologous sequences. Results Nonlinear sequence similarity between human and gecko embryogenesis lncRNAs We hypothesized that functionally related lncRNAs involved in embryonic development could harbor related motif contents, although lacking linear sequence similarity. To test this, we used the SEEKR (sequence evaluation from k-mer representation) standalone6 to find development associated lncRNAs (with already described functions in humans) that could be functionally convergent between human and gecko. Using the All human lncRNA (Gencode v41) set as a normalization set, we calculated the k-mer profile from k=3 to k=6 k-mer lengths of human NEAT1, MEG3 and EVX1AS against all sequences from P. picta genome v1 in order to find equivalent lncRNAs17. For MEG3 and EVX1AS, candidate noncoding transcripts were selected based on their noncoding nature, according to the Coding Potential Calculator18, and its high score in all kmer analyses (Pp-MEG3-like and Pp-EVX1AS-like from now on). Both candidates were present in the highest percentiles (>99 percentile) in all k-3 to k-6 k-mer analyses (Figure 1a). However, no such candidate was found for NEAT1. To verify the validity of this approach, we also compared the k-mer content distribution between human and mouse lncRNAs as they are known to be functionally conserved. These comparisons resulted in comparable percentile values to those found between human and gecko lncRNAs with 100 percentile for MEG3 and >97 percentile for EVX1AS. To evaluate the specificity of our results, we also performed the reverse analysis and compared our gecko candidates with the human set of all human lncRNAs. The results of this reverse SEEKR analysis, showed that our lncRNA candidates maintain a high k-mer correlation and that they are in the highest percentile (>96) in all k-mer analyses (from k=3 to k=6) (Supplementary Figure 1a). We also confirmed the balanced GC-content of our transcripts19 (58% and 67% for human MEG3 and gecko MEG3-like, 57% and 52% for human EVX1AS and gecko EVX1AS-like). All together, these results suggested these selected gecko noncoding RNAs are the most likely functional equivalents to human lncRNAs and ruled out the possibility of the candidates having been selected by chance. Given that the gecko noncoding transcripts were 6-to-8 times longer than their human counterparts (4448 bp vs 554bp for EVX1AS and 3751 bp vs 632bp for MEG3), we used 5 CROSSalign9 to identify regions of structural similarity between the different length profiles (Figure 1b). Both lncRNA pairs presented a structural distance lower than 0.095 (where 0 means identical structural profiles) and a 90% structure correlation. The normalized structural distance between the secondary structure profiles of EVX1AS was calculated as 0.087 (p-value=0.001) with a correlation of 90%. In the case of MEG3 the distance was 0.09 (p-value=0.01) with a correlation of 89%. Dinucleotide shuffled sequences of the gecko lncRNAs showed structural distances higher than 0.09 in both cases (0.097 for EVX1AS and 0.092 for MEG3). Interestingly, we did not find any specific regions of conventional linear sequence homology between the pairs of structurally equivalent lncRNAs by mVISTA using the 70% conservation over a 100bp window criteria, neither by a less stringent DotPlot analysis (with a 20% similarity) (Figure 1c and Supplementary Figure 1b)20,21. Thus, our analysis revealed the existence of a pair of lncRNAs with high levels of k-mer similarity and significant structural equivalency between two distant vertebrate species. Together, it suggested that these lncRNAs could play similar functions in embryonic development despite a lack of linear sequence similarity. Human and gecko EVX1AS are principally expressed in brain In order to verify if the gecko lncRNA candidates were actually expressed, we quantified the expression of the Pp-MEG3-like and Pp-EVX1AS-like in different tissues from pre-hatching geckos. Pp-MEG3-like was only expressed in the tail and the carcase (Figure 2a), with a very similar expression level in the two tissues. Expression analysis of Pp-EVX1AS-like demonstrated that this lncRNA is widely and tissue-specifically expressed in this species, showing the highest expression in the brain, 9 to 280 times higher than in the rest of the tissues (Figure 2a, Supplementary Figure 2a). LncRNA Evx1as has been described to transcriptionally regulate its nearby coding gene Evx122. Analysis of gecko Evx1 expression in the embryonic tissues showed the wide characteristic expression of this coding gene, with the highest expression levels in the heart and the carcase and the lowest in the lung (Figure 2a, Supplementary Figure 2a). We also analyzed the expression of both EVX1AS and EVX1 in an RNA pool of different human tissues purchased from Clontech. As previously observed in gecko embryos, EVX1AS was widely expressed among the different tissues with the highest expression also present in the brain (with 1.5-12 times higher values) and the lowest in the colon and the thymus (Figure 2b, Supplementary 6 Figure 2b). Regarding the coding gene EVX1, it likewise presented a broad tissue expression with the highest levels in the kidney and the lowest also in the lung (Figure 2b, Supplementary Figure 2b). Given that human brain tissue is not an accessible material, we also evaluated the expression of EVX1AS and EVX1 in the human neuroepithelial cell line SHSY5Y. Both genes were expressed in this cell line confirming the validity of this culture model for further EVX1AS analysis (Figure 2c). Given the multi-level similarities, from tissue expression and structure to motif content, we selected Pp-EVX1AS-like lncRNA for further analysis of its potentially shared function. Human and gecko EVX1AS are likely to have evolved from independent ancestor sequences. In contrast to coding genes, it has been described that there are ortholog lncRNAs with limited linear sequence similarity due to molecular divergence through evolution. However, these ortholog lncRNAs share common synteny and some k-mers or regions are conserved to retain its functionality23. Hence, despite no significant linear sequence similarities could be detected between human EVX1AS and gecko EVX1AS-like (Figure 1c, Supplementary Figure 1b), we further evaluated the potential evolutionary relationship of our gecko lncRNA. To this aim, we first used UCSC genome browser to visualize the conservation in the vicinity of the EVX1 gene and to analyze the evolutionary origin and dynamics of EVX1AS transcript. As observed in Supplementary Figure 3a, while sequence conservation of the EVX1 gene is conserved in the different species, EVX1AS alignment starts decreasing in the mouse and is totally lost from the lizard on. Moreover, we also compared the colocalization of genetic loci (“synteny”) between human and gecko. Initial analysis of the genomic neighborhood (Figure 3a, Supplementary Figure 3b) of both lncRNAs in each species revealed great differences in synteny of both transcripts. However, to further validate the appearing independent origin of our lncRNA, we employed Satsuma2 and GENESPACE pipeline24, which includes orthofinder25 and MCScanX26. As P. picta proteome is based on genome-annotation translation (Transdecoder, Haas, BJ https://github.com/TransDecoder/TransDecoder), Satsuma227 allowed us to compare large and complex DNA sequences (whole-genome comparison). Synteny results showed no cross- 7 correlation between human EVX1AS chromosome 7 and gecko EVX1AS-like genomic loci (Figure 3b). Nonetheless, we carried out GENESPACE pipeline to assess the existence of orthologous gecko EVX1 orthogroup loci in gecko. Both orthofinder and MCScanX define pp47247 as the orthologous of human EVX1, which is in scaffold 11. This is a distinct chromosome from ppEVX1AS-like and does present cross-correlation with the EVX1 region in human chromosome 7 by both Satsuma2 and MCScanX (Figure 3b). Among the derived transcripts of pp47247, we have not detected any long non-coding transcript (CPC2) that could reproduce EVX1/EVX1AS model. However, it would be interesting to obtain further transcriptomic evidence by deep noncoding RNA sequencing. As Paroedura picta’s genome is not based on long read data, we wanted to further confirm the genomic sequence of EVX1AS as well as its synteny. We amplified and re-sequenced our gecko EVX1AS-like candidate and confirmed that the gecko locus corresponding to our candidate lncRNA is correctly assembled in the scaffold 15. Thus, gecko EVX1AS-like localizes in scaffold 15, while gecko Evx1 ortholog is in scaffold 11. To sum up, all our synteny analyses indicate that both human EVX1AS and gecko EVX1AS-like evolved separately, from independent ancestor sequences, but converged into similar characteristics in parallel. Thus, we delved into these analogue function-related characteristics. Human and gecko EVX1AS share function-related characteristics RNA structure, subcellular localization and abundance of lncRNAs are generally related to their function and molecular roles. Thus, to experimentally assess their functional convergence, we evaluated several function-related characteristics (structure, localization and abundance) for both human and gecko EVX1AS. These two lncRNA forms showed an almost identical migration pattern when in vitro transcribed lncRNAs corresponding to the structural analogue region were migrated in a non-denaturing agarose gel (Figure 4a). Conversely, an unrelated in vitro transcribed lncRNA, which was used as a control, showed a totally different migration pattern, confirming the equivalency of the secondary structures of human and gecko EVX1AS predicted in silico (Figure 1b). 8 For subcellular localization assessment, the amounts of EVX1AS were quantified in whole and nuclear fractions of human neuroepithelial cells and the gecko brain. As previously described22, EVX1AS was found in both cellular compartments suggesting that apart from regulating EVX1 in cis, this lncRNA also exerts other yet undescribed biological functions. The percentage of EVX1AS present in the nucleus was very similar in both species (approximately 30%) (Figure 4b), further supporting its functional convergence. Evaluation of EVX1AS abundance was performed using human and gecko brain cDNA and a reference plasmid of each lncRNA (Supplementary Figure 4a). We had the limitation that while the gecko brains tested were embryonic, human brain samples were from adult individuals. However, the amount of EVX1AS molecules per cell in both cases was lower than 1 (Figure 4c), pointing to a cell type specific expression of the lncRNA within the brain. The higher amount of copies per cell present in the gecko brain may be suggestive of a high peak of activity of this lncRNA regulating EVX1 during embryonic development. The regulation of mouse Evx1 expression by Evx1as has been described to be mediated by the binding of the lncRNA to the mediator complex22. The mediator complex is a multiprotein complex that functions as a transcriptional coactivator and interacts with a wide range of proteins. MEME motif analysis of our lncRNAs revealed that despite differentially located within their sequence (even in opposite sense) there are several common enriched protein binding motifs within both lncRNAs (Supplementary Figure 4b). We also used lncLOOMv228, which searches for synteny conserved k-mers. However, none of the ordered 6-9 nucleotide motifs in both lncRNA sequences (Supplementary Figure 4c) matched with a TargetScan functional motif, supporting our previous synteny results. Interestingly, alignment of the MEME enriched motifs using TOMTOM Motif comparison tool showed enrichment of several brain function related protein binding motifs, as MATR3, HNRNPA1 and PTBP1 (Supplementary Figure 4d), that have been reported to interact with members of the mediator complex according to GENEMANIA server29. Thus, to analyze if both our lncRNA forms also bind the mediator complex, we performed an RNA immunoprecipitation experiment of MED1 protein using lysates from SHSY5Y cells. Immunoprecipitation of human MED1 was able to retrieve both, endogenous and IVT human EVX1AS. Additionally, we also observed the lncRNAMED1 interaction when IVT Pp-EVX1AS-like was added to the protein lysates. Conversely, our IVT lncRNA control did not interact with MED1, confirming the specificity of the binding (Figure 4d). 9 Equivalence in structure, subcellular localization and abundance of EVX1AS from human and EVX1AS-like from gecko, together with the interaction with the mediator complex, strongly suggested a shared cellular function. Human EVX1AS and gecko EVX1AS-like are functionally convergent To analyze whether the previous results were representative of a functional in-cell equivalence, we conducted overexpression experiments of both lncRNAs in two experimental paradigms: in vitro, in SHSY5Y cells; and in vivo, in chick embryonic brains. For this approach, we used the CRISPR-Display technique30, in order to specifically localize EVX1AS into the promoter of human or chicken EVX1. For in vitro experiments, we cloned either the human or gecko EVX1AS lncRNA preceded by two different guide RNAs specific for human EVX1 promoter (HS1 and HS2 being sgRNAs targeting human EVX1 TSS1 and TSS2 respectively). Otherwise, for in vivo experiments, gRNAs for chicken Evx1 promoter were designed (GG1 and GG2 as sgRNAs targeting chicken Evx1 TSS1 and TSS2 respectively). These constructs were followed by a 3’ box and cloned into a CMV driven vector to guide and tether the EVX1AS RNA to the EVX1 promoter (Supplementary Figure 5a, b). The transfection of Hs-EVX1AS vectors into the SHSY5Y cells showed a 4-13x106 fold statistically significant overexpression (p<0.01 for HS1 and p<0.001 for HS2) (Figure 5a). Using Pp-EVX1AS-like plasmids, 1.5-3x106 fold overexpression could be achieved, being significant only with HS2 gRNA (p<0.001) (Figure 5b). The directed overexpression of the human EVX1AS was able to induce the expression of EVX1 mRNA in the SHSY5Y human neuroblastoma cell line about 2 times (p<0.01 for HS2) (Figure 5a). When gecko EVX1AS-like was overexpressed in the human cells, we also observed a 1.5-2-fold increase in the mRNA levels of EVX1 (p<0.01 for HS2) (Figure 5b), similar to what it was previously observed using this same approach22. No induction of GAPDH negative control could be observed in neither of the cases (Supplementary Figure 5c,d) confirming that these two lncRNAs functionally converged in their role of controlling EVX1 mRNA expression. Interestingly, the tethering mediated by the sgRNA2 (named HS2) (closer to the EVX1 start codon) performed more efficiently, independently of the origin of the lncRNA. To further confirm the functional convergence of human and gecko EVX1AS, we tethered an unrelated human lncRNA to EVX1 promoter using the most efficient sgRNA (HS2). We confirmed that tethering of this lncRNA does not induce the expression of EVX1 (Supplementary Figure 5e,f), supporting that only EVX1AS homologs lead to upregulation of EVX1. 16 extracted as usual. The amount of specific nuclear RNA measured by RT-QPCR was compared to the total amount of the RNA in the whole cell. Quantification of molecules per cell. In order to determine the EVX1AS copy number in human and gecko cells, a reference plasmid incorporating the cDNA sequence of each EVX1AS form was used. Absolute quantification was performed using 5 ten-fold serial dilutions of the reference standard. Ct versus the dilution factor was plotted in a base-10 semi-logarithmic graph, fitting the data to a straight line. Plot was then used as a standard curve for extrapolating the number of molecules of EVX1AS in the cells. Gecko tissue dissection Three pre-hatching geckos at embryonic days E50 to E55 were anesthetized by hypothermia. After the opening of the egg, the embryos were sacrificed by cervical dislocation and a collection of tissues were obtained and flash frozen in liquid nitrogen: brain, heart, eyes, gut, lungs, liver, tail and the remaining carcase (containing mainly bones, muscle and skin). RNA tethering CRISPR-display was performed as previously described30. Briefly, human and gecko lncRNAs fused to a U1 3’box at their 3’end and to a scaffold at their 5’ end were ordered as gBlocks (IDT). Subsequently, sgRNAs targeting the human or chicken EVX1 TSS (available under request) were introduced by PCR and the whole construct was cloned into a pCMV plasmid. Then, the fusion lncRNA-sgRNA constructs were co-transfected with the catalytically inactive dCas9 into SHSY5Y cells or chicken embryos (see following section). All generated plasmids are available upon request Cells The neuroepithelial SHSY5Y cell line (CRL-2266) was purchased from ATCC (Manassas, VA, US). Cells were cultured in 50% EMEM and 50% F12 medium supplemented with 10 % FBS (Millipore, Burlington, MA, USA #S0115), 100 units/ml penicillin and 100 μg/ml streptomycin (Lonza, #17-602E). For overexpression experiments 250 ng of each plasmid were used. 150000 cells/well were seeded and transfected using X-TremeGENE HP DNA transfection reagent (Sigma-Aldrich, #6366546001), cells were harvested after 48 h. 17 Animals All animal experiments were approved by a local ethical review committee and conducted in accordance with personal and project licenses in compliance with the current normative standards of the European Union (Directive 2010/63/EU) and the Spanish Government (Royal Decrees 1201/2005 and 53/2013, Law 32/107). Fertilized hen eggs (Gallus gallus), obtained from Granja Santa Isabel (Córdoba, Spain), were incubated at 37.5 °C in a humidified atmosphere until the required stages (Bellairs and Osmond, 2014). The day when eggs were incubated was considered embryonic day E0. Gecko eggs were harvested from a local breeder colony of Madagascar ground geckos (Paroedura pictus) at Achucarro (based on the colony at the Department of Ecology of Charles University, Czech Republic). Adult geckos were maintained on a 12/12-h light/dark and temperature cycle (8 a.m. lights on; 28 ° C diurnal temperature, 23 ° C nocturnal temperature) and provided with ad libitum access to food and water. Eggs were incubated at 28 ° C in a lowhumidified atmosphere until the required stages 12. The day when eggs were found in the terrarium was considered E0. In ovo electroporation Electroporation of chick embryos was performed as previously described52. Briefly, eggs were incubated in a vertical position at 38°C. Plasmids were injected with a volume of less than 1 μl into the fourth ventricle of E3 chick embryos using a fine pulled glass needle. Four electric pulses (14-17 V, 15 ms pulses with a 950 ms interval: BTX electroporator ECM) were then applied to the brain between insulated silver 40mm x 0.8mm wire electrodes with flattened pole (Intracel). Drops of Ringer’s solution supplemented with antibiotics (penicillin/streptamycin: Sigma) were added to the egg. Embryos were incubated until E6, when tissue was harvested for further research. RNA extraction and RT-QPCR RNA from all samples was extracted using Direct-zol RNA miniprep kit (Zymo research, Irvine, USA, #R2053) with DNAse treatment. For the extraction of RNA from gecko and chicken tissues samples were homogenized with a pellet pestle prior to extraction. 18 500-1000 ng of RNA were used for the retrotranscription reaction using iScript cDNA Synthesis Kit (BioRad, CA, USA, #1708890). Expression values were determined by qPCR using Sybr Green (iTaq SYBR Green Supermix, Bio-Rad, #1725124) and specific primers. RPLP0 gene was used as endogenous control in human and in gecko and Rplp7 in chicken samples. Reactions were run in a BioRad CFX384 and melting curves were analyzed to ensure the amplification of a single product. All qPCR measurements were performed in duplicate and expression levels were analyzed using the 2–∆∆Ct method. To reduce the variability across experiments, we normalized the relative expression as follows: all values from the same experiment were normalized to the highest value, hence we obtained values ranging from 0 to 1. In the case of in vivo experiments, relative expression values were normalized using z-score to allow the combination of values from different qPCRs. All primer sequences are listed in Supplementary Table 1. RNA immunoprecipitation (RIP) For RIP experiments, SHSY5Y cells were lysed in RIP buffer (150 mM KCl, 25 mM Tris, 0.5 mM DTT, 0.5 % NP-40, PI), kept on ice for 15 minutes and homogenized using a syringe. IVT lncRNAs were incubated with RNA secondary structure buffer and added to the lysates. The mixes were pre-cleared with proteinG dynabeads (ThermoFisher, Waltman, MA, USA) for 1 h in a wheel shaker at 4ºC. Pre-cleared lysates were incubated with 1ug of MED1 antibody (Santa Cruz Biotechnologies, Dallas, TX, USA) for 1h at room temperature. After incubation dynabeads were added and further incubated for 30 min. The immunoprecipitation was washed three times with RIP buffer, three times with low salt buffer (50 mM NaCl, 10 mM Tris-HCl, 0.1 % NP-40) and three times with high salt buffer (500 mM NaCl, 10 mM Tris-HCl, 0.1 % NP-40). After the washes, 70 % of beads were resuspended in RNA extraction buffer and 30 % was used for WB. Western Blot Laemmli buffer (62 mM Tris-HCl, 100 mM dithiothreitol (DTT), 10 % glycerol, 2 % SDS, 0.2 mg/ml bromophenol blue, 5 % 2-mercaptoethanol) was added to the protein samples and were denatured by heat. Proteins were migrated on 8% SDS-PAGE gels. Following electrophoresis, proteins were transferred onto nitrocellulose membranes using a Transblot-Turbo Transfer System (Biorad) and blocked in 5 % non-fatty milk diluted in TBST (20 mM Tris, 150 mM NaCl and 0.1 % Tween 20) at room temperature for 1 h. The membranes were incubated overnight at 4ºC with primary antibodies diluted 1:500 in TBST. Immunoreactive bands were revealed using 19 the Clarity Max ECL Substrate (BioRad, #1705062) after incubation with a horseradish peroxidase-conjugated anti-mouse (1:10000 dilution in 2.5 % non-fatty milk) secondary antibody for 1 h at room temperature. The immunoreactive bands were detected using a Bio-Rad Molecular Imager ChemiDoc XRS and quantified using the ImageJ software (BioRad). The following antibodies were used for Western Blotting: HSP90 (Cell Signaling; #4874) and MED1 (sc-74475). Immunohistochemistry Embryonic chick brains were fixed by immersion in PFA (4% paraformaldehyde, PFA, diluted in phosphate buffered saline 0.1M– PBS, pH 7.3). Brains were transferred to PBS 6h after fixation. Brains were sectioned in the coronal plane at 50-70 μm thickness in a vibrating microtome (Leica VT1000S). Single and double immunohistochemical reactions were performed as described previously53 using the following primary antibodies: rabbit antibody to histone H3 - phospho S10 (Abcam ab47297; 1:1000), mouse antibody to LHX5 (DSHB ;1:30) and mouse antibody to Meis2 (DSHB ; 1:30). For secondary antibodies (all 1:1000), we used Alexa 647 goat antibody to rabbit IgG (Molecular Probes, A32733) and Alexa 488 goat antibody to mouse IgG (Molecular Probes, A11001). Sections were counterstained with DAPI. Statistics and Reproducibility The data are represented as the mean ± standard error of the mean of at least three biological replicates. Mean comparisons were performed by Student’s t-test. The statistically significance level was set at p< 0.1. Code availability Although many of the analysis were carried out in web-page apps, we have coded a pipeline to unify all our bioinformatic analysis (Olazagoitia_pipeline.ipynb), to generate the boxplots for SEEKR results (seekr_output_explore.Rmd) and to carry out GENESPACE pipeline (Genespace_pipeline.Rmd). If you would like further explanations or find trouble re-running them, contact us for support. 20 Data availability The files used for the bioinformatic analyses can be found in our GitHub (https://github.com/rodrisenovilla/Olazagoitia-Garmendia). The raw data from the expression studies is available as Supplementary Data. References: 1. Rinn, J. L. & Chang, H. Y. Genome regulation by long noncoding RNAs. Annu Rev Biochem 81, (2012). 2. Fico, A., Fiorenzano, A., Pascale, E., Patriarca, E. J. & Minchiotti, G. Long non-coding RNA in stem cell pluripotency and lineage commitment: functions and evolutionary conservation. Cellular and Molecular Life Sciences vol. 76 Preprint at https://doi.org/10.1007/s00018-018-3000-z (2019). 3. Johnsson, P., Lipovich, L., Grandér, D. & Morris, K. V. Evolutionary conservation of long noncoding RNAs; sequence, structure, function. Biochim Biophys Acta 1840, 1063 (2014). 4. Diederichs, S. The four dimensions of noncoding RNA conservation. Trends in Genetics vol. 30 Preprint at https://doi.org/10.1016/j.tig.2014.01.004 (2014). 5. Szcześniak, M. W., Kubiak, M. R., Wanowska, E. & Makalowska, I. Comparative genomics in the search for conserved long noncoding RNAs. Essays Biochem 65, 741 (2021). 6. Kirk, J. M. et al. Functional classification of long non-coding RNAs by k-mer content. Nat Genet 50, (2018). 7. Sprague, D. et al. Nonlinear sequence similarity between the Xist and Rsx long noncoding RNAs suggests shared functions of tandem repeat domains. RNA 25, (2019). 8. Smith, M. A., Gesell, T., Stadler, P. F. & Mattick, J. S. Widespread purifying selection on RNA structure in mammals. Nucleic Acids Res 41, (2013). 9. Ponti, R. D., Armaos, A., Marti, S. & Tartaglia, G. G. A method for RNA structure prediction shows evidence for structure in lncRNAs. Front Mol Biosci 5, (2018). 10. Statello, L., Guo, C. J., Chen, L. L. & Huarte, M. Gene regulation by long non-coding RNAs and its biological functions. Nature Reviews Molecular Cell Biology vol. 22 Preprint at https://doi.org/10.1038/s41580-020-00315-9 (2021). 11. Much, C. et al. Evolutionary divergence of Firre localization and expression. RNA (2022) doi:10.1261/rna.079070.121. 12. Noro, M., Uejima, A., Abe, G., Manabe, M. & Tamura, K. Normal developmental stages of the Madagascar ground gecko Paroedura pictus with special reference to limb morphogenesis. Developmental Dynamics 238, (2009). 13. Nomura, T., Kawaguchi, M., Ono, K. & Murakami, Y. Reptiles: A new model for brain evo-devo research. Journal of Experimental Zoology Part B: Molecular and Developmental Evolution vol. 320 Preprint at https://doi.org/10.1002/jez.b.22484 (2013). 14. Itasaki, N., Bel-Vialar, S. & Krumlauf, R. ‘Shocking’ developments in chick embryology: Electroporation and in ovo gene expression. Nature Cell Biology vol. 1 Preprint at https://doi.org/10.1038/70231 (1999). 15. Nakamura, H. & Funahashi, J. Electroporation: Past, present and future. Development Growth and Differentiation vol. 55 Preprint at https://doi.org/10.1111/dgd.12012 (2013). 21 16. García-Moreno, F. et al. Absence of Tangentially Migrating Glutamatergic Neurons in the Developing Avian Brain. Cell Rep 22, (2018). 17. Hara, Y. et al. Madagascar ground gecko genome analysis characterizes asymmetric fates of duplicated genes. BMC Biol 16, (2018). 18. Kang, Y. J. et al. CPC2: A fast and accurate coding potential calculator based on sequence intrinsic features. Nucleic Acids Res 45, (2017). 19. GC Content Calculator - Online Analysis and Plot Tool - BiologicsCorp. https://www.biologicscorp.com/tools/GCContent/#.ZFueB3ZBw2x. 20. Frazer, K. A., Pachter, L., Poliakov, A., Rubin, E. M. & Dubchak, I. VISTA: Computational tools for comparative genomics. Nucleic Acids Res 32, (2004). 21. Rice, P., Longden, L. & Bleasby, A. EMBOSS: The European Molecular Biology Open Software Suite. Trends in Genetics vol. 16 Preprint at https://doi.org/10.1016/S01689525(00)02024-2 (2000). 22. Luo, S. et al. Divergent lncRNAs regulate gene expression and lineage differentiation in pluripotent cells. Cell Stem Cell 18, (2016). 23. Quinn, J. J. et al. Rapid evolutionary turnover underlies conserved lncRNA–genome interactions. Genes Dev 30, 191 (2016). 24. Lovell, J. T. et al. GENESPACE tracks regions of interest and gene copy number variation across multiple genomes. Elife 11, (2022). 25. Emms, D. M. & Kelly, S. OrthoFinder: Phylogenetic orthology inference for comparative genomics. Genome Biol 20, 1–14 (2019). 26. Wang, Y. et al. MCScanX: a toolkit for detection and evolutionary analysis of gene synteny and collinearity. Nucleic Acids Res 40, (2012). 27. Grabherr, M. G. et al. Genome-wide synteny through highly sensitive sequence alignment: Satsuma. Bioinformatics 26, 1145–1151 (2010). 28. Ross, C. J. et al. Uncovering deeply conserved motif combinations in rapidly evolving noncoding sequences. Genome Biol 22, 1–31 (2021). 29. Warde-Farley, D. et al. The GeneMANIA prediction server: Biological network integration for gene prioritization and predicting gene function. Nucleic Acids Res 38, (2010). 30. Shechner, D. M., Hacisuleyman, E., Younger, S. T. & Rinn, J. L. Multiplexable, locusspecific targeting of long RNAs with CRISPR-Display. Nat Methods 12, (2015). 31. Agarwala, S., Sanders, T. A. & Ragsdale, C. W. Sonic hedgehog control of size and shape in midbrain pattern formation. Science (1979) 291, (2001). 32. Fogel, J. L., Chiang, C., Huang, X. & Agarwala, S. Ventral specification and perturbed boundary formation in the mouse midbrain in the absence of hedgehog signaling. Developmental Dynamics 237, (2008). 33. Statello, L., Guo, C. J., Chen, L. L. & Huarte, M. Gene regulation by long non-coding RNAs and its biological functions. Nat Rev Mol Cell Biol 22, 96 (2021). 34. Constanty, F. & Shkumatava, A. lncRNAs in development and differentiation: From sequence motifs to functional characterization. Development (Cambridge) 148, (2021). 35. Mattick, J. S. et al. Long non-coding RNAs: definitions, functions, challenges and recommendations. Nature Reviews Molecular Cell Biology 2023 17, 1–17 (2023). 36. Wang, J. et al. Neutral evolution of ‘non-coding’ complementary DNAs. Nature 2004 431:7010 431, 1–2 (2004). 22 37. Zhou, Y. et al. Activation of paternally expressed genes and perinatal death caused by deletion of the Gtl2 gene. Development 137, (2010). 38. Stadtfeld, M. et al. Aberrant silencing of imprinted genes on chromosome 12qF1 in mouse induced pluripotent stem cells. Nature 465, (2010). 39. Bell, C. C. et al. The Evx1/Evx1as gene locus regulates anterior-posterior patterning during gastrulation. Sci Rep 6, (2016). 40. Hezroni, H. et al. Principles of Long Noncoding RNA Evolution Derived from Direct Comparison of Transcriptomes in 17 Species. Cell Rep 11, 1110–1122 (2015). 41. Pegueroles, C. et al. Transcriptomic analyses reveal groups of co-expressed, syntenic lncRNAs in four species of the genus Caenorhabditis. RNA Biol 16, 320 (2019). 42. Szcześniak, M. W., Kubiak, M. R., Wanowska, E. & Makalowska, I. Comparative genomics in the search for conserved long noncoding RNAs. Essays Biochem 65, 741 (2021). 43. Cabili, M. N. et al. Localization and abundance analysis of human lncRNAs at single-cell and single-molecule resolution. Genome Biol 16, (2015). 44. González-Moro, I. et al. A long non-coding RNA that harbors a SNP associated with type 2 diabetes regulates the expression of TGM2 gene in pancreatic beta cells. Front Endocrinol (Lausanne) 14, 216 (2023). 45. Yamaguchi, K. et al. Technical considerations in Hi-C scaffolding and evaluation of chromosome-scale genome assemblies. in Molecular Ecology vol. 30 (2021). 46. Altschul, S. F., Gish, W., Miller, W., Myers, E. W. & Lipman, D. J. Basic local alignment search tool. J Mol Biol 215, (1990). 47. Camacho, C. et al. BLAST+: Architecture and applications. BMC Bioinformatics 10, (2009). 48. Robinson, J. T. et al. Integrative genomics viewer. Nature Biotechnology vol. 29 Preprint at https://doi.org/10.1038/nbt.1754 (2011). 49. Cantalapiedra, C. P., Hern randez-Plaza, A., Letunic, I., Bork, P. & Huerta-Cepas, J. eggNOG-mapper v2: Functional Annotation, Orthology Assignments, and Domain Prediction at the Metagenomic Scale. Mol Biol Evol 38, (2021). 50. Bailey, T. L., Johnson, J., Grant, C. E. & Noble, W. S. The MEME Suite. Nucleic Acids Res 43, (2015). 51. Castellanos-Rubio, A. et al. A long noncoding RNA associated with susceptibility to celiac disease. Science (1979) 352, 91–95 (2016). 52. García-Moreno, F., Vasistha, N. A., Begbie, J. & Molnár, Z. CLoNe is a new method to target single progenitors and study their progeny in mouse and chick. Development (Cambridge) 141, (2014). 53. Rueda-Alanã, E. & Garciá-Moreno, F. Time in Neurogenesis: Conservation of the Developmental Formation of the Cerebellar Circuitry. Brain Behav Evol (2021) doi:10.1159/000519068. Acknowledgements This work was supported by grants PGC2018-097573-A-I00 from the Spanish Ministry and COLLAB20/02 from the UPV-EHU to ACR; PGC2018-096173-A-I00 grant from the Spanish Ministry and PIBA 2020_1_0057 grant from the Basque Government to FGM. 23 AOG was supported by a predoctoral fellowship from the Basque Government and a post-doctoral fellowship from the University of the Basque Country (UPV/EHU). RSG is supported by a predoctoral fellowship by Fundación Tatiana Pérez de Guzmán el Bueno. Author contributions Conceptualization: FGM, ACR Methodology: AOG, RSG, FGM, ACR Software: RSG Investigation: AOG, FGM, ACR Resources: FGM, ACR Data Curation: RSG Writing - Original Draft: AOG, RSG, FGM, ACR Writing - Review & Editing: AOG, RSG, FGM, ACR Supervision: FGM, ACR Project administration: ACR Funding acquisition FGM, ACR Competing Interests Authors declare no competing interests. Figure legends Figure 1: Nonlinear sequence similarity between human and gecko embryogenesis lncRNAs. (a) k=3 to k=6 k-mer content analysis of MEG3 and EVX1AS lncRNAs using SEEKR tool. Candidate lncRNAs were compared against the Madagascar gecko genome. The boxplot represents all gecko genes, and the dots highlight the position of the selected candidate genes as well as its correlation values (Pearson, SEEKR) and its percentile (%ile). Blue and red dots correspond to the MEG3 and EVX1AS genes respectively. (b) Secondary structure similarity study using CROSSALIGN tool. Optimal matching region by OBE-DTW algorithm (left) and overall structural similarity by Standard-DTW (right). Structural profiles are obtained with CROSS Global Score for the two RNAs (score >0 means a double-stranded nucleotide; <0 single-stranded). (c) Linear conservation analysis of MEG3 (top) and EVX1AS (bottom) performed using mVISTA. Human sequence is shown on the x-axis and percentage similarity to the corresponding gecko sequence on the y-axis. The graphical plot is based on sliding-window analysis of the underlying genomic alignment. A 100-bp sliding window is at 25-bp nucleotide increment is used. 24 Figure 2: Human and gecko EVX1AS are principally expressed in the brain. Heat map showing relative expression (from 0 to 1) of (a) gecko MEG3-like, EVX1AS-like and Evx1 and (b) human EVX1AS and EVX1 in a range of tissues. Human and gecko RPLP0 were used as housekeeping control and normalizations were done to the highest value in each experiment. Elements on the image were added using Biorender. (c) EVX1AS and EVX1 expression analysis in the human neuroepithelial cell line SHSY5Y. Human RPLP0 was used as housekeeping control. Data represents the media and standard error of 3 independent analyses. Figure 3: Human and gecko EVX1AS are likely to have evolved from independent ancestor sequences. (a) Genome localization of human EVX1AS and gecko EVX1AS-like genes. Top: Genomic locus for EVX1-AS human gene (green). Bottom: Locus of its functionally equivalent gecko’s gene, pp0020455 ≈ pp-EVX1AS-like (blue). Genome visualization is carried out by Integrative Genome Viewer (IGV) with Human (GRCh38/hg38, RefSeq.gtf) and P. picta (v2, BRAKER.gtf) genomes, respectively. Except for pp0020455, which retains the V1 annotation. (b) Synteny cross-correlation (Satsuma2). Human chromosome 7 sequence alignment against gecko chromosome scaffolds. Only EVX1AS-like (no match) and scaffolds with match alignments have been shown for simplicity. The colour code represents the chromosome aligning in the coloured regions (e.g. all gecko genes display dark green coloured regions due to its cross-alignment with human chromosome 7, dark green). EVX1AS related genes are indicated with a bar where they are located in the chromosome. Figure 4: Human EVX1AS and gecko EVX1AS-like share function-related characteristics. (a) Mobility assay of in vitro transcribed (IVT) human (Hs) and gecko (Pp) EVX1AS-like. A non related lncRNA was used as a control (Ctrl). (b) Cellular localization of human EVX1AS and gecko EVX1AS-like in human SHSY5Y cells (n=3) and gecko brains (n=2). PO and EVX1 were used as cytoplasmic controls. (c) Analysis of EVX1AS and EVX1AS-like copy numbers in human and gecko brains, respectively. Data represents the mean and standard error of the experiments. (d) Analysis of endogenous and in vitro transcribed EVX1AS variants levels bound to MED1 protein after immunoprecipitation. An in vitro transcribed unrelated lncRNA was used as a control. (n=3) A representative western blot of the immunoprecipitation is shown. (+p<0.1, *p<0.05, ****p<0.0001 according to one tailed Student’s t-test). Elements on the image were added using Biorender. Figure 5: EVX1AS variants are functionally convergent among species. (a) Relative RNA expression of human EVX1AS lncRNA and EVX1 upon human EVX1AS overexpression in SHSY5Y cells. (b) Relative RNA expression of gecko EVX1AS-like lncRNA and human EVX1 upon pp-EVX1AS-like overexpression in SHSY5Y cells. Human RPLP0 was used as housekeeping control and normalizations were done to the highest value in each experiment. (c) Immunohistochemistry of the chicken embryonic midbrain on coronal 25 section showing the successful transfection of the lncRNA plasmid. Transfected cells are shown in red due to their ectopic expression of red fluorescent protein RFP. LHX5 and MEIS2 patterns of expression (in green) indicate the mesencephalic region that was actually transfected. DAPI counterstain in blue. (d) Relative RNA expression of human EVX1AS lncRNA and chicken Evx1 upon human EVX1AS overexpression in chicken mesencephalon. Scale bars represent 1mm, 500 micras and 100 micras from top to bottom. (e) Relative RNA expression of gecko EVX1AS-like lncRNA and chicken Evx1 upon ppEVX1AS-like overexpression in chicken mesencephalon. Chicken Rplp7 was used as housekeeping control. Data represents the mean and standard error of 3-6 independent experiments. (+p<0.1, *p<0.05, **p<0.01, ***p<0.001, ****p<0.0001 according to one tailed Student’s t-test). HS1 and HS2: sgRNAs targeting human EVX1 TSS1 and TSS2 respectively. GG1 and GG2: sgRNAs targeting chicken Evx1 TSS1 and TSS2 respectively.