Yu Hamba,Takashi Kamatani,Fuyuki Miya,Keith A. Boroevich,Tatsuhiko Tsunoda
Topologically associating domain underlies tissue specific expression of long intergenic non-coding RNAs
Computational bioinformatics,Expression study,Genomic analysis
Summary Accumulating evidence indicates that long intergenic non-coding RNAs (lincRNAs) show more tissue-specific expression patterns than protein-coding genes (PCGs). However, although lincRNAs are subject to canonical transcriptional regulation like PCGs, the molecular basis for the specificity of their expression patterns remains unclear. Here, using expression data and coordinates of topologically associating domains (TADs) in human tissues, we show that lincRNA loci are significantly enriched in the more internal region of TADs compared to PCGs and that lincRNAs within TADs have higher tissue specificity than those outside TADs. Based on these, we propose an analytical framework to interpret transcriptional status using lincRNA as an indicator. We applied it to hypertrophic cardiomyopathy data and found disease-specific transcriptional regulation: ectopic expression of keratin at the TAD level and derepression of myocyte differentiation-related genes by E2F1 with down-regulation of LINC00881. Our results provide understanding of the function and regulation of lincRNAs according to genomic structure.
Topologically associating domain underlies tissue specific expression of long intergenic non-coding RNAs

Accumulating evidence indicates that long intergenic non-coding RNAs (lincRNAs) show more tissue-specific expression patterns than protein-coding genes (PCGs). However, although lincRNAs are subject to canonical transcriptional regulation like PCGs, the molecular basis for the specificity of their expression patterns remains unclear. Here, using expression data and coordinates of topologically associating domains (TADs) in human tissues, we show that lincRNA loci are significantly enriched in the more internal region of TADs compared to PCGs and that lincRNAs within TADs have higher tissue specificity than those outside TADs. Based on these, we propose an analytical framework to interpret transcriptional status using lincRNA as an indicator. We applied it to hypertrophic cardiomyopathy data and found disease-specific transcriptional regulation: ectopic expression of keratin at the TAD level and derepression of myocyte differentiation-related genes by E2F1 with down-regulation of LINC00881. Our results provide understanding of the function and regulation of lincRNAs according to genomic structure.

A large proportion of the human genome is transcribed into RNA, yielding many non-coding RNAs (ncRNAs) in addition to mRNAs. Among ncRNAs, long non-coding RNAs (lncRNAs) are broadly defined as transcripts longer than 200 nucleotides that lack coding capacity., Accumulating evidence indicates that without being translated into proteins, lncRNAs are functional in many cellular processes, such as gene imprinting, development and immune response. Among those various mechanisms reported, many lncRNAs are associated with chromatin-modifying complexes and guide conformational changes in nuclear domains or regulation of transcriptional enhancer activity. Moreover, lncRNA dysfunction has been associated with a variety of human diseases, including cancer, cardiovascular disease, and neurodegenerative disorders. lncRNAs are classified based on their genomic origin: Antisense lncRNAs are transcribed from the opposite strand of protein-coding genes (PCGs), bidirectional eRNAs are transcribed from the enhancer region, and lincRNAs have transcription units independent of promoters and enhancers of PCGs. lincRNAs possess many characteristics of mRNAs: lincRNA transcription is regulated by key transcription factors and typical histone modifications, and lincRNA transcripts are processed by the canonical spliceosomal machinery.,,, Compared to protein-coding genes, lincRNA expression is highly specific to cell types and tissues.,, Compared to mRNA, which requires translation into protein to exert gene function, lincRNA transcripts exert their function directly on their own. This means that lincRNA transcriptional regulation serves as a functional control, emphasizing the importance of elucidating the molecular basis of specific expression. Although the mechanisms that cause tissue-specific gene expression patterns have been extensively investigated in terms of promoter sequences and epigenomic modifications, the mechanisms that can comprehensively explain the nature of lincRNAs and this difference in expression specificity between mRNAs and lincRNAs are still unclear. One of the highest order regulatory mechanisms that determine transcription patterns is genomic structure that allows for genome-wide regulation; chromatin compaction, as seen in heterochromatin regions, limits access to gene promoters. In the nucleus, chromosomes form higher-order structures that govern gene regulation. Recent advances in chromatin conformation capture technology have revealed the three-dimensional structure of the genome. Adjacent chromosomal regions were shown to fold into topological associating domains (TADs) ranging in size from a few hundred kb to 1–2 Mb. Chromatin structures that comprehensively regulate the tissue-specific expression patterns of lincRNAs must be tissue and genome-wide regulatory systems; TADs fulfill this requirement by confining enhancer-promoter interactions within each TAD, thus enabling gene expression regulation. Literature shows that lincRNAs expression are positively correlated with those of nearby genes, suggesting that regulation by enhancers insulated within TADs determine when and in which tissues lincRNAs are transcribed with neighboring genes. In this study, first, we combined TAD coordinates from 13 primary human tissues with expression data from the GTEx project to investigate into the distribution of lincRNA loci in the human genome in relation to TADs and clarify how they are associated with tissue-specific expression patterns. Next, by defining and utilizing two groups of PCGs: Co-expressed cis PCGs spatially close to the lincRNA locus and trans PCGs regulated by post-transcriptional lincRNA molecules, we proposed an analytical framework to infer biological processes relevant to traits such as tissue phenotypes and diseases. Because TAD structure is maintained throughout the cell lineage, our method should be applicable to a broad range of human gene expression data. Last, as an example for its medical applications, we analyzed hypertrophic cardiomyopathy (HCM) data, and identified lincRNA signatures specific to HCM. Overall, our findings demonstrate that interpreting the genome at the 3D level can provide insight into lincRNA functionality.

