Long Non-Coding RNAs Might Regulate Phenotypic Switch of Vascular Smooth Muscle Cells Acting as ceRNA: Implications for In-Stent Restenosis

Coronary in-stent restenosis is a late complication of angioplasty. It is a multifactorial process that involves vascular smooth muscle cells (VSMCs), endothelial cells, and inflammatory and genetic factors. In this study, the transcriptomic landscape of VSMCs’ phenotypic switch process was assessed under stimuli resembling stent injury. Co-cultured contractile VSMCs and endothelial cells were exposed to a bare metal stent and platelet-derived growth factor (PDGF-BB) 20 ng/mL. Migratory capacity (wound healing assay), proliferative capacity, and cell cycle analysis of the VSMCs were performed. RNAseq analysis of contractile vs. proliferative VSMCs was performed. Gene differential expression (DE), identification of new long non-coding RNA candidates (lncRNAs), gene ontology (GO), and pathway enrichment (KEGG) were analyzed. A competing endogenous RNA network was constructed, and significant lncRNA–miRNA–mRNA axes were selected. VSMCs exposed to “stent injury” conditions showed morphologic changes, with proliferative and migratory capacities progressing from G0-G1 cell cycle phase to S and G2-M. RNAseq analysis showed DE of 1099, 509 and 64 differentially expressed mRNAs, lncRNAs, and miRNAs, respectively. GO analysis of DE genes showed significant enrichment in collagen and extracellular matrix organization, regulation of smooth muscle cell proliferation, and collagen biosynthetic process. The main upregulated nodes in the lncRNA-mediated ceRNA network were PVT1 and HIF1-AS2, with downregulation of ACTA2-AS1 and MIR663AHG. The PVT1 ceRNA axis appears to be an attractive target for in-stent restenosis diagnosis and treatment.


Introduction
Coronary in-stent restenosis (SR) continues to challenge the scientific community, despite the emergence of new biomedical materials [1] and its association with medications (drug eluting stents) (DES) [2]. SR is a multifactorial process that involves the active participation of vascular smooth muscle cells (VSMCs) and endothelial cells as well as inflammatory and genetic factors [3]. VSMCs have great plasticity, with the ability to modify their synthetic machinery and cyto-architecture in response to various stimuli [4].
VSMCs are typically found in a quiescent state as well-differentiated cells with a functional contractile apparatus, requiring little protein turnover [5]. However, in pathological situations such as stent induced stress their molecular mechanisms are triggered, and VSMCs thereby transform into cells with a high proliferative and migratory capacity that secrete extracellular matrix proteins [6]. This phenotypic switch is the therapeutic target of DES via drugs that inhibit the cell cycle to control the expansion of VSMCs and the formation of neo-intima [7]. genetic biomarker discovery to allow for early diagnosis of post-angioplasty coronary stent restenosis.

Phenotypic Change Induction by Stent and PDGF-BB
HUASMC grown under different conditions are shown in Figure 1. Differences in growth pattern were observed by inducing differentiation of the contractile phenotype, with elongated, spindle-shaped cells. HUASMC of the stent-induced injury model showed a reticular organization with clear nuclear definition, especially in the vicinity of the stent ( Figure 1A). Furthermore, the stent injury model significantly promoted the proliferation of HUASMC ( Figure 1B). Contractile cells were quiescent, most in G0-G1 cell cycle phases. Cells entered the S and G2-M phases under stent and PDGF stimuli ( Figure  1C).

Stent Injury Model Promotes Cells Migration and the Expression of Genetic Markers of Phenotypic Switch
Cell migration was evaluated by wound healing assay, which revealed that the migratory capacity of HUASMC increased in the stent injury group (Figure 2A,B). Genes associated with a contractile or synthetic/proliferative phenotype are plotted in Figure 2C [3,22]. In the stent-stimulated group, there was lower expression of contractile phenotyperelated genes such as smooth muscle α-actin (ACTA2), calponin (CNN1), and collagen type IV (COL4A1), whereas genes associated with the synthetic/proliferative phenotype, such as osteopontine (SPP1), matrix metalloproteinase 1 (MMP1), and fibronectine (FN1), were upregulated.

