BACKGROUND
long non-coding rnas  are non-coding transcripts longer than 200 nt, which are often multiexonic and polyadenylated  <cit> .
compared to protein coding genes, lncrnas are transcribed at lower levels and tend to do so in a tissue-specific manner, which hampers their study and identification  <cit> .
so far, every search for lncrnas in a metazoan genome has resulted in hundreds to thousands of predicted lncrnas, with little overlap between studies.
to date, most predicted lncrnas remain without a known function.
nevertheless, there is a relatively small but steadily growing set of functionally characterized transcripts.
lncrnadb v <dig>  <cit> , a reference database for functionally validated lncrnas, lists  <dig> experimentally characterized human lncrnas, and for some of them, the function and molecular mechanism are well characterized.
for instance, xist is involved in x chromosome inactivation for dosage compensation  <cit> , hotair interacts with the chromatin remodeling complex mediating epigenetic modifications of dna  <cit> , h <dig> acts as a trans-regulator of imprinted genes  <cit> , and malat <dig> regulates alternative splicing and has been implicated in cancer  <cit> .
other lncrnas are only indirectly and loosely associated with a possible biological function.
for instance, a recent study listed lncrnas differentially expressed in normal and tumor samples but, for most of them, a direct implication in a biological process remains unclear  <cit> .
the lack of a clear function for most lncrnas, as well as their low levels of expression and sequence conservation, has led some authors to suggest that most lncrnas may actually represent transcriptional “noise,” i.e., the result of non-specific transcription  <cit> .
validating this interpretation requires the assessment of selective constraints acting on human lncrnas with a validated function.
however, most previous studies have considered lncrnas as a whole.
generally, these studies have found that, at the sequence level, lncrnas are overall much less conserved than protein coding genes in all studied organisms  <cit> .
hallmarks of selection have been found in some organisms when comparing patterns of sequence variation in introns and exons of lncrnas.
for instance, a recent study detected selective pressures acting on lncrnas of drosophila melanogaster using both polymorphism and inter-specific conservation data  <cit> .
for humans, by contrast, differences were weak or not significant   <cit> .
the authors suggested that due to the small human effective population size, selection is not strong enough to efficiently purge mutations on lncrnas.
despite this, other studies have found that exons are more conserved than introns in human lncrnas  <cit> .
finally, some studies have noted that the lack of conservation is not constant across the entire sequence and that some lncrnas contain highly conserved regions present across distant species .
a recent study showed that >85 % of lncrnas had conserved splice sites that can be dated back to the divergence of placental mammals, despite a fast turnover of exons and introns  <cit> .
it has been argued that these and other highly conserved elements may be related with the function of lncrnas.
alternatively, however, these elements may play a role at the dna level.
secondary structure may be key for the function of lncrnas, as supported by several independent analyses of some of the functionally characterized lncrnas.
for instance, in malat <dig>  a highly conserved uracil-rich region contributes to rna stability through the formation of a triple helix  <cit> .
it has also been shown that the tumor suppressor function of the lncrna meg <dig> can be attributed to two secondary fold motifs  <cit> .
some studies have found that specific folds in some lncrnas, such as sra  <cit>  and hotair  <cit> , are conserved in distant species as a result of compensatory mutations.
at the large scale, a genome-wide study based on  <dig> mammals detected that roughly 14 % of the homo sapiens genome can fold into structures that are evolutionarily conserved and that most of them  fall in regions of low sequence conservation  <cit> .
in addition, lncrnas have been found to be stable as measured by their half-life, suggesting widespread functionality  <cit> .
finally, it has been observed that lncrnas have a higher degree of secondary folding than predicted by chance, despite the fact that, surprisingly, lncrnas seem to be less structured than mrnas  <cit> .
taken together, there is accumulating evidence that structure may play an important role in lncrna functionality.
however, it remains to be established on a genome-wide scale whether the patterns of secondary structure can effectively constrain sequence evolution in lncrnas, particularly in species, such as human, with a low effective population size.
in conclusion, we still have a very poor understanding of how selective pressures may act on lncrnas at the sequence and structural levels.
several key questions remain open that are central to the understanding of the evolution and function of lncrnas.
for instance, what are the signatures of selection in those lncrnas which are known to have a function?
what role does lncrna secondary structure play in shaping sequence variation?
and, finally, what fraction of annotated human lncrnas is functional?
to address these questions and gain further insights into what evolutionary pressures may be acting on lncrnas, it is essential to combine evolutionary analyses at different levels.
firstly, inter- and intra-species level comparisons provide different degrees of resolution and are differentially affected by typical confounding factors such as the difficulties in aligning non-coding sequences.
secondly, given the lower sequence complexity of rnas and their ability to maintain conserved structures despite high sequence variation, we consider it important to account for possible constraints at the structural level.
finally, given that a set of truly functional human lncrnas exists, this can be exploited as a golden reference for establishing relationships between evolutionary constraints and functionality, thereby avoiding misleading comparisons with protein coding genes, whose functionality is achieved by decoding their sequence into proteins.
in this study, we focused on human intergenic lncrnas to ensure that the observed sequence constraints were not influenced by overlapping protein coding genes.
the studied lncrnas were derived from gencode  <dig>  <cit>  and were filtered with stringent criteria.
we also used a control data set of truly functional and intergenic lncrnas, consisting of 39 h. sapiens lncrnas with an experimentally characterized biological function  <cit> .
we analyzed patterns of sequence divergence, patterns of sequence polymorphism in different populations, and structural properties of these lncrnas.
in line with several previous studies, overall sequence conservation and single nucleotide polymorphism  density did not provide evidence of selection when comparing introns and exons.
finer and unprecedented analyses, however, revealed compelling evidence for purifying selection acting on functional lncrnas in all human populations studied.
firstly, conserved elements were enriched in exons as compared to introns.
secondly, using population genetics parameters, we found that exons have an excess of low frequency polymorphisms as compared to introns.
finally, we found that snps are depleted in structured regions with low accessibility.
this finding provides the first direct evidence of the impact of secondary structure in lncrnas sequence variation.
importantly, these findings were also apparent for the bulk of predicted lncrnas that remain uncharacterized, suggesting that the fraction of functional lncrnas under selective constraint in this set is not negligible.
RESULTS
exons in lncrnas are enriched in conserved elements but do not show overall higher conservation than introns
to provide a common background with previous studies using different sets of human lncrnas, we first analyzed phastcons scores in exonic and intronic regions of lncrnas and flanking protein coding genes, as well as in flanking intergenic regions.
since most human lncrnas seem to be primate-specific  <cit> , we based our analysis on scores computed using an evolutionary model specific for primates .
strikingly, in the set of predicted lncrnas  we observed that exons are significantly less conserved than introns and have similar levels of conservation as intergenic regions .
thus, compared to a previous study using a 46-vertebrate model  <cit> , we detected even fewer constraints, which may be due to the relatively poor quality of some primate genomes.
this reinforces the idea that predicted human lncrnas are in general very poorly conserved through evolution.
however, this result may be due to the presence of noisily transcribed, non-functional transcripts in the broad set, and we expect larger constraints in functionally characterized lncrnas.
indeed, a recent study using mouse  found that functional lncrnas have levels of sequence constraint similar to those observed in protein coding genes  <cit> .
however, according to the authors, some lncrnas of their functional set overlapped with protein coding genes or were classified as “protein coding” in a previous study  <cit> , which may have resulted in an overestimation of their conservation.
here we assessed conservation for the  <dig> human lncrnas with an experimentally determined function , which has been strictly filtered for any potential overlap with protein coding genes.
we found that the functional and the broad sets show different distributions of phastcons score ratios in exons and introns .
in contrast to the broad set, for functional lncrnas we observed the expected pattern that exons are more conserved than introns, although these differences are not significant.
since divergence estimates may be influenced by the presence of repeated elements, we calculated their abundance using the repeatmasker software  <cit> .
the percentage of sequences having repeats is quite similar when comparing the functional and broad sets, being slightly higher for the functional  than for the broad set .
however, for those sequences having repeats, the percentage of sequences covered by interspersed repeats is higher for the broad  than for the functional  set.
to evaluate whether these repeats are affecting our estimates, we also plotted phastcons scores for the best match  subset of sequences having the same amount of mapped repeats .
in this later subset, differences between exons and introns were also significant, confirming previous results obtained using the entire broad set .
thus, differences between the functional and the broad sets do not arise from different levels of repeated elements.
overall our results show that, contrary to what may be expected, conservation in lncrnas proven to be functional is also very weak.
this result implies that lack of inter-species conservation, as measured with this standard approach, cannot be taken as evidence of lack of functionality.
as mentioned above, it has been suggested that short and highly conserved sequence elements may be involved in the function of lncrnas, but it is as yet unclear whether these elements may play a role at the dna level  <cit> .
other authors have proposed that conservation in lncrnas is limited to splice-related motifs and that conservation in exon cores should be rare  <cit> .
these models are compatible with observations of overall low sequence conservation.
indeed, if functionality of lncrnas is conferred by short elements separated by largely unconstrained sequences, one could expect overall low conservation scores.
in addition, if the observed conserved elements are indeed involved in lncrna function, and not acting solely at the dna level, one would expect them to specifically associate with exonic regions, thereby forming part of the mature lncrna transcript.
we compared the abundance of conserved elements, which are discrete regions having high conservation scores as predicted by phastcons, in both functional and broad human lncrnas and using a multiple genome alignment of  <dig> vertebrates  <cit> .
we observed that the percentage of lncrnas covered by conserved elements is significantly higher in exons than in introns in both functional and broad data sets .
these results support the idea that selective constraints may be limited to the maintenance of a few clusters of positions, which may be involved in lncrna function by participating in structure or binding motifs present in the mature transcript.fig.
<dig> boxplots showing the percentage of exonic and intronic sequences covered by conserved elements in the functional and broad human data sets.
horizontal lines inside boxes represent the median, boxes show the interquartile range , vertical lines correspond to the highest and lowest value within  <dig> *iqr, and dots represent outliers
human lncrna exons show signatures of selection at the population level
considering the low conservation of lncrnas across species, it has been suggested that these molecules may have a high turnover and a short evolutionary lifespan  <cit> .
if that is the case, selective constraints in functional lncrnas may be stressed at the species or population level.
we first focused on differences in snp densities in exonic and intronic regions, which have been assessed before in the human african  population without finding significant differences  <cit> .
we computed the snp density in exons and introns in this and four other major human populations , european , east asian , and south asian ), which are roughly fourfold smaller than the afr population in terms of effective size  <cit> , and focused on differences between populations and between the broad and functional sets.
the observed snp density is fairly variable between populations, with the afr and sas populations having the highest and the lowest snp density, respectively , which is consistent with previous studies showing the highest genetic diversity in african populations  <cit> .
lncrnas and intergenic regions have higher snp densities, as compared to protein coding genes, and differences between them are generally not significant .
the distributions of snp densities in the functional and broad sets are not significantly different .
in the two sets, we observed that exons tend to accumulate fewer snps than introns, but differences were only significant for some populations in the broad set .
thus, our results are generally in line with those of a previous study restricted to the afr population  <cit> .
however, our results reveal that lncrnas with a known function display similarly low differences in snp densities between exons and introns; therefore, this feature cannot be used as evidence for a lack of functionality.
to gain a deeper insight into the selective pressures acting on human lncrnas, we performed a more thorough analysis by estimating several population genetics parameters, including nucleotide diversity , derived allele frequency , and tajima's d. nucleotide diversity  is defined as the average number of pairwise nucleotide differences per site  <cit> .
figure 2a shows the nucleotide diversity of the two sets of human lncrnas, as well as that of surrounding protein coding genes and intergenic regions.
we made three major observations.
first, nucleotide diversity levels are different between the four categories: intergenic regions and protein coding exons show the highest and lowest levels of genetic diversity, respectively, and the broad set of lncrnas has higher values than the functional set.
second, levels of nucleotide diversity vary among populations, and they can be ordered from highest to lowest levels , and the order is the same in the four categories studied.
of note, the lowest levels of snp density in the sas population are not related with the lowest π levels, since sas has higher π levels than eur and eas populations.
third, we observed, for the first time in human populations, that nucleotide diversity is significantly smaller in exons than in introns in both functional and broad lncrna sets.
we also evaluated whether the differential levels of repeats in the functional and broad sets are biasing our results, computing π for a subset of broad lncrnas having the same amount of mapped repeats .
the levels of π are similar to those for the whole set and are significantly lower in exons compared to introns, indicating that the differential composition of repeats in the sets is not biasing our results .
overall, in human lncrnas, snp density and nucleotide diversity seem to be subjected to different degrees of constraint, and only nucleotide diversity has robust detectable differences between exonic and intronic sequences.fig.
<dig>
a mean nucleotide variability  for exonic and intronic regions of the functional and broad lncrna sets in human and nearby protein coding genes and intergenic regions.
error bars represent the standard error of the mean.
b tajima's d values for exonic and intronic regions for the functional and broad lncrna in human and nearby protein coding genes and intergenic regions.
tajima's d values were computed for observed values  and coalescence simulations .
horizontal lines inside boxes represent the medians, boxes show the interquartile range , vertical lines correspond to the highest and lowest values within  <dig> *iqr, and dots represent outliers
to further evaluate whether the observed genetic diversity patterns deviate from neutrality expectations, we performed tajima's d tests  <cit> .
tajima's d is calculated as the difference between two measures of genetic diversity: the mean number of pairwise differences and the number of segregating sites, each scaled so that they are expected to be equal in a neutrally evolving population of constant size.
tajima's d was calculated for each data set  and for coalescent simulations that were computed using the observed population mutation rate value  for each region and a basic model  with the ms program  <cit> .
tajima's d values were negative in the four data sets  and in all five populations studied .
tajima's d values in the broad_bm subset were similar to those for the whole broad set, indicating that the differential composition of repeats in the sets is not biasing our estimates .
the observed tajima's d values are different from those obtained in the coalescence simulations, supporting the hypothesis that observed values deviate from neutral expectations due to an excess of polymorphism at low frequency.
the bias towards low frequency variants in lncrnas was confirmed in both exonic and intronic regions when evaluating the daf .
deviations from neutral expectations may be interpreted as the consequence of a recent population bottleneck and/or purifying selection.
human populations are known to have undergone a recent expansion  <cit> , which may contribute to the negative tajima's d values detected in all regions studied, including intronic and surrounding intergenic regions.
however, we also detected that π is not uniformly distributed in exonic and intronic regions and also not between lncrnas, protein coding genes, and intergenic regions.
thus, selective constraints contribute to the observed deviations from neutral expectations.
taken together, our results suggest that purifying selection may be acting on human lncrnas to prevent the accumulation of deleterious mutations, in both the functional and broad sets.
secondary structure constrains sequence variation in lncrnas
it has been proposed that some lncrnas may function through the adoption of specific secondary structure folds  <cit> .
in a previous study, the presence of a high number of correlated positions on multiple alignments was interpreted as evidence of evolutionary conservation of rna secondary structures  <cit> .
we evaluated the secondary structure of human lncrnas, rrna, mrna, and intergenic regions using accessibility scores calculated with two independent methods, which indicate the probability that each site belongs to an unpaired region according to an ensemble of computationally predicted secondary structures .
rrnas should be considered as a positive control, since their functionality is known to depend on their secondary structure.
by contrast, intergenic regions should be considered as a negative control, since their function  is not expected to be driven by their rna secondary structure.
although the function of mrnas depends primarily on the encoded protein, protein coding transcript sequences have been shown to be constrained at the structural level  <cit> .
regardless of the method used to calculate accessibilities, all data sets had similar distributions of residue accessibility, in which non-accessible residues likely to be paired or close to paired residues constitute the largest fraction .
firstly, we evaluated whether conserved positions  and non-conserved positions have different accessibilities.
the distributions of accessibilities in conserved and non-conserved positions are significantly different in the functional set  but not in the broad set.
however, when computing the median accessibilities for conserved and non-conserved positions for each lncrna, differences remain significant only for the sfold method .
these results suggest that conserved elements may be enriched in secondary structure folds, which in turn may be related to their function.
secondly, to evaluate whether the secondary structure influences the location of snps, we calculated the prevalence of polymorphic sites at positions with different accessibilities.
we observed that positions of low accessibility showed lower probabilities of having snps .
importantly, in the rrna, functional, broad, and mrna data sets, the differences between the distributions of positions with snp or without them were significant and larger in the range of positions with very low accessibilities  than in the rest of the accessibility ranges, independent of the method used to calculate accessibilities .
these low accessibility positions are likely to be paired or close to paired residues and constitute the largest fraction .
note that accessibilities independently computed using the two different softwares behave in the same way for all sets, the only exception being the intergenic regions.
according to the rnafold program intergenic regions do not show a tendency to prevent the accumulation of snps in low accessibilities, while according to the sfold program the behavior of the intergenic regions is similar to that of the broad and mrna regions.
these results suggest that the secondary structures predicted in the intergenic regions should be considered with caution.
importantly, both programs show that the differences between this particular range of accessibilities and others are particularly stressed in both the rrna and the functional sets.
this indicates that, overall, snps are prevented from accumulating in positions of low accessibility, that is, positions in paired regions that participate in the formation of secondary structure folds, and therefore may be key in achieving their function.fig.
<dig> density plots showing the accessibility distribution for positions containing or not containing an snp in the five major populations: african , ad mixed american , european , east asian , and south asian .
accessibility was computed using the sfold  and the rnafold  programs
fig.
<dig> difference between accessibility distributions of positions with or without snp within a given range of accessibilities .
probabilities within ranges were calculated using the integrate.xy function on a density distribution .
vertical lines represent the confidence intervals estimated using a bootstrapping after  <dig> replicates.
accessibility was computed using the sfold  and the rnafold  programs using snps from the african  population
some of the lncrnas may be partially annotated, and this may affect the predictions of the secondary structure.
thus, we selected a subset of putative full-length transcripts by keeping those that had the same length in gencode  <dig> and  <dig>  which is the latest release.
the subsets resulted in  <dig> out of  <dig> for the functional lncrna set and  <dig> out of  <dig> for the broad lncrna set.
in both cases we detected the same trend as obtained when using the whole data set, with snps prevented from accumulating in regions with low accessibility .
thus, the presence of partially annotated genes does not seem to affect our estimates of accessibility.
to evaluate whether our results are biased due to the nucleotide composition of the sequence context, we compared gc content  with the mean number of snps and the accessibility scores .
the three parameters  were calculated for non-overlapping windows of five nucleotides.
as expected, we observed a negative correlation between % gc and accessibility, confirming previous results  <cit> .
importantly, the mean number of snps remains similar for different values of % gc, indicating that the observed depletion of snps in low accessibility sites does not depend on gc content.
previous studies showed that purifying selection is maintaining a splice-related motif, i.e., an exonic splicing enhancer , near exon boundaries to ensure an efficient splicing of multiexonic lncrna  <cit> .
schüler et al.
<cit>  concluded that purifying selection acts to maintain ese motifs but not necessarily rna folding, since they failed to find a correlation between evolutionary rate and secondary structure stability.
in our study we detected that snp density is lower in ese motifs than in non-ese regions, and differences were significant for the broad set in the five populations studied , providing additional support to the idea that constraints are larger in ese than in non-ese regions.
we wanted to test whether the observed relationship between accessibility and snp density is due to the presence of ese motifs, which may point to splicing as the main factor driving the observed relationships between conservation and structure.
to this end we classified the positions of lncrnas according to the presence or not of ese motifs, and we compared the accessibility distributions for positions not having and having snps .
overall the behavior of the sites with or without annotated eses is similar for both the sfold and rnafold programs, although in the broad set differences are higher for the ese positions in all populations studied.
thus, the reduction of snps in positions of low accessibility cannot be solely explained by the presence of ese elements.
altogether, our results suggest that secondary structure constrains ancient and recent sequence variation in lncrnas, and that this is largely independent of the presence of known motifs involved in splicing.
finally, an alternative way to measure whether snps that impair folding are purged by natural selection is to estimate the impact of the variation on the energetic stability of the fold.
we did so by comparing the minimal gibbs free energy  of the reference structure and the structure of the lncrna having a certain snp, as reported in the lncrnasnp database  <cit>  .
the density plots are significantly different in the two sets .
notably, in the functional data set, median values of the change in minimal energy are narrowly centered around zero, suggesting that snps located in functional lncrnas do not generally affect the stability of the secondary structure.
conversely, in the broad set, energy changes are shifted to positive values, suggesting that snps accumulated in these lncrnas may result in less stable secondary structures.
to our best knowledge, this is the first study that provides compelling evidence for an impact of secondary structure on lncrna sequence variation.fig.
<dig>
a diagram showing how median Δg was calculated for each lncrna, which is based on the Δg of the native structure and the structure with snps .
b density plot showing the median values of Δg for the functional  and broad  human sets
CONCLUSIONS
we have found evidence of selection acting on lncrnas at both sequence and structural levels.
when evaluating divergence data, which include ancient events, we observed that exons are observably but not significantly more conserved in exons compared to introns in the functional set.
interestingly, in both functional and broad sets, we observed a significant enrichment of conserved elements in exonic regions which may be related with lncrna functionality.
when evaluating more recent events using sequence polymorphisms, we found evidence that purifying selection prevents increases in the frequency of slightly deleterious mutations, especially in exonic regions, in both functional and broad sets.
furthermore, in lncrnas with an experimentally characterized function we found that structural regions with low accessibility are more likely to be conserved.
in addition, we observed that in lncrnas, mrnas, and rrnas, segregating sites are prevented from accumulating in low accessibility, paired regions, and snps in functional lncrnas had little impact on the stability of the secondary structure.
importantly, these results are independent of the gc content, the presence of ese motifs, and the presence of partial sequences.
taken together, these results suggest that, overall, lncrna structure introduces constraints on the evolution of its sequence.
we have observed that functional and broad human lncrnas have different evolutionary constraints, although in both sets nucleotide diversity is driven by recent purifying selection.
the functional set is generally more conserved, especially in exons, and secondary structure may be maintained through constraints on snp location.
in the broad set, selective constraints are generally weaker at both the sequence and secondary structure levels.
despite these overall differences, it is difficult to predict the functionality of an individual lncrna based on the observed sequence or structural constraints, since there is a great variation in each of these single values.
this indicates that the set of functionally characterized human lncrnas is a heterogeneous group, with respect to their evolutionary signatures.
heterogeneity in the functional set may be a consequence of the different functions in which they are involved.
note that, for most parameters studied, the functional and broad sets have overlapping distributions, suggesting that numerous lncrnas of the broad set may be functional.
in summary, our study provides new evidence that lncrnas are subjected to purifying selection in human populations, and therefore numerous predicted lncrnas are potentially functional.
in addition we found first evidence that secondary structure of lncrnas shapes recent sequence variation.
in general, conservation is low in lncrnas exons but remains detectable in short, discrete regions, which have a higher tendency to participate in structural folds.
altogether our results support a model in which the functionality of lncrnas can be maintained despite large sequence divergence, probably by maintaining the presence of short elements, likely involved in folding and other forms of functionality, which are surrounded by loosely constrained regions that may act as spacers.
future experimental analyses are needed to determine whether those short conserved regions are actually functional in the mature lncrna.

 long non-coding rnas are non-coding transcripts longer than 200 nt, which are often multiexonic and polyadenylated. compared to protein coding genes, lncrnas are transcribed at lower levels and tend to do so in a tissue-specific manner, which hampers their study and identification. so far, every search for lncrnas in a metazoan genome has resulted in hundreds to thousands of predicted lncrnas, with little overlap between studies. to date, most predicted lncrnas remain without a known function. nevertheless, there is a relatively small but steadily growing set of functionally characterized transcripts. lncrnadb v, a reference database for functionally validated lncrnas, lists experimentally characterized human lncrnas, and for some of them, the function and molecular mechanism are well characterized. for instance, xist is involved in x chromosome inactivation for dosage compensation, hotair interacts with the chromatin remodeling complex mediating epigenetic modifications of dna, h acts as a trans-regulator of imprinted genes, and malat regulates alternative splicing and has been implicated in cancer. other lncrnas are only indirectly and loosely associated with a possible biological function. for instance, a recent study listed lncrnas differentially expressed in normal and tumor samples but, for most of them, a direct implication in a biological process remains unclear. the lack of a clear function for most lncrnas, as well as their low levels of expression and sequence conservation, has led some authors to suggest that most lncrnas may actually represent transcriptional “noise,” i.e., the result of non-specific transcription. validating this interpretation requires the assessment of selective constraints acting on human lncrnas with a validated function. however, most previous studies have considered lncrnas as a whole. generally, these studies have found that, at the sequence level, lncrnas are overall much less conserved than protein coding genes in all studied organisms. hallmarks of selection have been found in some organisms when comparing patterns of sequence variation in introns and exons of lncrnas. for instance, a recent study detected selective pressures acting on lncrnas of drosophila melanogaster using both polymorphism and inter-specific conservation data. for humans, by contrast, differences were weak or not significant. the authors suggested that due to the small human effective population size, selection is not strong enough to efficiently purge mutations on lncrnas. despite this, other studies have found that exons are more conserved than introns in human lncrnas. finally, some studies have noted that the lack of conservation is not constant across the entire sequence and that some lncrnas contain highly conserved regions present across distant species. a recent study showed that >85 % of lncrnas had conserved splice sites that can be dated back to the divergence of placental mammals, despite a fast turnover of exons and introns. it has been argued that these and other highly conserved elements may be related with the function of lncrnas. alternatively, however, these elements may play a role at the dna level. secondary structure may be key for the function of lncrnas, as supported by several independent analyses of some of the functionally characterized lncrnas. for instance, in malat a highly conserved uracil-rich region contributes to rna stability through the formation of a triple helix. it has also been shown that the tumor suppressor function of the lncrna meg can be attributed to two secondary fold motifs. some studies have found that specific folds in some lncrnas, such as sra and hotair, are conserved in distant species as a result of compensatory mutations. at the large scale, a genome-wide study based on mammals detected that roughly 14 % of the homo sapiens genome can fold into structures that are evolutionarily conserved and that most of them fall in regions of low sequence conservation. in addition, lncrnas have been found to be stable as measured by their half-life, suggesting widespread functionality. finally, it has been observed that lncrnas have a higher degree of secondary folding than predicted by chance, despite the fact that, surprisingly, lncrnas seem to be less structured than mrnas. taken together, there is accumulating evidence that structure may play an important role in lncrna functionality. however, it remains to be established on a genome-wide scale whether the patterns of secondary structure can effectively constrain sequence evolution in lncrnas, particularly in species, such as human, with a low effective population size. in conclusion, we still have a very poor understanding of how selective pressures may act on lncrnas at the sequence and structural levels. several key questions remain open that are central to the understanding of the evolution and function of lncrnas. for instance, what are the signatures of selection in those lncrnas which are known to have a function? what role does lncrna secondary structure play in shaping sequence variation? and, finally, what fraction of annotated human lncrnas is functional? to address these questions and gain further insights into what evolutionary pressures may be acting on lncrnas, it is essential to combine evolutionary analyses at different levels. firstly, inter- and intra-species level comparisons provide different degrees of resolution and are differentially affected by typical confounding factors such as the difficulties in aligning non-coding sequences. secondly, given the lower sequence complexity of rnas and their ability to maintain conserved structures despite high sequence variation, we consider it important to account for possible constraints at the structural level. finally, given that a set of truly functional human lncrnas exists, this can be exploited as a golden reference for establishing relationships between evolutionary constraints and functionality, thereby avoiding misleading comparisons with protein coding genes, whose functionality is achieved by decoding their sequence into proteins. in this study, we focused on human intergenic lncrnas to ensure that the observed sequence constraints were not influenced by overlapping protein coding genes. the studied lncrnas were derived from gencode and were filtered with stringent criteria. we also used a control data set of truly functional and intergenic lncrnas, consisting of 39 h. sapiens lncrnas with an experimentally characterized biological function. we analyzed patterns of sequence divergence, patterns of sequence polymorphism in different populations, and structural properties of these lncrnas. in line with several previous studies, overall sequence conservation and single nucleotide polymorphism density did not provide evidence of selection when comparing introns and exons. finer and unprecedented analyses, however, revealed compelling evidence for purifying selection acting on functional lncrnas in all human populations studied. firstly, conserved elements were enriched in exons as compared to introns. secondly, using population genetics parameters, we found that exons have an excess of low frequency polymorphisms as compared to introns. finally, we found that snps are depleted in structured regions with low accessibility. this finding provides the first direct evidence of the impact of secondary structure in lncrnas sequence variation. importantly, these findings were also apparent for the bulk of predicted lncrnas that remain uncharacterized, suggesting that the fraction of functional lncrnas under selective constraint in this set is not negligible.  exons in lncrnas are enriched in conserved elements but do not show overall higher conservation than introns to provide a common background with previous studies using different sets of human lncrnas, we first analyzed phastcons scores in exonic and intronic regions of lncrnas and flanking protein coding genes, as well as in flanking intergenic regions. since most human lncrnas seem to be primate-specific, we based our analysis on scores computed using an evolutionary model specific for primates. strikingly, in the set of predicted lncrnas we observed that exons are significantly less conserved than introns and have similar levels of conservation as intergenic regions. thus, compared to a previous study using a 46-vertebrate model, we detected even fewer constraints, which may be due to the relatively poor quality of some primate genomes. this reinforces the idea that predicted human lncrnas are in general very poorly conserved through evolution. however, this result may be due to the presence of noisily transcribed, non-functional transcripts in the broad set, and we expect larger constraints in functionally characterized lncrnas. indeed, a recent study using mouse found that functional lncrnas have levels of sequence constraint similar to those observed in protein coding genes. however, according to the authors, some lncrnas of their functional set overlapped with protein coding genes or were classified as “protein coding” in a previous study, which may have resulted in an overestimation of their conservation. here we assessed conservation for the human lncrnas with an experimentally determined function, which has been strictly filtered for any potential overlap with protein coding genes. we found that the functional and the broad sets show different distributions of phastcons score ratios in exons and introns. in contrast to the broad set, for functional lncrnas we observed the expected pattern that exons are more conserved than introns, although these differences are not significant. since divergence estimates may be influenced by the presence of repeated elements, we calculated their abundance using the repeatmasker software. the percentage of sequences having repeats is quite similar when comparing the functional and broad sets, being slightly higher for the functional than for the broad set. however, for those sequences having repeats, the percentage of sequences covered by interspersed repeats is higher for the broad than for the functional set. to evaluate whether these repeats are affecting our estimates, we also plotted phastcons scores for the best match subset of sequences having the same amount of mapped repeats. in this later subset, differences between exons and introns were also significant, confirming previous results obtained using the entire broad set. thus, differences between the functional and the broad sets do not arise from different levels of repeated elements. overall our results show that, contrary to what may be expected, conservation in lncrnas proven to be functional is also very weak. this result implies that lack of inter-species conservation, as measured with this standard approach, cannot be taken as evidence of lack of functionality. as mentioned above, it has been suggested that short and highly conserved sequence elements may be involved in the function of lncrnas, but it is as yet unclear whether these elements may play a role at the dna level. other authors have proposed that conservation in lncrnas is limited to splice-related motifs and that conservation in exon cores should be rare. these models are compatible with observations of overall low sequence conservation. indeed, if functionality of lncrnas is conferred by short elements separated by largely unconstrained sequences, one could expect overall low conservation scores. in addition, if the observed conserved elements are indeed involved in lncrna function, and not acting solely at the dna level, one would expect them to specifically associate with exonic regions, thereby forming part of the mature lncrna transcript. we compared the abundance of conserved elements, which are discrete regions having high conservation scores as predicted by phastcons, in both functional and broad human lncrnas and using a multiple genome alignment of vertebrates. we observed that the percentage of lncrnas covered by conserved elements is significantly higher in exons than in introns in both functional and broad data sets. these results support the idea that selective constraints may be limited to the maintenance of a few clusters of positions, which may be involved in lncrna function by participating in structure or binding motifs present in the mature transcript.fig. boxplots showing the percentage of exonic and intronic sequences covered by conserved elements in the functional and broad human data sets. horizontal lines inside boxes represent the median, boxes show the interquartile range, vertical lines correspond to the highest and lowest value within *iqr, and dots represent outliers human lncrna exons show signatures of selection at the population level considering the low conservation of lncrnas across species, it has been suggested that these molecules may have a high turnover and a short evolutionary lifespan. if that is the case, selective constraints in functional lncrnas may be stressed at the species or population level. we first focused on differences in snp densities in exonic and intronic regions, which have been assessed before in the human african population without finding significant differences. we computed the snp density in exons and introns in this and four other major human populations, european, east asian, and south asian ), which are roughly fourfold smaller than the afr population in terms of effective size, and focused on differences between populations and between the broad and functional sets. the observed snp density is fairly variable between populations, with the afr and sas populations having the highest and the lowest snp density, respectively, which is consistent with previous studies showing the highest genetic diversity in african populations. lncrnas and intergenic regions have higher snp densities, as compared to protein coding genes, and differences between them are generally not significant. the distributions of snp densities in the functional and broad sets are not significantly different. in the two sets, we observed that exons tend to accumulate fewer snps than introns, but differences were only significant for some populations in the broad set. thus, our results are generally in line with those of a previous study restricted to the afr population. however, our results reveal that lncrnas with a known function display similarly low differences in snp densities between exons and introns; therefore, this feature cannot be used as evidence for a lack of functionality. to gain a deeper insight into the selective pressures acting on human lncrnas, we performed a more thorough analysis by estimating several population genetics parameters, including nucleotide diversity, derived allele frequency, and tajima's d. nucleotide diversity is defined as the average number of pairwise nucleotide differences per site. figure 2a shows the nucleotide diversity of the two sets of human lncrnas, as well as that of surrounding protein coding genes and intergenic regions. we made three major observations. first, nucleotide diversity levels are different between the four categories: intergenic regions and protein coding exons show the highest and lowest levels of genetic diversity, respectively, and the broad set of lncrnas has higher values than the functional set. second, levels of nucleotide diversity vary among populations, and they can be ordered from highest to lowest levels, and the order is the same in the four categories studied. of note, the lowest levels of snp density in the sas population are not related with the lowest π levels, since sas has higher π levels than eur and eas populations. third, we observed, for the first time in human populations, that nucleotide diversity is significantly smaller in exons than in introns in both functional and broad lncrna sets. we also evaluated whether the differential levels of repeats in the functional and broad sets are biasing our results, computing π for a subset of broad lncrnas having the same amount of mapped repeats. the levels of π are similar to those for the whole set and are significantly lower in exons compared to introns, indicating that the differential composition of repeats in the sets is not biasing our results. overall, in human lncrnas, snp density and nucleotide diversity seem to be subjected to different degrees of constraint, and only nucleotide diversity has robust detectable differences between exonic and intronic sequences.fig. a mean nucleotide variability for exonic and intronic regions of the functional and broad lncrna sets in human and nearby protein coding genes and intergenic regions. error bars represent the standard error of the mean. b tajima's d values for exonic and intronic regions for the functional and broad lncrna in human and nearby protein coding genes and intergenic regions. tajima's d values were computed for observed values and coalescence simulations. horizontal lines inside boxes represent the medians, boxes show the interquartile range, vertical lines correspond to the highest and lowest values within *iqr, and dots represent outliers to further evaluate whether the observed genetic diversity patterns deviate from neutrality expectations, we performed tajima's d tests. tajima's d is calculated as the difference between two measures of genetic diversity: the mean number of pairwise differences and the number of segregating sites, each scaled so that they are expected to be equal in a neutrally evolving population of constant size. tajima's d was calculated for each data set and for coalescent simulations that were computed using the observed population mutation rate value for each region and a basic model with the ms program. tajima's d values were negative in the four data sets and in all five populations studied. tajima's d values in the broad_bm subset were similar to those for the whole broad set, indicating that the differential composition of repeats in the sets is not biasing our estimates. the observed tajima's d values are different from those obtained in the coalescence simulations, supporting the hypothesis that observed values deviate from neutral expectations due to an excess of polymorphism at low frequency. the bias towards low frequency variants in lncrnas was confirmed in both exonic and intronic regions when evaluating the daf. deviations from neutral expectations may be interpreted as the consequence of a recent population bottleneck and/or purifying selection. human populations are known to have undergone a recent expansion, which may contribute to the negative tajima's d values detected in all regions studied, including intronic and surrounding intergenic regions. however, we also detected that π is not uniformly distributed in exonic and intronic regions and also not between lncrnas, protein coding genes, and intergenic regions. thus, selective constraints contribute to the observed deviations from neutral expectations. taken together, our results suggest that purifying selection may be acting on human lncrnas to prevent the accumulation of deleterious mutations, in both the functional and broad sets. secondary structure constrains sequence variation in lncrnas it has been proposed that some lncrnas may function through the adoption of specific secondary structure folds. in a previous study, the presence of a high number of correlated positions on multiple alignments was interpreted as evidence of evolutionary conservation of rna secondary structures. we evaluated the secondary structure of human lncrnas, rrna, mrna, and intergenic regions using accessibility scores calculated with two independent methods, which indicate the probability that each site belongs to an unpaired region according to an ensemble of computationally predicted secondary structures. rrnas should be considered as a positive control, since their functionality is known to depend on their secondary structure. by contrast, intergenic regions should be considered as a negative control, since their function is not expected to be driven by their rna secondary structure. although the function of mrnas depends primarily on the encoded protein, protein coding transcript sequences have been shown to be constrained at the structural level. regardless of the method used to calculate accessibilities, all data sets had similar distributions of residue accessibility, in which non-accessible residues likely to be paired or close to paired residues constitute the largest fraction. firstly, we evaluated whether conserved positions and non-conserved positions have different accessibilities. the distributions of accessibilities in conserved and non-conserved positions are significantly different in the functional set but not in the broad set. however, when computing the median accessibilities for conserved and non-conserved positions for each lncrna, differences remain significant only for the sfold method. these results suggest that conserved elements may be enriched in secondary structure folds, which in turn may be related to their function. secondly, to evaluate whether the secondary structure influences the location of snps, we calculated the prevalence of polymorphic sites at positions with different accessibilities. we observed that positions of low accessibility showed lower probabilities of having snps. importantly, in the rrna, functional, broad, and mrna data sets, the differences between the distributions of positions with snp or without them were significant and larger in the range of positions with very low accessibilities than in the rest of the accessibility ranges, independent of the method used to calculate accessibilities. these low accessibility positions are likely to be paired or close to paired residues and constitute the largest fraction. note that accessibilities independently computed using the two different softwares behave in the same way for all sets, the only exception being the intergenic regions. according to the rnafold program intergenic regions do not show a tendency to prevent the accumulation of snps in low accessibilities, while according to the sfold program the behavior of the intergenic regions is similar to that of the broad and mrna regions. these results suggest that the secondary structures predicted in the intergenic regions should be considered with caution. importantly, both programs show that the differences between this particular range of accessibilities and others are particularly stressed in both the rrna and the functional sets. this indicates that, overall, snps are prevented from accumulating in positions of low accessibility, that is, positions in paired regions that participate in the formation of secondary structure folds, and therefore may be key in achieving their function.fig. density plots showing the accessibility distribution for positions containing or not containing an snp in the five major populations: african, ad mixed american, european, east asian, and south asian. accessibility was computed using the sfold and the rnafold programs fig. difference between accessibility distributions of positions with or without snp within a given range of accessibilities. probabilities within ranges were calculated using the integrate.xy function on a density distribution. vertical lines represent the confidence intervals estimated using a bootstrapping after replicates. accessibility was computed using the sfold and the rnafold programs using snps from the african population some of the lncrnas may be partially annotated, and this may affect the predictions of the secondary structure. thus, we selected a subset of putative full-length transcripts by keeping those that had the same length in gencode and which is the latest release. the subsets resulted in out of for the functional lncrna set and out of for the broad lncrna set. in both cases we detected the same trend as obtained when using the whole data set, with snps prevented from accumulating in regions with low accessibility. thus, the presence of partially annotated genes does not seem to affect our estimates of accessibility. to evaluate whether our results are biased due to the nucleotide composition of the sequence context, we compared gc content with the mean number of snps and the accessibility scores. the three parameters were calculated for non-overlapping windows of five nucleotides. as expected, we observed a negative correlation between % gc and accessibility, confirming previous results. importantly, the mean number of snps remains similar for different values of % gc, indicating that the observed depletion of snps in low accessibility sites does not depend on gc content. previous studies showed that purifying selection is maintaining a splice-related motif, i.e., an exonic splicing enhancer, near exon boundaries to ensure an efficient splicing of multiexonic lncrna. schüler et al. concluded that purifying selection acts to maintain ese motifs but not necessarily rna folding, since they failed to find a correlation between evolutionary rate and secondary structure stability. in our study we detected that snp density is lower in ese motifs than in non-ese regions, and differences were significant for the broad set in the five populations studied, providing additional support to the idea that constraints are larger in ese than in non-ese regions. we wanted to test whether the observed relationship between accessibility and snp density is due to the presence of ese motifs, which may point to splicing as the main factor driving the observed relationships between conservation and structure. to this end we classified the positions of lncrnas according to the presence or not of ese motifs, and we compared the accessibility distributions for positions not having and having snps. overall the behavior of the sites with or without annotated eses is similar for both the sfold and rnafold programs, although in the broad set differences are higher for the ese positions in all populations studied. thus, the reduction of snps in positions of low accessibility cannot be solely explained by the presence of ese elements. altogether, our results suggest that secondary structure constrains ancient and recent sequence variation in lncrnas, and that this is largely independent of the presence of known motifs involved in splicing. finally, an alternative way to measure whether snps that impair folding are purged by natural selection is to estimate the impact of the variation on the energetic stability of the fold. we did so by comparing the minimal gibbs free energy of the reference structure and the structure of the lncrna having a certain snp, as reported in the lncrnasnp database. the density plots are significantly different in the two sets. notably, in the functional data set, median values of the change in minimal energy are narrowly centered around zero, suggesting that snps located in functional lncrnas do not generally affect the stability of the secondary structure. conversely, in the broad set, energy changes are shifted to positive values, suggesting that snps accumulated in these lncrnas may result in less stable secondary structures. to our best knowledge, this is the first study that provides compelling evidence for an impact of secondary structure on lncrna sequence variation.fig. a diagram showing how median Δg was calculated for each lncrna, which is based on the Δg of the native structure and the structure with snps. b density plot showing the median values of Δg for the functional and broad human sets  we have found evidence of selection acting on lncrnas at both sequence and structural levels. when evaluating divergence data, which include ancient events, we observed that exons are observably but not significantly more conserved in exons compared to introns in the functional set. interestingly, in both functional and broad sets, we observed a significant enrichment of conserved elements in exonic regions which may be related with lncrna functionality. when evaluating more recent events using sequence polymorphisms, we found evidence that purifying selection prevents increases in the frequency of slightly deleterious mutations, especially in exonic regions, in both functional and broad sets. furthermore, in lncrnas with an experimentally characterized function we found that structural regions with low accessibility are more likely to be conserved. in addition, we observed that in lncrnas, mrnas, and rrnas, segregating sites are prevented from accumulating in low accessibility, paired regions, and snps in functional lncrnas had little impact on the stability of the secondary structure. importantly, these results are independent of the gc content, the presence of ese motifs, and the presence of partial sequences. taken together, these results suggest that, overall, lncrna structure introduces constraints on the evolution of its sequence. we have observed that functional and broad human lncrnas have different evolutionary constraints, although in both sets nucleotide diversity is driven by recent purifying selection. the functional set is generally more conserved, especially in exons, and secondary structure may be maintained through constraints on snp location. in the broad set, selective constraints are generally weaker at both the sequence and secondary structure levels. despite these overall differences, it is difficult to predict the functionality of an individual lncrna based on the observed sequence or structural constraints, since there is a great variation in each of these single values. this indicates that the set of functionally characterized human lncrnas is a heterogeneous group, with respect to their evolutionary signatures. heterogeneity in the functional set may be a consequence of the different functions in which they are involved. note that, for most parameters studied, the functional and broad sets have overlapping distributions, suggesting that numerous lncrnas of the broad set may be functional. in summary, our study provides new evidence that lncrnas are subjected to purifying selection in human populations, and therefore numerous predicted lncrnas are potentially functional. in addition we found first evidence that secondary structure of lncrnas shapes recent sequence variation. in general, conservation is low in lncrnas exons but remains detectable in short, discrete regions, which have a higher tendency to participate in structural folds. altogether our results support a model in which the functionality of lncrnas can be maintained despite large sequence divergence, probably by maintaining the presence of short elements, likely involved in folding and other forms of functionality, which are surrounded by loosely constrained regions that may act as spacers. future experimental analyses are needed to determine whether those short conserved regions are actually functional in the mature lncrna. 