First, we investigated tissue specificity of lincRNA expression by using expression data from 54 human tissues of GTEx. t-distributed stochastic neighbor embedding (t-SNE) analysis of the expression data showed a clear clustering of samples by tissue of origin (Figure S1A), and the expression pattern of 7345 lincRNAs appeared to be at least as tissue-specific as that of 18358 PCGs (Figure S1B). To quantitatively compare tissue specificity of expression between lincRNAs and PCGs, we calculated tissue specificity scores (0 = ubiquitous; 1 = specific to a single tissue) based on the Shannon entropy specificity of each gene (Figure 1A). Overall, lincRNA has higher tissue specificity than PCG, consistent with previous findings.,, The distribution of the tissue specificity scores showed bimodal distributions with high and low specificity groups for both lincRNA and PCG. The percentage of the completely tissue-specific group of lincRNAs was larger than that of PCGs (8.81% versus 1.06%), and the tissue specificity score of lincRNAs was significantly higher than that of PCGs (p< 10−4; permutation test; Figure 1A right). This suggests that there is likely a genome-wide regulatory mechanism that results in a tissue-specific expression pattern of lincRNAs that is differentiated from the pattern of expression of PCGs. To characterize lincRNA expression in the human tissues, we conducted an agglomerative clustering analysis on the lincRNA expression profiles from the GTEx data. This resulted in the generation of three major clusters: class 1, 2, and 3 (Figure 1B). The distribution of the tissue specificity scores was lowest for class 1 and highest for class 3 (Figure 1C). Next, we examined the lincRNA expression intensities and found that class 1 was the highest and class 3 was the lowest (Figure 1D). We also calculated correlation coefficients between expression levels of lincRNA and those of neighboring genes for each cluster, and its result showed that lincRNA expression levels were positively correlated with those of neighboring genes in most cases (Figure 1E).