Stent Injury Model Promotes Cells Migration and the Expression of Genetic Markers of Phenotypic Switch
Cell migration was evaluated by wound healing assay, which revealed that the migratory capacity of HUASMC increased in the stent injury group (Figure 2A,B). Genes associated with a contractile or synthetic/proliferative phenotype are plotted in Figure 2C [3,22]. In the stent-stimulated group, there was lower expression of contractile phenotype-related genes such as smooth muscle α-actin (ACTA2), calponin (CNN1), and collagen type IV (COL4A1), whereas genes associated with the synthetic/proliferative phenotype, such as osteopontine (SPP1), matrix metalloproteinase 1 (MMP1), and fibronectine (FN1), were upregulated.

Figure 2.
A wound healing assay was performed to compare the migratory capacity in both con tions; selected images of both conditions are depicted in (A). Panel (B) shows the quantitative pression of this experiment. Contractile cells were mostly quiescent, while stent-induced prolife tive cells covered almost 60% of the scratch at 36 h. Panel (C) shows the expression of selected ge known to be markers of phenotypic switch; blue and red color intensity accounts for down-a upregulation, respectively. In contractile cells there was an overexpression of ACTA2, CNN1, a COL4A1; the same genes were downregulated in the stent-induced proliferative phenotype, wh FN1, MMP1, and SPP1 were uniformly upregulated. Cont_ = Replicates with contractile HUASM Stent_ = Replicates with stent-injured HUASMC; ** p < 0.01; *** p < 0.001.

Stent Injury Model Drives Activation of ncRNA Genes along with mRNA Downregulatio
The transcriptomic profile of the stent injury model cells revealed activation and su pression of both coding and non-coding genes (ncRNA) ( Figure 3A). The relative abu dance of ncRNA was superior to mRNA, representing 65.7 % of the mapped transcri ( Figure 3B). Using a threshold of absolute Fold Change ≥ 2 and an FDR ≤ 0.05, we iden fied 3517 differentially expressed (DE) mRNA and 2911 DE ncRNA. There was no diff ence in expression level among DE ncRNA, and DE mRNAs were mostly downregulat (Table 1)  . Panel (B) shows the quantitative expression of this experiment. Contractile cells were mostly quiescent, while stent-induced proliferative cells covered almost 60% of the scratch at 36 h. Panel (C) shows the expression of selected genes known to be markers of phenotypic switch; blue and red color intensity accounts for down-and upregulation, respectively. In contractile cells there was an overexpression of ACTA2, CNN1, and COL4A1; the same genes were downregulated in the stent-induced proliferative phenotype, while FN1, MMP1, and SPP1 were uniformly upregulated. Cont_ = Replicates with contractile HUASMC; Stent_ = Replicates with stent-injured HUASMC; ** p < 0.01; *** p < 0.001.

Stent Injury Model Drives Activation of ncRNA Genes along with mRNA Downregulation
The transcriptomic profile of the stent injury model cells revealed activation and suppression of both coding and non-coding genes (ncRNA) ( Figure 3A). The relative abundance of ncRNA was superior to mRNA, representing 65.7 % of the mapped transcripts ( Figure 3B). Using a threshold of absolute Fold Change ≥ 2 and an FDR ≤ 0.05, we identified 3517 differentially expressed (DE) mRNA and 2911 DE ncRNA. There was no difference in expression level among DE ncRNA, and DE mRNAs were mostly downregulated (Table 1).

Identification of Differentially Expressed Genes in HUASMC Exposed to PDGF and Stent
Differential expression analysis confirmed substantial differences in mRNA and ncRNA expression between the contractile and stimulated cells. Using more stringent criteria, FDR ≤ 0.01 and absolute log2FC ≥ 2 to declare significance and FPKM > 1 to confirm quantifiable expression, we identified 1099, 509, and 64 differentially expressed mRNA, lncRNA, and microRNA, respectively ( Figure 4). Table 2 shows that most of the protein coding genes (72.1%) were silenced during the phenotypic switch of HUASMCs, while lncRNAs were overexpressed in this process (59.3%). MicroRNAs did not show any differences under the experimental conditions. Expression pattern of all transcripts (A) and non-coding transcripts (ncRNA) (B). Labeled red dots indicate differential expression (≥ two-fold change and FDR ≤ 0.05); volcano plots reflect number, significance, and reliability of differentially expressed transcripts; red dots; and green dots indicate upregulation and downregulation, respectively. The x-axis represents the value of log2 Fold change and the y-axis represents the adjusted FDR.

