Analysis of lncRNA-miRNA-mRNA Interactions in Hyper-proliferative Human Pulmonary Arterial Smooth Muscle Cells

We previously reported enhanced proliferation of smooth muscle cells on the combined exposure of HIV proteins and cocaine leading to the development of HIV-pulmonary arterial hypertension. Here, we attempt to comprehensively understand the interactions between long noncoding RNAs (lncRNAs), mRNAs and micro-RNAs (miRNAs) to determine their role in smooth muscle hyperplasia. Differential expression of lncRNAs, mRNAs and miRNAs were obtained by microarray and small-RNA sequencing from HPASMCs treated with and without cocaine and/or HIV-Tat. LncRNA to mRNA associations were conjectured by analyzing their genomic proximity and by interrogating their association to vascular diseases and cancer co-expression patterns reported in the relevant databases. Neuro-active ligand receptor signaling, Ras signaling and PI3-Akt pathway were among the top pathways enriched in either differentially expressed mRNAs or mRNAs associated to lncRNAs. HPASMC with combined exposure to cocaine and Tat (C + T) vs control identified the following top lncRNA-mRNA pairs, ENST00000495536-HOXB13, T216482-CBL, ENST00000602736-GDF7, and, TCONS_00020413-RND1. Many of the down-regulated miRNAs in the HPASMCs treated with C + T were found to be anti-proliferative and targets of up-regulated lncRNAs targeting up-regulated mRNAs, including down-regulation of miR-185, -491 and up-regulation of corresponding ENST00000585387. Specific knock down of the selected lncRNAs highlighted the importance of non-coding RNAs in smooth muscle hyperplasia.

down-regulate the expression of bone morphogenetic protein receptors (BMPR), HIV-Tat demonstrated maximum effect in the presence of cocaine in human PASMCs 22 .
Multiple molecular mechanisms including pathways regulated by noncoding RNAs have been reported to contribute to the development of PAH 23,24 . The significance of noncoding RNAs can be realized by the facts that only 2% of the human genome codes for protein coding genes and 90% of the rest of 98% of the genome codes for these non-protein coding genes 25 . Non-coding RNAs such as micro RNAs (miRNAs) and long non-coding RNAs (lncRNAs) play critical roles in regulating the fundamental aspects of cellular homeostasis and tissue development 26,27 . miRNAs are generally short RNAs with 21-23 nucleotides (nts) that regulate gene expression post transcriptionally by binding to mRNAs in a sequence specific manner. We had earlier characterized a novel role for miR-216 and -301 in modulating HIV-protein Tat and cocaine mediated SMC hyperplasia via targeting BMPR-2 mRNA 10 . Our studies also revealed the role of miRNAs carried in extracellular vesicles derived from HIV-infected and morphine exposed macrophages in regulating proliferation of PASMCs 11 .
LncRNAs, novel class of non-coding RNAs that were recently discovered with the advancement of RNA-sequencing techniques are longer than 200 nts 28 . While the roles of miRNAs in PAH have been elaborately studied, functions of lncRNAs in smooth muscle dysfunction and PAH remain less explored. lncRNAs are reported to regulate gene expression by multiple mechanisms. Some function as decoy, guide or scaffold to a diverse set of transcription factors and chromatin silencing proteins and therefore regulate expression of protein coding genes at the transcriptional and post-transcriptional levels 28 . The lncRNAs sharing common miRNA bindings sites as target mRNA can also act as competitive endogenous RNA (ceRNA) by functioning as miRNA sponge or they can also act as precursors and encode miRNAs therefore regulating the levels and downstream functions of miRNAs 29 . Since multiple ways of interaction between miRNAs to lncRNAs, miRNA to mRNA, lncRNA to mRNA (cis and trans) have been reported to play crucial roles in determining the cellular fate, it becomes essential to study these interactions in an integrated fashion to comprehensively decipher the factors underpinning the pathology of PAH.
Here, we report the first integrated analysis of altered lncRNA, miRNA and mRNA expression in hyper-proliferative human primary PASMCs treated with cocaine and Tat uncovering the critical interactions between these molecules that may eventually help us not only in understanding the cocaine and Tat induced hyper-proliferation but also in general regulation of smooth muscle dysfunction associated with all types of pulmonary hypertension.