Using TAD coordinates detected by Hi-C on 13 primary human tissues, we tested the hypothesis that tissue-specific expression and function of lincRNAs is regulated in a manner dependent on the higher-order structure of the genome. When PCGs and lincRNAs were compared, lincRNAs were found to have a significantly higher percentage of loci within TADs across multiple tissues (Figures 2A and S2A). More interestingly, the localization of lincRNAs in TADs was more significant than those of other ncRNA subtypes such as antisense, miRNA, and rRNA. This result suggests that beyond just being a set of untranslated transcripts, ncRNAs, particularly lincRNAs, can be clearly differentiated in the context of the chromatin environment, such as TADs. To rule out the possibility that the differences of the gene locus enrichment for TADs were because of differences in gene length, we conducted a permutation test in which gene loci were randomly replaced in the genome (Figure S2B). The results showed that although the PCG loci were less enriched in TADs relative to the random distributions in almost all tissues (p< 10−4; permutation test; Figure S3A), lincRNA loci were significantly more enriched in TADs than compared to the random distributions in all tissues (p< 10−4; permutation test; Figure S3B). Incidentally, antisense ncRNAs showed a wide variation and no clear trends in their localization across tissues (Figure S3C). To examine the relationship between gene loci enrichment in TADs and tissue specificity of their expressions, we defined coordinates of TADs that could be analyzed across tissues. First, we extracted TAD regions common across all 13 tissues (Figure 2B). Although the common TAD regions had a lower overall genome coverage (total coverage: 60.1%) compared to the TAD regions in each tissue, they had a median length of 0.64M bp per region and median gene count of 10, which did not deviate significantly from the general definition for a TAD (Figures S4A and S4B). To confirm whether these regions function as effective TAD insulators, we examined the differences in correlations between gene expression levels according to the gene positions in relation to the TAD regions by using the GTEx expression data. More specifically, we calculated the moving average of the expression correlations according to the intergenic distance for three conditions: gene pairs belonging to the same TAD, those belonging to TAD and TAD boundary, and those belonging to different TADs. Gene pairs belonging to the same TAD maintained higher correlations than other pairs across long distances ∼1M bp (Figure 2C), indicating that the gene expression regulation by TADs was observed on a genome scale, supporting previous studies.,, Next, we defined the expression level of each TAD as the average of the expression levels of the included genes, and examined the correlations of the expression levels among the TADs of each chromosome using all of the GTEx expression data. They showed almost entirely positive correlations except for a small portion that had differences in TAD transcription levels among tissues (Figures S5A and S5B). These results confirmed that the common TAD regions derived from the TAD coordinates across the 13 tissues maintain the basic characteristics of TAD, which previous studies have suggested.,, Similar to the TAD regions identified from the primary tissues, the lincRNA loci were found to be strongly enriched within the common TAD regions against the random distribution by permutation, whereas PCG and antisense loci were less enriched in TAD regions (p< 10−4; permutation test; Figures 2D–2F and S3A–S3C). To confirm these results, we performed a permutation test using circular randomization that preserves the size and number of gene regions as well as their structure in genomic relationship to each other (p< 10−4; permutation test; Figure S6A). We also removed possible confounders by excluding housekeeping (HK) genes, setting a threshold for expression levels (Excluding genes with a median TPM less than 0.5 in all tissue), and restricting lincRNAs to those with experimentally proven functional roles. Because most of the HK genes are PCGs and are known to be located at the TAD boundary,, exclusion of HK genes could have affected patterns of PCG localization and expression specificity. However, the absence of HK genes did not significantly affect the trends in genomic location and specificity patterns for each gene type (p< 10−4; permutation test; Figure S6B). Filtering out low expression genes eliminated groups with particularly high tissue specificity scores (Figure S6E), which includes a subset of class 3 (Figures 1B–1E), genes with low expression levels and exclusive expression patterns. When restricted to lincRNAs with experimental functional annotations, permutation pvalue was higher than in the other conditions (p = 5.0 × 10−4; permutation test; Figure S6F) and there was no statistically significant difference in the tissue specificity scores (Figure S6G). It should be noted that the number of experimentally proven functional lincRNA genes is limited to a few hundred, and there may be an artificial bias that makes lincRNAs with large condition-specific expression variations, such as human disease, more likely to be selected for functional evaluation. In either condition, lincRNA loci were consistently localized to TADs and PCG loci were localized to TAD boundaries, with higher tissue specificity for gene groups in TADs. We also examined the tissue specificity of gene expressions inside and outside the TAD regions. Because the broad genomic environment in testis, where dramatic chromatin architectural reorganization throughout spermatogenesis and temporary disappearance of TADs are observed, does not ensure maintenance of TAD structure and regulation, we excluded the expression data of testis from the following analyses. Without testis expression data, lincRNA expression maintained higher tissue specificity than those of PCGs (Figure S7A). When lincRNAs and PCGs were divided into groups of those which were located inside and outside TADs, expressions of the groups included in TADs showed higher tissue specificities for both types of genes (Figure 2G). The relationship between the TAD localization and tissue specificity of lincRNA expressions was also examined for the three previously described classes (class 1, 2, and 3) of lincRNAs identified from their expression patterns. In class 1 and 2, which had relatively low overall tissue specificity, the expression of lincRNAs inside TADs showed significantly higher tissue specificity compared to those outside TADs (Figure S7B; class 1: p = 5.078 × 10−9, class 2: p = 3.16 × 10−3). The difference in the expression specificities between groups within and outside TADs was the largest in class 2. On the other hand, in class 3, which has the highest overall tissue specificity, there was no significant difference in expression specificities between lincRNAs inside and outside TADs. This implies that the direct contribution of localization in TADs to the tissue specificity is limited and that additional regulatory mechanisms are required to explain the increased tissue specificity. Comparison of the proportion of genes present in TADs between classes showed that class 2 had the highest ratio, class 3 had a ratio as high as class 2, and class 1 had the lowest ratio (Figure S7C; class 1: ratio = 0.594, class 2: ratio = 0.684, class 3: ratio = 0.676), suggesting that localization of loci to TADs possibly serves as a basis for downstream regulation of expression. We examined the correlation of gene expression between lincRNAs and PCGs with intergenic distance and position in the genome structure, and found that lincRNA-PCG pairs belonging to the same TAD maintained a slightly higher correlation than other positional gene pairs at distances greater than 10K bp. Of interest, contrary to the expectation that gene pairs located at the TAD boundaries would inevitably have higher correlation coefficients because of their ubiquitous expression patterns, in-TAD gene pairs located a few hundred bp closer together tended to have higher correlation coefficients than those at the TAD boundaries (Figure 2H). The distribution of the distance from the loci of lincRNAs and PCGs inside TADs to the nearest TAD boundary, respectively, showed that lincRNAs were located farther from the boundary and more internal to TADs than PCGs (p = 3.90 × 10−134; Figure 2I).