Identification of Differentially Expressed Genes in HUASMC Exposed to PDGF and Stent
Differential expression analysis confirmed substantial differences in mRNA and ncRNA expression between the contractile and stimulated cells. Using more stringent criteria, FDR ≤ 0.01 and absolute log2FC ≥ 2 to declare significance and FPKM > 1 to confirm quantifiable expression, we identified 1099, 509, and 64 differentially expressed mRNA, lncRNA, and microRNA, respectively ( Figure 4). Table 2 shows that most of the protein coding genes (72.1%) were silenced during the phenotypic switch of HUASMCs, while lncRNAs were overexpressed in this process (59.3%). MicroRNAs did not show any differences under the experimental conditions.     Table 3 shows the top twenty dysregulated lncRNA during phenotypic switch induced by stent and PDGF. Among the top ten upregulated genes six have not previously been annotated, and there are two antisense transcripts and two long intervening noncoding RNAs (lincRNAs) as well. Meanwhile, in the downregulated group there is only one novel LncRNA along with four antisense, two lincRNAs, and two sense-intronic RNAs.   Table 3 shows the top twenty dysregulated lncRNA during phenotypic switch in duced by stent and PDGF. Among the top ten upregulated genes six have not previousl been annotated, and there are two antisense transcripts and two long intervening noncod ing RNAs (lincRNAs) as well. Meanwhile, in the downregulated group there is only on novel LncRNA along with four antisense, two lincRNAs, and two sense-intronic RNAs. Table 3. Top ten upregulated and downregulated lncRNAs in HUASMCs exposed to stent injury.

Gene Enrichment Analysis of DE mRNA and Cis-Target mRNA of lncRNAs
GO analysis of DE mRNA showed significant enrichment for terms related to collagen fibril organization (GO: 0030199), extracellular matrix organization (GO: 0030198), regulation of chemokine-mediated signaling pathway (GO: 0070099), and regulation of smooth muscle cell proliferation (GO: 0048660). Moreover, notable KEGG terms included protein digestion and absorption, IL-17 signaling pathway, AGE-RAGE signaling pathway in diabetic complications, and TGF-beta signaling pathway, among others ( Figure 5A,B).
The complementary pairing of lncRNA and mRNA bases as well as target gene prediction of LncRNA was performed using the LncTar tool. LncTar uses the complementary sequence between LncRNA and mRNA to calculate the free energy of the pairing site. Additionally, Cis-target gene prediction analysis was performed; thus, the protein coding genes of the 100K region upstream and downstream of DE lncRNA were screened out as its target genes. Finally, the of Cis and LncTar target genes were combined for subsequent analysis.
There were 1095 target mRNAs identified, 544 for downregulated lncRNAs and 551 for upregulated lncRNAs. GO analysis of lncRNA-mRNA target genes showed enrichment for terms related to wound healing, spreading of epidermal cells (GO:0035313), negative regulation of pathway-restricted SMAD protein phosphorylation (GO:0060394), collagen biosynthetic process (GO:0032964), and blood vessel endothelial cell proliferation involved in sprouting angiogenesis (GO:0002043) ( Figure 5C). On the other hand, the main enriched KEGG pathways included extracellular matrix-receptor interaction, fatty acid biosynthesis, and protein digestion and absorption ( Figure 5D).

Competing Endogenous RNA (ceRNA) Network Construction
According to the ceRNA theory, a pair of co-expressed lncRNA and mRNA compete for miRNA target sites; therefore, lncRNAs serve as endogenous sponges that target one or more miRNAs, affecting the post-transcriptional regulation of mRNAs. The ceRNA network was constructed with 1241 interactions comprising 9 miRNAs, 39 lncRNAs, and 158 mR-NAs; 4 downregulated miRNAs interact with 15 upregulated lncRNAs and 49 upregulated mRNAs ( Figure 6A). Meanwhile, 5 upregulated miRNAs interact with 24 downregulated lncRNAs and 109 downregulated mRNAs ( Figure 6B).