Results
Differential expression of lncRNAs in Cocaine and Tat protein treated HPASMCs. As shown in the schematic diagram in Fig. 1, microarrays were used to capture the expression profile of lncRNAs in Tat and cocaine treated HPASMCs. This analysis revealed a number of lncRNAs with significant differential expression ( Fig. 2A). A total of 1,291 lncRNAs were significantly up-regulated and a total 2,155 were down-regulated in the combined treatment of cocaine and Tat (C + T) compared to control. LncRNAs have been classified into many different groups based on their genomic location with relevance to the nearby protein coding gene sequences 30 which helps to predict some of their possible mode of action in regulating other cellular RNA molecules. Based on this classification, significantly differentially expressed lncRNAs in the C + T group compared to the untreated control were found to be 72% (2,509) intergenic, 12% (415) natural antisense and 9% (321) intronic antisense regions (Fig. 2B) in the genome. These lncRNAs ranged from 100 to 43,681 nts in size (Fig. 2C) and were distributed across all chromosomes (Fig. 2D). In-silico functional association between differentially expressed lncRNAs and mRNAs in cocaine and tat treated HpAsMCs. We looked for associations between the significantly differentially expressed mRNAs and lncRNAs in the relevant databases mentioned in Methods, to find out the reported association between them along with their known functional significance. We investigated the association of all significantly altered lncRNA and mRNAs in the C + T group compared to control using the LncDisease database and found that 37 significantly dys-regulated lncRNAs, were related to multiple differentially expressed mRNAs associated with cardiovascular diseases. The up-regulated lncRNA ENST00000495536 was one of the top candidates in the list and was associated with down-regulated HOXB13 mRNA in cocaine + Tat group (Fig. 4A).
We next performed a correlation analysis of the significantly differentially expressed lncRNA and mRNA in the C + T group in conjunction with the lncRNA Cancer database. Significantly differentially expressed lncRNAs were associated with significantly differentially expressed mRNAs that showed a differential expression pattern www.nature.com/scientificreports www.nature.com/scientificreports/ up-regulated lncRNA T151531. Figure 4A, lists some of the top lncRNA and mRNA associations retrieved from this analysis and suggests underlining interactions between these RNAs in regulating critical proliferation signals in HPASMC.
Given that lncRNAs are known to regulate the expression of proximal protein coding genes in cis regulatory mechanisms 33 , we explored the lncRNA and mRNA microarray data to identify lncRNAs that may regulate in cis nearby mRNAs. A functional analysis of these mRNAs using the DAVID annotation tool revealed their association to multiple proliferation related pathways such as Hippo, PI3-Akt and cAMP signaling (Fig. 4B). A list of, long intergenic noncoding RNAs (lincRNAs) or antisense lncRNAs, along with their associated mRNA was prepared for the significantly differentially expressed lincRNAs and antisense lncRNAs in the C + T group based on their genomic proximity (<300 kbp from each other). A total of 345 pairs of lincRNA -mRNA associations were found. The lncRNA TCONS_00020413 was one of the top candidates in this list and was up-regulated in C + T compared to control, while expression of its nearby gene, RND1, was significantly down regulated (Fig. 4A). RND1 is a tumor suppressor protein that belongs to the Rho family GTPases. A biological functional analysis of these significantly differentially expressed mRNAs, in the lincRNA/mRNA pairs, done using ingenuity pathway analysis (IPA), listed cancer as one of the top diseases, with 208 mRNAs out of 345 being associated with it. GDF7 was one of the top candidate genes among the cancer related molecules. It was associated with three lncR-NAs namely, T185733, NR_046836 and ENST00000602736 (Fig. 4A). The analysis of significantly differentially expressed antisense lncRNA and significantly differentially expressed mRNA, that were in close proximity to each other, yielded 21 lncRNA/mRNA pairs. The up-regulated antisense lncRNA T041024 in C + T group was in proximity to the down-regulated ALOX5 gene implicated in hypertension.
In order to validate the microarray expression data of the identified lncRNA/mRNA pairs, we performed qRT-PCR on the top ranked lncRNAs that had a raw intensity above 100 (raw intensity is a surrogate measure of gene expression in microarrays) and their associated mRNAs from each of the three data sets mentioned above (Fig. 4A). From the lncRNA/mRNA associations established using the vascular disease database, we confirmed www.nature.com/scientificreports www.nature.com/scientificreports/ the up-regulation of the lncRNA, ENST00000495536, which was 2.4-fold up-regulated in cocaine + Tat treated HPASMCs compared to untreated cells. This lncRNA was associated with the HOXB13 gene, which was 2.1-fold down-regulated in cocaine + Tat. These alterations in the cocaine + Tat treated cells were also significant when compared to HPASMCs treated with only cocaine or Tat mono-treatments. Furthermore, mono-treatment of cocaine or Tat resulted in no significant change in either ENST00000495536 or HOXB13 mRNA levels compared to un-treated controls (Fig. 4C). We also observed that the lncRNA and mRNA pair NR_108042-BDKR1 was up or down-regulated ~2 -fold, respectively, in the cocaine + Tat group compared to control. From the lncRNA and mRNA associations obtained by cancer co expression analysis, we validated by qRT-PCR, the differential expression of lncRNA T216482 (2 fold down-regulated) and its associated mRNA, CBL (1.9 fold down-regulated) in C + T treated cells compared to control (Fig. 4C). These changes were also similarly significant when C + T treated cells were compared to cells with the mono-treatment of cocaine or Tat.
We also validated by qRT-PCR, the expression of lincRNA TCONS_00020413 identified from the genome proximity analysis (Fig. 4A) and its associated gene, RND1. This lncRNA was significantly up-regulated by 2-fold in C + T treated HPASMC compared to control, whereas mono-treated cells with either cocaine or Tat resulted in only a non-significant 1.3 fold and 1.5 fold up-regulation, respectively (Fig. 4D). Its associated mRNA, RND1, was 2-fold down-regulated in cells treated with the cocaine + Tat (Fig. 4D). Further, we validated the mRNA expression of GDF7 (BMP12) and the three lincRNAs associated with it that was identified by their close genomic proximity. Validation by qRT-PCR established a 2-fold up-regulation of GDF7 while its associated lncRNA T185733 was 2.8-fold up-regulated and the other two lncRNAs, NR_046836 and ENST00000602736, were 1.7-fold and 1.8-fold down-regulated respectively, in cocaine + Tat compared to control (Fig. 4D). One of the lncRNA -mRNA pairs that was identified from the microarray analysis and selected for validation by qRT-PCR, NR_027022 -GOLTIA, showed no significant difference in expression in cocaine + Tat treated cells when compared to untreated control (Fig. 4D).  www.nature.com/scientificreports www.nature.com/scientificreports/