Because the expression pattern of lincRNA is more cell- and tissue-specific than that of PCG, lincRNAs could be useful and sensitive biomarkers; their expression profiles can be used for identifying expression states specific to particular biological phenomena. Therefore, we proposed the following analysis scheme for biomedical applications (Figure S8). First, lincRNAs that vary specifically for a biological phenomenon of interest are identified with a differential expression (DE) analysis. Next, the cis- and trans-target genes of lincRNAs are analyzed. Based on the results above, genes that have loci in the same TAD and show similar expression patterns to lincRNAs are identified as cis-target genes. Because lincRNAs and cis-target genes are specifically regulated by regulatory sequences insulated within TADs, gene ontology (GO) term enrichment analysis are used to identify trait-specific functions common across the cis-target genes. For trans-target genes, under consideration of their functionality, genome-wide genes that have expression patterns negatively correlated with the lincRNAs are identified. lincRNAs often regulate trans-target genes with co-factors. Therefore, if there is a common upstream regulator of the trans-target genes, it is likely to be a co-factor of the lincRNA. Further information on the subcellular localization of the lincRNAs would enable identification of the upstream factors that interact directly with lincRNAs. Because lincRNAs are polyA-modified, albeit inefficiently, they are detectable by standard RNA-Seq as well as PCGs. TADs are also maintained in most tissues and cells. Therefore, this analysis framework is applicable to a wide range of samples profiled with conventional RNA sequencing methods.