PPI Network Analysis
We conducted a protein-protein interaction (PPI) network analysis to explore the most significant clusters of DE mRNA among the 150 mRNA in the ceRNA network. The STRING database version 11.5 was used to obtain significant interactions, and the resulting network was visualized in Cytoscape ( Figure 7A). The MCC (maximum clique centrality) method from the cytoHubba app in Cytoscape was used to screen for key proteins ( Figure 7B). The top upregulated genes were HIF1A, MMP1, and FGF2 while the top downregulated genes were TP53, ELN, and BCL2L11. Figure 5. Gene ontology analysis and functional enrichment of KEGG terms for differentially expressed mRNA (panels (A) and (B)) and mRNA cis-target of differentially expressed lncRNAs (panels (C) and (D)); color intensity depicts significance, while ellipse size represents the number of genes.

Competing Endogenous RNA (ceRNA) Network Construction
According to the ceRNA theory, a pair of co-expressed lncRNA and mRNA compete for miRNA target sites; therefore, lncRNAs serve as endogenous sponges that target one or more miRNAs, affecting the post-transcriptional regulation of mRNAs. The ceRNA network was constructed with 1241 interactions comprising 9 miRNAs, 39 lncRNAs, and 158 mRNAs; 4 downregulated miRNAs interact with 15 upregulated lncRNAs and 49 up-   E n rich ed term s Im p u t g en es D Figure 5. Gene ontology analysis and functional enrichment of KEGG terms for differentially expressed mRNA (panels (A,B)) and mRNA cis-target of differentially expressed lncRNAs (panels (C,D)); color intensity depicts significance, while ellipse size represents the number of genes.

PPI Network Analysis
We conducted a protein-protein interaction (PPI) network analysis to explore the most significant clusters of DE mRNA among the 150 mRNA in the ceRNA network. The STRING database version 11.5 was used to obtain significant interactions, and the resulting network was visualized in Cytoscape ( Figure 7A). The MCC (maximum clique centrality) method from the cytoHubba app in Cytoscape was used to screen for key proteins ( Figure 7B). The top upregulated genes were HIF1A, MMP1, and FGF2 while the top downregulated genes were TP53, ELN, and BCL2L11.

PPI Network Analysis
We conducted a protein-protein interaction (PPI) network analysis to explore the most significant clusters of DE mRNA among the 150 mRNA in the ceRNA network. The STRING database version 11.5 was used to obtain significant interactions, and the resulting network was visualized in Cytoscape ( Figure 7A). The MCC (maximum clique centrality) method from the cytoHubba app in Cytoscape was used to screen for key proteins ( Figure 7B). The top upregulated genes were HIF1A, MMP1, and FGF2 while the top downregulated genes were TP53, ELN, and BCL2L11. The network was created by STRING and visualized in Cytoscape. Hub protein network (panel (B)) according to CytoHubba MCC coefficient. Red and green represent upregulated and downregulated proteins, respectively; circle size represents MCC score. PPI, protein-protein interaction; ceRNA, competing endogenous RNA.
according to CytoHubba MCC coefficient. Red and green represent upregulated and downregulated proteins, respectively; circle size represents MCC score. PPI, protein-protein interaction; ceRNA, competing endogenous RNA.

Discussion
Vascular smooth muscle cell phenotypic switch is a dynamic process. Individual cells within a tissue can be found in intermediate steps between two opposite edges [22]. In our experimental model, HUASMC exposed to bare metal stent and PDGF-BB changed their cytoarchitecture and gained proliferative and migratory capacities.  (panel (B)). Pearson correlation coefficient for lnrRNA-mRNA interaction (r > 0.5). Red and green represent upregulated and downregulated RNAs, respectively, grey lines indicate interactions. ceRNA, competitive endogenous RNA; lncRNA, long non-coding RNA; miRNA, microRNA; mRNA, messenger RNA.