Dys-regulation of miRNAs in cocaine and tat protein treated HpAsMCs.
To determine the miR-NAs that are involved in Tat and cocaine mediated hyper-proliferation of smooth muscle cells, we performed small RNA sequencing on total RNA isolated from HPASMC treated with cocaine and Tat for 48 h. Analysis of significantly differentially expressed miRNAs revealed 78 miRNAs in the cocaine + Tat group, 48 in the Tat only group and 60 in the cocaine only group, when compared with untreated control cells (Fig. 5A). In order to determine the genes that are regulated by these miRNAs and ascertain their functional relevance, we used the IPA miRNA target filtering tool to identify all the known down-stream target genes of these miRNAs and performed a biological functional analysis on them using DAVID (Fig. 6B). We next filtered these genes to only include those genes that were significantly differentially expressed in our microarray data and performed a biological functional analysis on them (Fig. 5C). The above analysis was performed for miRNAs and mRNAs that were significantly differentially expressed in cocaine + Tat group compared to control. The functional analysis results from DAVID on all down-stream mRNA targets of the significantly differentially expressed miRNAs suggest that these miR-NAs may regulate mRNAs related to cancer and proliferation related pathways such as PI3-Akt, Ras, MAPK (Fig. 5B). Furthermore, the functional analysis results from DAVID on the significantly differentially expressed down-stream mRNA targets of the significantly differentially expressed miRNAs revealed associations to various proliferation related pathways such as PI3-Akt, Ras and calcium signaling (Fig. 5C).
In addition, a combined functional analysis of the significantly differentially expressed mRNA and miRNAs in cocaine + Tat treated HPASMCs using IPA highlighted several PAH related factors that support smooth muscle proliferation (Fig. 5D). Among the identified targets were the down-regulated genes, BMP7, CBL, TNFSF9, ITLN1, Nos3, TNRC6B, HDAC, and the up-regulated genes, CALCA, NOX1, KALRN, CXCL16, MYOCD and NOS1, which are known to induce smooth muscle cell proliferation. The analysis also revealed the significantly differentially expressed miRNAs in our data, miR-221, -10a, -122b, -21 and let-7, which are reported to activate proliferative signals in smooth muscle cells (Fig. 5D). Overall, the dysregulated miRNA and mRNA profile in cocaine + Tat smooth muscle cells shows that they may play a critical role in cocaine + Tat induced hyper-proliferation of HPASMC. www.nature.com/scientificreports www.nature.com/scientificreports/ Integrated analysis of significantly dys-regulated lncRNA, miRNA and mRNA in HPASMC exposed to combined treatment of cocaine and tat. Long noncoding RNAs are known to function as competing endogenous RNAs (ceRNAs) and affect mRNA levels by sequestering common miRNAs that target both the lncRNA and mRNAs 34 . We explored this function of ceRNA by analyzing the association of the up-regulated lncRNAs and www.nature.com/scientificreports www.nature.com/scientificreports/ mRNAs in our microarray data to the down-regulated miRNAs from our small RNAseq data for the HPASMCs exposed to the combined treatment of cocaine and Tat. Since lncRNAs with high expression levels are reported to function as sponge RNAs much effectively than lncRNAs with comparatively low expression 35 , we selected those lncR-NAs that were significantly up-regulated with raw signal intensity greater than 100, for this analysis. For reasons of annotation, only those lncRNAs that were listed in any one or more of the public databases, GENCODE and Ensemble or described by Cabili et al. 36 , were selected. We used LncBase v.2 available in DIANA Tools to retrieve information on miRNA binding to specific lncRNAs. The LncBase v.2 database holds information for over 70 k miRNA to lncRNA interactions, supported by either direct or indirect experiments derived from manually curated publications, and the analysis of 153 Argonaut-crosslinking immunoprecipitation (AGO CLIP)-Seq libraries 37 . We used IPA to ascertain the experimentally validated miRNA to mRNA interactions.
Out of the selected up-regulated lncRNAs, 10 of them had one or more binding sites for 12 of the down-regulated miRNAs. These 12 miRNAs targeted 134 of the up-regulated mRNAs. Using this information, we generated a putative molecular interactive network describing how lncRNAs acting as ceRNAs could potentially regulate mRNA expression by modulating the level of its up-stream miRNA regulator(s) ( Table 3). This analysis indicates that some of the up-regulated lncRNAs may function as sponge or ceRNAs of the down-regulated anti-proliferative miRNAs in hyper-proliferative smooth muscle cells.
Of the analyzed up-regulated lncRNAs, TCONS_00001909 had 22, 7, and 6 predicted binding sites for the down-regulated miRNAs, miR-376b, -193a and -370, respectively. This connection can potentially explain the role of this lncRNA in the up-regulation of mRNAs such as ZNF74, SCN11A, C10orf25, MAGEA3/MAGEA6, UBXN10 and TSNARE1 that are regulated by these miRNAs ( Table 3). The lncRNA, TCONS_00028198 contains 15, 13 and 21 binding sites for the down-regulated miRNAs, miR-376b, -33a, and -185 respectively, which in turn is known to regulate several of the up-regulated mRNAs implicated in proliferation such as, VAT1L, UBE2QL1, SSX7, SLAMF6, S100A7A, RNASE2, PLN, PCDHA10, PABPC1L2A, CXCL16, LRRC38, FXYD4, ETV3L, ZNF74, SCN11A and C10orf25 (Table 3). The lncRNA, ENST00000585387 has two predicted binding sites each for the miRNAs, miR-491 and -185 and thus may act as a ceRNA for these miRNAs which are known to regulate the up-regulated genes, NOS1, KCNIP2, IGF2, ELF5, EFNA3, CHAD, C17orf107 and BAK1. Overall, evidence from our bioinformatic analysis suggests that many lncRNAs that are significantly up-regulated have the potential to sequester many anti-proliferative miRNAs in HPASMCs exposed to combined treatment of cocaine and Tat and thereby regulate the expression of their downstream target genes.
Our in silico analysis using down-regulated mRNA/lncRNA with up-regulated miRNAs also identified an inter-regulatory relationship between them. We identified 21 down-regulated lncRNAs that may potentially target 7 up-regulated miRNAs in cells treated with combined C + T treatment when compared with untreated controls. Based on prediction and experimentally validated findings, these miRNAs had a potential to target 131 down-regulated mRNAs in C + T vs control comparison (Table 4).
Next, we validated the expression of the lncRNAs, TCONS_00001909, TCONS_00028198 and ENST00000585387 along with four of their corresponding predicted miRNA targets: miR-376, -193, -185 and -491 (Fig. 6A), using qRT-PCR. These three lncRNAs were up-regulated in qRT-PCR by 2, 1.8 and 3.6 -fold, respectively, which was in concordance with the microarray data. The expression of these selected CeRNAs in C + T treated cells was significantly up regulated not only when compared to untreated cells, but also when compared to mono-treated cells with cocaine or Tat. All four target miRNAs validated by qRT-PCR were down regulated in cocaine + Tat treated cells compared to other treatments (Fig. 6B).