As one of biomedical applications of the above analysis scheme, we applied it with a dataset of expression profiles from HCM to identify pathologically relevant molecular profiles common to patients with heterogeneous genetic backgrounds. We identified 6,095 DEGs in human HCM by using RNA-Seq data of myocardial tissues from healthy human volunteers and myocardial tissues from HCM patients with multiple genetic backgrounds (Table S1, FDR<0.05, Figure S9A). The results of GO analysis for all DEGs showed that GO enrichment related to translational processes was detected at a top level in HCM, with decreased expression of translation-related genes consistent with a recent report, and GO terms such as mitochondria and energy production, which are important for myocardial function, also showed significant enrichment, indicating that the overall DEGs captured the characteristics of HCM (Figure S9B). Next, we identified lincRNAs with differing expression levels in HCM. We found that 132 lincRNAs were significantly up-regulated and 140 lincRNAs were significantly down-regulated in the myocardium tissues from patients compared to those from healthy volunteers (Figure 3A). We found that 113 of these lincRNAs resided within 11 TADs with altered expression levels, indicating HCM-specific cis-association by lincRNAs (Figure 3B). We performed GO enrichment analysis of these 113 lincRNAs and DEGs in the same TADs according to three patterns: DEGs in all TADs with altered expression levels, DEGs in TADs with up-regulated expression levels, and DEGs in TADs with down-regulated expression levels. Significant GO enrichment was observed only for DEGs in TADs with down-regulated expression levels, which indicated a relevance of a GO category in biological pathways related to keratinization (Figures 3C, S10A, and S10B). Among the genes used in the GO analysis, keratin associated proteins were found to be downregulated in HCM (Figure 3D). On the other hand, we could observe that some genes of the keratin family were markedly up-regulated (Figures 3E and S10C). It has been reported that the K8/K18 pair is ectopically expressed downstream of TNF-α in a desmin-deficient mouse model of heart failure and exerts a cardioprotective effect in the myocardium after tissue injury, and that K8/K18 is also ectopically expressed in human failing myocardium, where TNF-α expression is up-regulated. Because the symptoms of HCM are also caused by mechanical stress because of abnormal physiological power production in sarcomeres, it is possible that the ectopic expression of intermediate filaments acts to protect the myocardium. Ectopic expression of intermediate system filaments may have a protective effect on the myocardium. However, in HCM, we did not observe any changes in the expression of the K8/K18 pair as observed in the heart failure mouse model, and the activity of the upstream factor TNF was also suppressed compared to controls, whereas the expression of the K5/K14 and K1/K10 pairs was significantly upregulated (Figures 3D,3E, and S10C). This suggests that several different keratin pairs contribute to myocardial protection depending on the type of myocardial stress of different intensity and duration in humans. Finally, we analyzed genome-wide trans-regulations of the post-transcriptional lincRNAs in HCM. To identify potential pairs of lincRNAs functional in HCM and their trans-target genes, we calculated the correlation of gene expression between pairs of 272 lincRNAs and 1930 PCGs, both with significantly altered expression levels (Figures 3A and S9A). Among the pairs with significantly correlated expression levels (absolute value of correlation coefficient >0.8), 303 pairs were positively correlated whereas 22 pairs were negatively correlated (Figure 4A). Furthermore, more than half of the negatively correlated pairs (15 pairs) were between PCGs and LINC00881. LINC00881 is known to be expressed specifically in cardiomyocytes after differentiation, whereas it was significantly reduced in HCM (Figure 4B). Gene ontology analysis of the 15 PCGs that were negatively correlated with LINC00881 showed an enrichment of ontology terms related to myocyte development (Figures 4B and 4C). According to lncATLAS, a database of subcellular localization of lncRNAs, LINC00881 is relatively localized in the nucleus, suggesting that it functions in the nucleus (Figure 4D). Because the 15 negatively correlated PCGs are located on different chromosomes and not near to the locus of LINC00881 (chromosome 3), it is highly likely that LINC00881 regulates these genes through nuclear factors. Therefore, we applied oPOSSUM (ver.3.0) to the promoter region sequence data to search for nuclear factors that commonly regulate the 15 PCGs (Table S2). In addition, to narrow down the number of plausible factors with different approaches, we used Ingenuity Pathway Analysis (IPA) to search for nuclear factors that function upstream of the set of PCGs that vary in HCM, and 13 regulatory factors were commonly identified in these two analyses (Figure 4E). Of the three transcription factors with positive IPA Activation Zscore, E2F1 was shown to bind well to the promoter region of the target gene (Figures S11A–S11E). These results support the importance of E2F1 as a regulator of HCM traits and suggest that down-regulation of LINC00881 may be involved in the regulation of myocardial-specific pathways in HCM.