here we evaluate signatures of selection in human lncrnas using inter-specific data and intra-specific comparisons from five major populations, as well as by assessing relationships between sequence variation and predictions of secondary structure. we found that rna secondary structure constrains sequence variation in lncrnas, so that polymorphisms are depleted in paired regions with low accessibility and tend to be neutral with respect to structural stability. important implications of our results are that secondary structure plays a role in the functionality of lncrnas, and that the set of predicted lncrnas contains a large fraction of functional ones that may play key roles that remain to be discovered. lncrnapurifying selectiondivergencepolymorphismsecondary structurehttp://dx.doi.org/ /501100000781european research councilerc-2012-stg-310325issue-copyright-statement© the author 2016 altogether, our results show compelling evidence of recent purifying selection acting on both characterized and predicted lncrnas. this has raised the questions of what fraction of the predicted lncrnas is actually functional, and whether selection can effectively constrain lncrnas in species with small effective population sizes such as human populations.
 metazoans transcribe many long non-coding rnas that are poorly conserved and whose function remains unknown. this has raised the questions of what fraction of the predicted lncrnas is actually functional, and whether selection can effectively constrain lncrnas in species with small effective population sizes such as human populations.  here we evaluate signatures of selection in human lncrnas using inter-specific data and intra-specific comparisons from five major populations, as well as by assessing relationships between sequence variation and predictions of secondary structure. in all analyses we included a reference of functionally characterized lncrnas. altogether, our results show compelling evidence of recent purifying selection acting on both characterized and predicted lncrnas. we found that rna secondary structure constrains sequence variation in lncrnas, so that polymorphisms are depleted in paired regions with low accessibility and tend to be neutral with respect to structural stability.  important implications of our results are that secondary structure plays a role in the functionality of lncrnas, and that the set of predicted lncrnas contains a large fraction of functional ones that may play key roles that remain to be discovered. electronic supplementary material the online version of this article contains supplementary material, which is available to authorized users. keywords lncrnapurifying selectiondivergencepolymorphismsecondary structurehttp://dx.doi.org/ /501100000781european research councilerc-2012-stg-310325issue-copyright-statement© the author 2016 