Involvement of lncRNAs ENST00000585387 and ENST00000495536 in cocaine and
Tat mediated smooth muscle hyperplasia. Next, we selected the lncRNA ENST00000495536 that was associated with the gene, HOXB13 and the lncRNA ENST00000585387 premised to act as a ceRNA for miRNAs, miR-491, -185, for examining of their role in smooth muscle proliferation. We utilized gapmeR based LNAs to knockdown the expression of these two lncRNAs in HPASMCs. As shown in Fig. 7A, a 50% reduction in the levels of both lncRNA was observed, with a 40 nM concentration of antisense oligonucleotides (ASO) compared to HPASMCs transfected with scrambled gapmeRs (Fig. 7A). Both MTS and Cyquant proliferation assays revealed that the knockdown of either, ENST00000495536 or ENST00000585387, results in the inhibition of cocaine and www.nature.com/scientificreports www.nature.com/scientificreports/ Tat induced proliferation of HPASMCs compared to un-transfected or, scrambled gapmeR transfected C + T treated cells (Fig. 7E,F). These changes in proliferation were also clearly visible under the microscope as shown by the representative phase contrast images in Supplementary Fig. S1. Concomitant to these findings we observed that the knockdown of lncRNA 387 or 536 also prevented the C + T mediated down-regulation of miR-185 and -491 or HOXB13 mRNA (Fig. 7B,C). Together, these results suggest the role of lncRNAs in the cocaine and Tat mediated hyper proliferation of HPASMCs.