Discussion
Vascular smooth muscle cell phenotypic switch is a dynamic process. Individual cells within a tissue can be found in intermediate steps between two opposite edges [22]. In our experimental model, HUASMC exposed to bare metal stent and PDGF-BB changed their cytoarchitecture and gained proliferative and migratory capacities.
Santin et al. exposed VSMCs to stent fragments in vitro, resembling the stent restenosis process [23]. Bare metal stent increased PDGF-BB receptor expression in VSMCs; moreover, under this stimulus, VSMCs secreted growth factors, including PDGF-BB. Cells in direct  For the stent injury model, a piece of bare metal stent was set gently over the HUASMCs. Platelet-derived growth factor (PDGF-BB) 20ng/mL was added to conditioning media.

Cell Migration
Cell migration capability was measured by wound healing assay. When HUASMCs grew to 90% confluency, they were scratched with a sterile 100 μL pipette tip. Cell debris were removed by washing with Hank's balanced salt solution. Cells were observed after every 12 h under an inverted microscope until confluency was reached. The assay was performed in triplicate [75]. Images were processed with the software Image Processing and Analysis in Java (ImageJ v.1.48). The proportion of covered area was determined for each condition using the plug-in MRI_Wound_Healing_Tool [76].

Cell Cycle Analysis
HUASMCs in both experimental conditions were harvested, centrifuged and pelleted. Ethanol 70% at 4 °C was added drop by drop under continuous agitation. Cells were washed twice with PBS and 2% fetal bovine serum and resuspended in 300 μL medium containing 20 μL of RNasa A (Invitrogen, Waltham, MA, USA) and 2 μL propidium iodide (Sigma Aldrich, Burlington, MA, United States ). After mixing, cells were incubated for 30 min in the dark [77]. Cell cycle analysis was performed on a FACs Canto II (Becton-Dickinson, Minato, Tokyo, Japan).

RNA Extraction
Total RNA was extracted from cells using TRIzol reagent (Invitrogen) (Waltham, MA, USA) according to the manufacturer's protocol. RNA grade Glycogen 1 μg/μL (Thermo Fisher Scientific, Inc., Waltham, MA, USA) was added for enhanced RNA precipitation. RNA purity and quantification were assessed in NanoDrop one and RNA integrity (RIN) was assessed using an Agilent 4200. (Figure S1).

RNAseq Analysis
Ribosomal depleted unstranded libraries were prepared by CD Genomics (Shirley, NY, USA) and sequenced with an Illumina HiSeq X10 PE150 (paired end). High-quality clean data were obtained after a series of quality controls that included (1) removing reads Stent-

Cell Migration
Cell migration capability was measured by wound healing assay. When HUASMCs grew to 90% confluency, they were scratched with a sterile 100 µL pipette tip. Cell debris were removed by washing with Hank's balanced salt solution. Cells were observed after every 12 h under an inverted microscope until confluency was reached. The assay was performed in triplicate [75]. Images were processed with the software Image Processing and Analysis in Java (ImageJ v.1.48). The proportion of covered area was determined for each condition using the plug-in MRI_Wound_Healing_Tool [76].

Cell Cycle Analysis
HUASMCs in both experimental conditions were harvested, centrifuged and pelleted. Ethanol 70% at 4 • C was added drop by drop under continuous agitation. Cells were washed twice with PBS and 2% fetal bovine serum and resuspended in 300 µL medium containing 20 µL of RNasa A (Invitrogen, Waltham, MA, USA) and 2 µL propidium iodide (Sigma Aldrich, Burlington, MA, USA). After mixing, cells were incubated for 30 min in the dark [77]. Cell cycle analysis was performed on a FACs Canto II (Becton-Dickinson, Minato, Tokyo, Japan).

RNA Extraction
Total RNA was extracted from cells using TRIzol reagent (Invitrogen) (Waltham, MA, USA) according to the manufacturer's protocol. RNA grade Glycogen 1 µg/µL (Thermo Fisher Scientific, Inc., Waltham, MA, USA) was added for enhanced RNA precipitation. RNA purity and quantification were assessed in NanoDrop one and RNA integrity (RIN) was assessed using an Agilent 4200. (Figure S1).

RNAseq Analysis
Ribosomal depleted unstranded libraries were prepared by CD Genomics (Shirley, NY, USA) and sequenced with an Illumina HiSeq X10 PE150 (paired end). High-quality clean data were obtained after a series of quality controls that included (1) removing reads containing adapters; (2) removing reads containing n > 10%, where n represents the base which cannot be determined; (3) removing reads containing low quality-base representing over 50% of the total base. Between 55 million and 73 million paired reads were obtained for the samples, allowing a deep coverage. The filtered clean reads were mapped to the reference genome by HISAT2 [78]. (Table S1).