In this study, we demonstrated that lincRNA loci are significantly enriched within TAD regions and their localization is associated with tissue specificity at a genome-wide level. More specifically, PCG loci were enriched in TAD boundaries, whereas lincRNAs showed locus enrichment for TADs. lincRNAs enriched in TADs had higher tissue specificity than other lincRNAs or PCG, suggesting that TADs could be the basis for tissue-specific expression of lincRNAs. The increased tissue specificity associated with TAD localization complements the mechanisms of tissue-specific expression by promoters and epigenomes that have been investigated. Importantly, we validated the cis functions previously reported for individual lincRNA loci, such as ANRIL and ANRASSF1,, as whole gene type properties and demonstrated its significance. As shown in Figure S7C, lincRNAs expressed in multiple tissues have increased tissue specificity according to TAD localization, and these lincRNAs’ regulation or function can be associated with a sub-TAD state in which the structure within the TAD changes from tissue to tissue, whereas one-to-one recognition of transcription factors and regulatory sequences may be preferable for the highly specific lincRNA regulation. Often, the unknown function of lincRNAs has been inferred from the function of their neighboring genes because of the expression correlation between lincRNAs and their neighboring genes. We show at the whole genome level that expression correlation is higher for genes in the same TAD than for gene pairs at other locations regardless of intergenic distance, as has been suggested so far.,, Especially cis-association between lincRNAs and neighboring PCGs in TADs, which are highly tissue-specific genes, could suggest their functional relationship. Furthermore, among genes in TADs, lincRNAs may function as core TAD factors because lincRNAs are located farther from the boundary than PCGs. Recent study shows that lncRNA Eleanor activates ESR1 gene and neighboring genes in the same Eleanor-TAD in breast cancer cells. Such lincRNA cis-activation at TAD level may occur in a broad genomic region. PCGs were significantly enriched at TAD boundaries with low expression specificity; McArthur and Capra reported that house-keeping genes are abundant at the highly conserved TAD boundary across multiple cell types, which is consistent with our results. In addition to lincRNAs and PCGs, we also examined enrichment in TADs for various gene types registered in GENCODE and found that the rate of enrichment varies among gene types. As this study shows, there is potential to gain new understanding of function and its regulation in many gene types from the perspective of genome structure. Future reevaluation of various genomic elements from the perspective of higher order structure can be expected to yield new insights into genomic function and evolution. We further surmised that the identification of two groups of cis- and trans-regulated genes of lincRNAs would allow us to infer the mechanisms that produce differences in traits between samples. To demonstrate this, we analyzed data from human HCM to identify a group of genes associated with HCM-specific expression variation of lincRNAs and their upstream genes, which provided insight into the pathways relevant to HCM pathogenesis. HCM is often caused by mutations in sarcomere-related genes. However, there are many cases in which genetic abnormalities are not identified, and there is a need to understand the molecular mechanisms of the pathogenesis of HCM and to develop therapeutic methods. Here, based on our new approach and data from patients with multiple genetic backgrounds, we found elevated keratin gene expression and an HCM-specific transcriptional pathway by E2F1 with down-regulation of LINC00881. Ectopic expression of keratin in myocardium has been reported in studies of desmin-deficient mouse models of heart failure, although little is known about its relevance in human cardiac disease. Although only expression of the keratin K8/K18 pair was detected in the mouse model, the human HCM samples showed little expression of K8/K18 and elevated expression of the K5/K14 and K1/K10 pairs instead. In mice, K8/K18 is thought to function as an alternative cytoskeleton to compensate for desmin deficiency, whereas what the differences in keratin pairing mean in human stressed myocardium has not yet been investigated. Future work is needed to clarify its precise background mechanisms. In conclusion, this study demonstrated that, at a genome-wide level, the lincRNA loci were significantly localized to TADs compared to PCGs and that this localization was associated with the tissue-specificity of gene expression. Based on this, we proposed a new analysis framework for diseases/traits that are relevant to lincRNA expression, applied it to HCM data, and identified transcriptional regulations specific to HCM pathogenicity. This study will contribute to a fundamental understanding of the molecular basis underlying tissue-/disease-specific expression of lincRNAs.

It should be noted that there is some limitations of this study: it uses data from existing tissue-derived TADs as information on 3D structures, and it is not possible to refer to sub-TADs or local loop structures that are finer than information on a given TAD region. Although more TAD coordinates from primary human tissues are needed, there are often different algorithm-based tools used in Hi-C data analysis to identify TADs resulting in the introduction of considerable technical noise. This makes it challenging to ensure robustness, especially with multiple external TAD datasets. It is difficult to clearly identify molecular mechanisms by which lincRNA expression specificity is increased in TADs or that cause the expression correlation between lincRNA and coding gene pairs in TADs. Future accumulation of time-series expression data and experimental corroboration is awaited to elucidate those mechanisms. Furthermore, the relatively small sample size (n = 37) for the analysis of HCM-specific transcriptional states could potentially introduce bias, requiring further external validation and experimental corroboration.



Further information and requests for resources and data should be directed to and fulfilled by the Lead Contact, Tatsuhiko Tsunoda (tsunoda@bs.s.u-tokyo.ac.jp).

This study did not generate new unique materials.