Discussion
Dys-regulated expression of miRNAs, mRNAs and lncRNAs has been reported to play a critical role in the progression of PAH 38,39 . However, an integrated analysis of the potential interactions between these different types of RNA molecules is only beginning to emerge 34,40 . In this study, we predicted the inter-regulatory relationship between lncRNAs, mRNAs and miRNAs and identified their role in hyper-proliferation of pulmonary smooth muscle cells. We utilized high throughput techniques such as microarray and small RNAseq to quantify the expression levels of mRNAs, lncRNAs and miRNAs. Our analysis of dys-regulated mRNAs and lncRNAs, as described above, revealed a number of lncRNA -mRNAs associations involved in regulating hyperplasia of  www.nature.com/scientificreports www.nature.com/scientificreports/ HPASMCs on exposure to cocaine and Tat protein. The KEGG pathway analysis of the significantly dys-regulated mRNAs identified neuroactive ligand receptor signaling and Ras signaling among the top signaling pathways. The combined analysis of dys-regulated miRNAs, mRNAs and lncRNAs revealed the potential of lncRNAs to act as ceRNAs/sponges for many anti-proliferative miRNAs leading to up-regulation of pro-proliferative mRNAs in cocaine and Tat treated HPASMC. Overall, these analyses provided clues to the possible cis and trans acting roles of various up-regulated lncRNAs in hyper-proliferative smooth muscle cells by inter-regulating, either directly or indirectly, transcriptional or post-transcriptional expression of downstream proliferative genes.
Recent research on lncRNAs is proving to be a key to unlocking many underlying molecular mechanisms of development and pathology of diseases 41 . A vast number of lncRNAs have been identified that still remain to be characterized. Adding to the complexity of understanding lncRNAs is their versatile nature of cellular regulation, making them both intriguing and challenging to be characterized 42 . Lately, many studies have identified multiple lncRNAs implicated in PAH 39,[43][44][45][46] , however, an integrated study considering miRNAs, mRNAs and lncRNAs with relevance to HPASMC hyper-proliferation and PAH have been lacking.
The pathway analysis we performed on dys-regulated mRNAs identified many factors playing active role in neuroactive ligand receptors and ras signaling. It is interesting to note that genes such as C5AR1, CHRNA4 and HTR1A from the neuroactive ligand receptor pathway have been reported in association with PAH, lung cancer, COPD and idiopathic pulmonary fibrosis [47][48][49] . On the other hand, Ras signaling molecules such as FGFR2, PLA2G2A and PLA2G5 have been shown to be involved in vascular remodeling, idiopathic pulmonary fibrosis  www.nature.com/scientificreports www.nature.com/scientificreports/ (IPF) and chronic thromboembolic pulmonary hypertension (CTEPH) [50][51][52] . Analyzing lncRNAs associated with molecules in these highly enriched/activated pathways in cocaine and Tat induced HPASMC can help us understand their mechanistic roles in these diseases.
One of the top up-regulated vascular disease related lncRNAs in HPASMCS treated with both cocaine and Tat was ENST00000495536. It is classified as an intergenic lncRNA located around 12,000 bp, in the antisense strand, from the HOXB13 gene, PRAC2 small nuclear protein and miRNA 3185. The gene, PRAC2, was not expressed in our mRNA microarray data, however the expression of HOXB13 was significantly down-regulated in cocaine and Tat treated cells in the microarray data, and was confirmed by qRT-PCR. Loss of HOXB13 in colorectal tumor cells has been associated with the increase in proliferation 53,54 . The knock down of HOXB13 associated lncRNA, ENST00000495536 in HPASMC prevented the cocaine and Tat mediated decrease in HOXB13 mRNA expression with the corresponding inhibition of smooth muscle hyper-proliferation. The possibility of an association www.nature.com/scientificreports www.nature.com/scientificreports/ between the lncRNA, ENST00000495536 and the gene, HOXB13 is further strengthened by the observation that in most cases, the lincRNAs function in trans-while antisense lncRNAs function in cis, to modulate gene expression 55 . However, the mechanistic association of HOXB13 and its nearby lncRNA, ENST00000495536 in HPASMC and PAH needs to be further established.
Our analysis identified cancer as one of the top diseases associated with the significantly differentially expressed lncRNAs and mRNAs and revealed the integrated link between the gene/lncRNA pairs: CBL/T216482 and GDF7 (BMP12)/T185733, prominently associated in cancer pathways. The expression of both the lncRNA, T216482 and its associated mRNA, CBL were down-regulated in cocaine and Tat treated HPASMCs. The CBL gene has been reported to be decreased in primary bone tumors, and ectopic CBL expression has been reported to reduce bone tumorigenesis by promoting tyrosine kinase receptor degradation 56 . Enhanced CBL activity is known to result in the down-regulation of EGFR expression and inhibition of proliferation in colon tumor cells 57 . GDF7 also known as BMP12 is reported to be involved in tissue regeneration 58,59 . It's in proximity to not only lncRNA T185733 but also to two other lncRNAs, NR_046836 and ENST00000602736. Furthermore, RND1 located up-stream of the up-regulated intergenic lncRNA, TCONS_00020413 is known to act as a tumor suppressor by attenuating Ras/MAPK signaling 60 . Its expression was reduced in response to the combined treatment of cocaine and Tat.
Analysis of significantly dys-regulated miRNAs and their mRNA targets in cocaine and Tat treated cells revealed many up-regulated miRNAs including miR-133 and miR-125 that target mRNAs involved in proliferative signaling. Of these two miRNAs, miR-133 is known to play role in vascular stress, remodeling and cell survival 61 . whereas miR-125 has been reported to play a role in pulmonary hypertension with increased expression observed in the lungs of hypoxic animals 62 . LncRNAs are known to regulate gene expression by sponging miRNAs that regulate mRNAs 63 . Our combined bioinformatics analysis of the up-regulated lncRNA and mRNAs, along with the down-regulated miRNAs, in cocaine and Tat cells revealed a number of lncRNA that possessed binding site for these down-regulated anti-proliferative miRNAs. We found that the lncRNA, TCONS_00028198 contained biding sites for miRNAs, miR-376b, -33a, and -185, and the lncRNA, ENST00000585387 contained bind sites for miRNA, miR-185 and -491, making them potential regulators of these miRNAs. Of these miRNAs, miR-376b and miR-33a have been shown to be significantly down-regulated in PH 64,65 and miR-185 has been shown to mediate lung epithelial cell death post oxidative stress 66 . In addition, multiple studies suggest miR-185 and miR-491 as suppressors of cell-proliferation including vascular smooth muscle cells [67][68][69][70] . We observed that the knock down of lncRNA 397 results in up-regulation of miR-185 and miR-491 in HPASMCs exposed to the combined treatment of cocaine and Tat and this correlated with the decrease in C + T mediated hyper-proliferation of cells. Another significantly up-regulated lncRNA, TCONS_00001909 contained potential binding sites for miRNA, miR-193a, implicated in PAH 71 and miR-370, recently identified as a tumor suppressor 72 .
Many of the up-regulated mRNA targets of these miRNAs have not yet been characterized in either smooth muscle proliferation or in PAH however some of the mRNAs namely, CXCL16 and LRRC38 are known to play a positive role in proliferation 73,74 . The up-regulated gene, IGF2, which is a predicted target of miR-185, has been reported to be up-regulated in cells having a hypoxic environment 75 . Furthermore, it is intriguing to note that some of the alterations in the expression of lncRNAs, mRNAs and miRNAs may be attributed to a regulatory mechanism(s) operative at the DNA level considering that we observed the association between the down-regulated histone deacetylase-4 (HDAC4) mRNA and up-regulated miR-10 76 in C + T treated cells. It's plausible that miRNA dependent changes in the expression of HDAC proteins, known to play a critical role in the chromatin condensation and transcriptional repression 77,78 may in turn be involved in regulating the expression of lncRNA,miRNA and mRNA levels.
Long noncoding RNAs such as NEAT1, NRON, lincRNA-p21 and PANDA are known to modulate HIV infection and replication [79][80][81] . However to the best of our knowledge, this is the first attempt to report the effect of viral proteins in the interplay between lncRNAs, mRNAs and miRNAs. Although the functions of lncRNAs characterized here are based on their primary sequence, there are reports that suggest that the tertiary structure of lncR-NAs play a major role in deciding lncRNA function 82 . Therefore, further functional evaluation of the identified lncRNA, miRNA and mRNA associations is needed for a comprehensive mechanistic understanding of their role in smooth muscle hyper-proliferation and pulmonary vascular remodeling.