Gene Expression Quantitative Analysis
Expression levels were measured as Fragments per Kilobase of transcript per Million mapped reads (FPKM). The quantification and gene expression levels of the transcripts were assessed using the Cuffquant and Cuffnorm components of Cufflinks software, using the mapped reads' positional information on the gene [79].
Differential expression was performed utilizing DESeq2 [80]. We considered a threshold of absolute Fold Change (FC) ≥ 2 and a False Discovery Rate (FDR) < 0.05 to identify significant changes between two conditions. Hierarchical clustering based on the log2(FC) of differentially expressed genes was evaluated with DESeq2 Principal component analysis tool. (Dataset S1). Heatmaps were produced in R using the bioconductor package ComplexHeatmap [81].

Filtering of Candidate lncRNAs
To screen for lncRNA within the merged transcripts, we performed basic filtering including five steps: (1) Screening exons, filtering low-expression and low-quality single exon transcripts; (2) Selecting transcripts which were longer than 200 bp and had more than two exons; (3) Screening transcripts with known annotations using Cuffcompare; (4) Transcript expression level filtering by calculating the expression level of each transcript and choosing those with FPKM ≥ 0.5; (5) Coding potential filtering; by combining several mainstream coding potential analysis methods (CNCI [82], CPC [83], and Pfam [84]), it was determined which transcripts had no coding potential ( Figure S2).

Prediction of lncRNA Target Genes
The basic principle of Cis target gene prediction is that the function of lncRNA is related to its adjacent protein-coding gene; thus, the protein coding genes of the 100 Kb region upstream and downstream of a lncRNA were screened out as its target gene. Other criteria were based on the complementary pairing of lncRNA and mRNA bases. LncTar uses the complementary sequence between LncRNA and mRNA to calculate the free energy and normalized free energy of the pairing site. Thus, the target gene of the LncRNA is that wich normalized free energy is below the threshold. Finally, the target genes of Cis and LncTar were combined for subsequent analysis [85].

Functional Analysis
Gene ontology (GO; Gene Ontology Consortium, 2000) enrichment analysis of differentially expressed genes was carried out using topGO [86]; Fisher's exact test was used to calculate p-values. The specific genes that are involved in the major metabolic pathways and signal transduction pathways can be determined by Pathway significant enrichment analysis. The Kyoto Encyclopedia of Genes and Genomes (KEGG) is a database resource for understanding the high-level functions and utilities of biological systems. Therefore, we employed the KOBAS software to detect the enrichment of differentially expressed mRNA and lncRNA target genes. The pathways with FDR ≤ 0.05 were defined as the significantly enriched pathways [87].

Construction of lncRNA-miRNa-mRNA Competing Endogenous RNA (ceRNA) Network
We constructed a ceRNA network based on the hypothesis that lncRNAs regulate the activity of mRNAs by sequestering and binding miRNAs, thereby acting as miRNA sponges [88]. All transcripts were filtered according to expression profiles of fold change log2(FC) > 2.5 for lncRNAs and log2(FC) > 1.5 for miRNA and mRNA, with FDR < 0.01.
More stringent criteria were used in order to ascertain the biological significance of DE transcripts, as lncRNA's low abundance and nuclear localization may hamper its ceRNA ability [89]. To construct the ceRNA network, we predicted miRNA-mRNA and lncRNA-miRNA interactions based on miRTarBase [90], miRDB [91], and TargetScan [92]. The resulting matrix was intersected with our sequencing data, and miRNAs with opposite expression levels to the target lncRNAs and mRNAs were selected. Finally, the lncRNA-miRNA-mRNA network was reconstructed and visualized using Cytoscape software v3.7.2 (San Diego, CA, USA).

Protein-Protein Interaction Network
Target mRNAs included in the ceRNA network were analyzed for protein-protein interaction (PPI) with the STRING online search tool, setting a combined score of ≥0.4. Cytoscape v3.7.2 was used to visualize the PPI network, and significant enriched interactions were selected by CytoHubba using maximum clique centrality score (MCC).