GTEx expression data (gene transcripts per kilobase million (TPM)) used in this study was downloaded from the GTEx portal [https://gtexportal.org/home/datasets; accessed 12.20.2019, file “GTEx_Analysis_2017-06-05_v8_RNASeQCv1.1.9_gene_tpm.gct.gz”], which contained 17,382 postmortem samples. The editing level was called on the GTEx v8p release (study accession phs000424.v8.p2). Based on gene type classification of GENCODE Release 25 annotation file, expression values of 7345 lincRNA and 18358 protein-coding genes (PCGs) were extracted from GTEx data and clustering was performed for each gene types. For sample clustering, we performed a t-SNE and a Euclidean distance-based clustering with scikit-learn (v.0.22) using the following parameters: sklearn.manifold.TSNE(n_components = 2, init = ’random’, random_state = 101, method = ’barnes_hut’, n_iter = 1000, verbose = 2, perplexity = 40).fit_transform(sample_expression_data).

Clustering based on lincRNA expression patterns across tissued was performed using scikit-learn (v.0.22) agglomerative clustering method with default parameters. Median gene-level TPM by tissue was log2-transformed and clustered into three classes to avoid imbalances in the number of genes per cluster. Gene expressions based on the three clusters determined were plotted on a heatmap.

All gene loci were annotated with GENCODE Release 25 annotation file. The basic gene annotation file was converted to BED file format using a custom bash script. Neighboring genes of lincRNAs were identified by bedtools (v2.29.2) using the distance between gene bodies. Pearson's correlation coefficients for expression correlations were calculated between adjacent genes with the GTEx expression data.

To assess the specificity of gene expression across all tissues, an entropy-based specificity score was calculated using the Jensen-Shannon specificity variance metric (JSS DPM) in the tspex package. Using GTEx expression data, the score for each gene was quantified as a continuous value between 0 (ubiquitous expression) and 1 (tissue-specific expression).

We downloaded and used TAD coordinates on hg38 estimated from human Hi-C data from 3D Genome Browser [http://3dgenome.fsm.northwestern.edu/publications.html; accessed 02.09.2020, file “hg38”]. For 13 primary tissues (adrenal gland, aorta, bladder, small bowel, prefrontal cortex, hippocampus, liver, lung, ovary, pancreas, psoas muscle, spleen and left ventricle). For TADs from each of the primary tissues, we identified all genes that reside within TADs for lincRNA, protein coding, antisense, processed pseudogene, transcribed processed pseudogene, processed transcript, sense intronic, To be Experimentally Confirmed (TEC), unprocessed pseudogene, transcribed unprocessed pseudogene, unitary pseudogene, transcribed unitary pseudogene, polymorphic pseudogene, miRNA, snRNA, snoRNA, rRNA and miscRNA. Only those genes of which the entire gene body was contained within TAD regions were considered. We used a permutation scheme to obtain randomized distribution of gene positions for evaluating whether or not loci of each gene type are enriched within TADs: the gene loci were re-shuffled (10,000 permutations genome-wide) to count the number of genes in TADs and to estimate permutation p-values. Each permutation consists of the same TAD coordinates and gene sets, in which TAD regions were kept whereas the gene positions were randomly shuffled using the shuffle utility of bedtools (v.2.29.2). This was repeated for each of the 13 tissues. To obtain a list of common TAD regions, TAD coordinates which overlapped across all the 13 tissues were determined using the intersect utility of bedtools (v.2.29.2). The resulting list of 2,326 common regions were used as common TADs for the statistical framework described below. Genes contained in the common TADs were identified using the intersect utility of bedtools (v.2.29.2). To confirm the insulating effect of TADs on gene expression, we examined the moving average of the expression correlations according to the intergenic distance for gene pairs with different positions in relation to the TAD and the correlation of TAD expression levels for each chromosome. Pearson's correlation coefficients were calculated for gene expression values in GTEx, and the intergenic distances were calculated with GENCODE Release 25 annotation file. Moving averages were used to examine trends in expression correlations and smooth out short-term variation by calculating the average of the data over an interval of log10(bp+1) intergenic distance. Confidence intervals were calculated from standard error of correlation coefficients within the sliding window. To quantify the expression levels of each TAD region, GTEx TPM data were log2-transformed and standardized, and median expression levels were calculated among genes included in each TAD. Pearson's correlation coefficients were calculated for expression levels between TADs. For all lincRNAs and PCGs, in addition to the permutation tests and tissue specificity score calculations described above, a circular randomized permutation test was performed (10,000 permutations chromosome-wide) that preserves the size and number of gene regions plus their structure in genomic relationship to each other, using the regionR package (v.1.30.0) and R (v.3.6.2). Permutation tests and tissue specificity scores described above were performed for lincRNAs and PCGs under the following conditions: exclusion of HK genes, restriction to lincRNAs with experimentally proven functional roles, setting a threshold for expression levels (median TPM more than 0.5 at least in one tissue). In addition to the position patterns of lincRNAs and PCGs within and outside TADs, we compared the distribution of the shortest distances between genes and TAD boundaries in order to examine their distribution within TADs. The distance from each lincRNA and PCG within a TAD to the nearest TAD boundary was calculated using the distanceToNearest utility of regionR packages (v.1.30.0) and R(v.3.6.2).

Publicly available RNA-Seq datasets with accession code GSE130036, which contained 37 fastq files (myocardial tissues from 28 HCM patients and 9 healthy donors) were downloaded from the NCBI GEO database. Quality control of the reads was performed with FastQC (v.0.11.8) and MultiQC (v.1.9). The reads were trimmed with fastp (v.0.20.1) and aligned to the human hg38 reference (GENCODE Release 25 basic gene annotation) using STAR (v.2.7.3a) with default parameters. Quantification of gene expression was performed using RSEM (v.1.3.1) with default parameters. The raw count matrix containing all the samples was created using tximport package (v.1.14.2) and R (v.3.6.2). Normalization and differential expression analysis were carried out using DESeq2 package (v1.26.0) and R (v.3.6.2). Difference of gene expression was considered significant if their log2 fold change (HCM/control) was >0.58 or <−0.58 and its false discovery rate (FDR) was <0.05 upon correction for multiple testing with the Benjamini–Hochberg method.

For determining DEGs between myocardial tissues of HCM patients and healthy donors, after confirming that 41567 genes with observed background expression (TPM >0 in at least one sample) covered 20595 reference genes registered in PANTHER, GO term enrichment analysis was performed using the PANTHER Overrepresentation Test (release 02.24.2021) and the GO ontology database (release 02.01.2021).

We calculated median log2 fold change of genes within each TAD that had at least one differentially expressed lincRNA according to the common TAD coordinates, setting the log2 fold change of genes with FDR ≥0.05 to 0. TADs with median log2 fold change >0.58 or <-0.58 for DEGs within them were defined as having significant differences in expression levels. We performed GO analysis on the genes included in TADs with significant differences in expression levels according to three gene sets: all genes, up-regulated genes, and down-regulated genes.

Correlations of expression levels between differentially expressed lincRNAs and PCGs from all 37 samples were calculated. For PCG and lincRNA expression levels, we used normalized TPM values. These normalized values were obtained with RSEM (v.1.3.1). Only gene pairs with Pearson's correlation coefficient >0.8 were considered for further study.

To characterize the molecular function of LINC00881, we ascertained its subcellular localization using LncATLAS database (lncatlas.crg.eu). LncATLAS provides relative concentration index (RCI) to define and quantify RNA localization, which is defined as the log2-transformed ratio of FPKM (fragments per kilobase per million mapped) in the cytoplasm vs. nucleus.

To identify transcription factors predicted to be involved in regulatory function of LINC00881, we searched over-representation of predicted transcription factor binding sites (TFBS) in gene promoters and proximal 5′ regulatory regions (±5,000 bp from TSS) by Single Site Analysis (SSA) of oPOSSUM ver.3.0. The CORE collection of JASPAR database was used as the source of DNA binding profiles. A set of 15 coding genes negatively correlated with LINC00881 were input to SSA to obtain Z-scores and Fisher scores for each predicted transcription factor. To functionally filter the sequence-based predicted transcription factors, we selected nuclear factors involved in transcriptional regulation from those predicted to be upstream regulators (categories: “transcriptional regulator” and “ligand-dependent nuclear receptor”) by the use of Ingenuity Pathway Analysis (IPA, Qiagen). Z-Score (ZS), used to calculate the predicted activation state from IPA, was used to infer likely activation states of upstream regulators based on comparison with a model that assigns random regulation directions. Only those hits found in both analyses were considered for the overlap. To support the ability of the identified co-factor candidates to actually regulate their target genes, human ChIP-seq data registered in ChIP-Atlas, for 6 genes (OXO3, NKX2-5, E2F1, MYCN, MAX and NR3C1) were specified and the corresponding results were streamed to IGV (v.2.16.0) to obtain snapshots.

Student’s t-test was used to determine the percentage of each gene type in the TAD. Mann-Whitney-Wilcoxon test two-sided with Bonferroni correction was used to compare gene specificity between the two groups. ns: P> 0.05, ∗∗P< 0.01, ∗∗∗P< 0.001, ∗∗∗∗P< 10-4, Student’s t-test and Mann-Whitney-Wilcoxon test. P< 0.05 and FDR <0.05 were considered to be statistically significant. Unless otherwise stated, all statistical analyses were performed in Python 3.8.0 statistical environment, using the packages pandas (v.0.25.3), scikit-learn (v.0.22), scipy (v.1.5.2), numpy (v.1.17.4), statannot (v.0.2.3) and tspex (v.0.6.2) and all plots were generated in python using the packages seaborn (v.0.11.0) and matplotlib (v.3.1.2).