Materials and Methods
Cell culture and RNA isolation. Primary human pulmonary arterial smooth muscle cells (HPASMC) were grown on smooth muscle cell media (SMCM) with growth factors, 2% fetal bovine serum and penicillin/streptomycin until 70% confluency on 6 well plate and made quiescent for 48 h with serum free media and then treated with either cocaine (C) at 1 µM final concentration or with Tat (T) protein at 25 ng/mL concentration or with both (C + T) for 12, 24 and 48 h based on our previous findings 12,21 . Cells were harvested in TriZol reagent and the total RNA was isolated as per the manufacturer recommended protocol. Considering that HIV-PAH is more prevalent in males particularly among intravenous drug users 1,83 pulmonary arterial smooth muscle cells used for this study were from males.
Microarray analysis for lncRNA and mRNA expression. The total RNA isolated from HPASMC were quantified by the NanoDrop ND-1000 and RNA integrity was assessed by standard denaturing agarose gel electrophoresis. For microarray analysis, Agilent Array platform was employed. The sample preparation and microarray hybridization were performed based on the manufacturer's standard protocols with minor modifications. Each sample was amplified and transcribed into fluorescent cRNA along the entire length of the transcripts without 3' bias utilizing a random priming method (Array Star Flash RNA Labeling Kit, Arraystar Inc). The labeled www.nature.com/scientificreports www.nature.com/scientificreports/ cRNAs were hybridized onto the Human LncRNA Array v4.0 (8 × 60 K, Arraystar). After washing the slides, the arrays were scanned by the Agilent Scanner G2505C.
Agilent Feature Extraction software (version 11.0.1.1) was used to analyze the acquired array images. Quantile normalization and subsequent data processing were performed using GeneSpring GX v12.1 software package (Agilent Technologies). The statistical significance of the differentially expressed LncRNAs and mRNAs between cocaine and Tat (C + T) group and control group were calculated using the Student's two sample t-test. The cutoff for significance was set at an absolute fold-change ≥1.5 and p-value ≤0.05. Biological functions of the significantly differentially expressed genes were obtained from the KEGG (Kyoto Encyclopedia of Genes and Genomes; http://www.genome.jp/kegg) and GO (Gene Ontology; http://www.geneontology.org) databases/online-software. The data were also analyzed through the use of IPA (QIAGEN Inc., https://www.qiagenbioinformatics.com/products/ingenuitypathway-analysis) for the analysis of associated pathways. The statistical significance of the association of genes with a biological function or pathway was calculated using the right-tailed Fisher's Exact Test.
small RNA sequencing. The global unbiased miRNA profiles of HPASMCs treated with cocaine and/ or Tat were interrogated using illumina's small RNA sequencing technology (TruSeq Small RNA Sample Prep Kit). Samples were prepared in triplicates (n = 3, for each treatment group). The small RNA-seq samples were sequenced for 50-cycle single end reads using a HiSeq2500 Sequencing System (Illumina, San Diego, CA). The 12 samples were multiplexed in a single lane. The sequence yield of the lane was approximately 375 million reads providing enough sequencing depth for each small RNA-Seq sample.
The small RNA-Seq reads with non-canonical letters (e.g. N) were first removed from the samples. Adaptors were clipped from the remaining reads. Resulting reads that were shorter than 17 bp were discarded. The remaining reads were mapped to the human genome (GRCh38.rel77) using Bowtie 84 . The read abundance estimates of all known human miRNAs (mirBase v21) were computed using miRDeep2 85 .
The miRNA expression counts obtained from this process were analyzed using a negative binomial generalized linear model (NB-GLM) from the edgeR package 86 to identify statistically significant differentially expressed miRNAs between the control and the three treatment conditions, cocaine (C), Tat (T) and combined treatment of cocaine and Tat (C + T). The edgeR package implements advance empirical Bayes methods to estimate feature-specific biological variation under minimal levels of biological replication. Finally, Hierarchical Clustering (Euclidean distance and Ward's linkage) was performed to show the distinguishing clusters of miRNA expression patterns among different sample groups.
Integrative analysis of dys-regulated non coding RNAs and mRNAs. Differentially expressed mRNAs, miRNAs and lncRNAs that were more than 1.5 fold up or down regulated between treatments groups (C + T vs. Con) with a p-value less than 0.05 were considered for integrated analysis. The association between significantly differentially expressed miRNAs and mRNAs were explored using the IPA software. The analysis was performed, first by only considering experimentally validated miRNA/mRNA interactions without including microarray mRNA data (Fig. 5B); and then by including microarray mRNA data to also consider experimentally validated and highly predictive miRNA/mRNA interactions (Fig. 5C). The IPA software uses information on experimentally validated miRNA/mRNA interactions from TarBase and miRecords. It uses information on predicted miRNA/mRNA interactions from TargetScan. IPA also maintains a repository of miRNA related information from peer-reviewed literature that it uses for its analysis. The regulatory relationships between differentially expressed miRNAs and lncRNAs were obtained from LncBase v2.0 using the DIANA Tools software 37 . Only experimentally validated interactions were considered. Both IPA software and DAVID (Database for Annotation, Visualization and Integrated Discovery) software 87 were used to perform pathway enrichment analysis.
Significantly differentially expressed lncRNAs were also annotated for biological functions and diseases using several databases. In particular, the LncRNADisease database 88 provided a source of experimentally validated and predicted lncRNA-disease associations, including a separate list of predicted lncRNA association with cardiovascular function and disease, and a list of lincRNA disease associations predicted by co-expression analysis. This database also provided a curated list of experimentally supported lncRNA interactions with mRNA. We used data from the lncRNA-Cancer gene co-expression network (lnCaNet) 89 to identify genes that are putatively co-expressed with lncRNAs based on expression correlation. Information from the Lnc2Cancer database 90 was used to obtain experimentally validated associations between lncRNA and human cancer.
Quantitative Real time (qRt)-pCR. The cDNA was prepared using RT 2 First strand synthesis kit (Cat#330404, Qiagen) as per the manufacturer's protocol. qRT-PCR was performed using RT 2 SYBR Green ROX qPCR Master Mix kit (Cat#330522, Qiagen). The qRT-PCR primers for the selected mRNA and lncRNAs were either custom designed using PrimerQuest tool and obtained from Integrated DNA technologies or ordered form Qiagen. Tables S1 and S2 (Supplementary File II) lists details of all the primers used in this study.
GapmeR antisense oligonucleotide transfections. Antisense oligonucleotide (ASO) based gapmeRs were designed against lncRNA using Qiagen, gapmer designer tool. The top ranked ASO were ordered along with positive and negative controls. HPASMC were reverse transfected while seeding in 6 well plate at 2 × 10 5 cells/well with specific gapmeR ASO at 40 nM using HiPerfect transfection reagent (Cat#301704; Qiagen). After 24 h of transfection, the cells were serum starved for 24 h and then treated for 48 h with cocaine and Tat. RNA isolation and qRT-PCR were performed as described above to evaluate the knockdown of specific lncRNAs. For proliferation assay cells were reverse transfected in a 96 well plate with a cell density of 10 4 cells/well. Proliferation assay was performed as described earlier 9 . Proliferation of the cells was quantified using CellTiter 96 ® AQueous (2019) 9:10533 | https://doi.org/10.1038/s41598-019-46981-4 www.nature.com/scientificreports www.nature.com/scientificreports/ One Solution Cell Proliferation Assay (MTS) (Promega,G3582) and CyQUANT ® Cell Proliferation Assay Kit (Invitrogen, C7026) as recommended by the manufacturers. statistical analysis. The statistical analyses of the qRT-PCR data as shown in Figs 4C,D, 6A,B and 7A-D were carried out using the normalized data obtained relative to control. We used Kruskal Wallis test to assess if differences exist among the treatment groups followed by Dunn's test for pairwise comparisons. The multiple testing adjustments were carried out using Benjamini and Hochberg's method. Since our interest was also to assess the treatment effects of cocaine and Tat separately as well as combined cocaine and Tat (as a single group) with respect to untreated control we later used one factor analyses model (Kruskal Wallis test) where the Factor is the treatment effects. In this way, we were able to assess the individual treatment effects as well as the synergistic effect of the combined treatments. We did not design the study as a two factor model because we were not assessing the interaction effect of the treatments. The statistical analyses for the proliferation assay data as shown in Fig. 7E,F was carried out similarly on the background corrected raw readings. After the completion of the statistical analyses, for the better visual clarity of the outcomes, the data in the figures were presented in the relative scale as compared to control by calculating the fold changes. Since the control group was the reference for the fold change calculation, its value was set to 1 in the bar plots and hence lack the error bars. The test results were considered significant if the multiple tests adjusted P values were ≤ 0.05.

Data Availability
The datasets generated during and/or analyzed during the current study are available from the corresponding author on reasonable request.