Increasing Upstream Chromatin Long–Range Interactions May Favor Induction of Circular RNAs in LysoPC-Activated Human Aortic Endothelial Cells

Circular RNAs (circRNAs) are non-coding RNAs that form covalently closed continuous loops, and act as gene regulators in physiological and disease conditions. To test our hypothesis that proatherogenic lipid lysophosphatidylcholine (LPC) induce a set of circRNAs in human aortic endothelial cell (HAEC) activation, we performed circRNA analysis by searching our RNA-Seq data from LPC-activated HAECs, and found: (1) LPC induces significant modulation of 77 newly characterized cirRNAs, among which 47 circRNAs (61%) are upregulated; (2) 34 (72%) out of 47 upregulated circRNAs are upregulated when the corresponding mRNAs are downregulated, suggesting that the majority of circRNAs are upregulated presumably via LPC-induced “abnormal splicing” when the canonical splicing for generation of corresponding mRNAs is suppressed; (3) Upregulation of 47 circRNAs is temporally associated with mRNAs-mediated LPC-upregulated cholesterol synthesis-SREBP2 pathway and LPC-downregulated TGF-β pathway; (4) Increase in upstream chromatin long-range interaction sites to circRNA related genes is associated with preferred circRNA generation over canonical splicing for mRNAs, suggesting that shifting chromatin long-range interaction sites from downstream to upstream may promote induction of a list of circRNAs in lysoPC-activated HAECs; (5) Six significantly changed circRNAs may have sponge functions for miRNAs; and (6) 74% significantly changed circRNAs contain open reading frames, suggesting that putative short proteins may interfere with the protein interaction-based signaling. Our findings have demonstrated for the first time that a new set of LPC-induced circRNAs may contribute to homeostasis in LPC-induced HAEC activation. These novel insights may lead to identifications of new therapeutic targets for treating metabolic cardiovascular diseases, inflammations, and cancers.


INTRODUCTION
Cardiovascular diseases remain some of the most prevalent health challenges today. Predominant among these, atherosclerosis and resultant ischemic heart disease, stroke, and peripheral artery disease are estimated to cause 8.9 million deaths per year, making it the single most significant cause of death worldwide (Wang H. et al., 2016). We and others have reported that hyperlipidemia, together with other CVD risk factors, such as hyperglycemia, chronic kidney disease, obesity, and hyperhomocysteinemia (HHcy), promotes atherosclerosis development via several mechanisms. These mechanisms include endothelial cell (EC) activation and injury (Shao et al., 2014;Yin et al., 2015;Li et al., 2016aLi et al., , 2018a; monocyte recruitment and differentiation (Combadière et al., 2008;Fang et al., 2014); decreased function of regulatory T cells (Tregs) (Ait-Oufella et al., 2006;Yan et al., 2008;Li et al., 2017); transdifferentiation of Tregs into antigen-presenting cell (APC) like Tregs ; impaired vascular repairability of bone marrow-derived progenitor cells (Du et al., 2012;Barbier et al., 2015;Li et al., 2016b); increased migration and proliferation of vascular smooth muscle cells (Monroy et al., 2015;Ferrer et al., 2016), and high fatinduced adipocyte hypertrophy and metabolic healthy obesity (Virtue et al., 2017). Human aortic endothelial cell (HAEC) activation is widely understood to contribute substantively to the early development of atherosclerosis by increasing leukocyte recruitment through upregulation of adhesion molecules such as ICAM1 (Tabas et al., 2015;Vestweber, 2015;Gimbrone and García-Cardeña, 2016). We have recently reported that exposure to proatherogenic stimuli, such as the proinflammatory lysophospholipid lysophosphatidylcholine (LPC), can induce prolonged HAEC activation as part of early atherogenesis (Li et al., 2016a(Li et al., ,c, 2018b; this phenomenon is one of a larger class of conditional danger-associated molecular patterns (DAMPs) that trigger the onset of oxidative, autoimmune, and inflammatory tissue damages (Wang X. et al., 2016;Shao et al., 2017). However, the regulatory mechanisms underpinning these physiological changes remain poorly understood; and whether and how noncoding RNAs (ncRNAs) mechanisms control HAEC activation remains poorly characterized.
For many years, proteins have represented the primary functional end product of genetic information, though the genes that encode them account for less than 2% of the genome (Anastasiadou et al., 2018). As new sequencing technology advances, as many as six types of ncRNAs have been identified, including three small ncRNA types [i.e., RNA interference (RNAi), tRNA-derived stress-induced RNA (tiRNA), and small nuclear RNA (snRNA)] and three long RNAs types [i.e., transfer RNA (tRNA), ribosomal RNA (rRNA), and long non-coding RNAs (lncRNA)] (Gomes et al., 2018). RNAi can be categorized into five subgroups including small interfering RNA (siRNA), piwi-interacting RNA (piRNA), and microRNA (miRNA), repeat associated siRNA (rasiRNA), and agotron. LncRNAs can be further classified into six subtypes: sense lncRNA, intronic lncRNA, antisense lncRNA, intergenic lncRNA, enhancer RNA (eRNA), and circular RNA (circRNA). These ncRNAs together constitute almost 60% of the transcriptional output in human cells, and they have been shown to regulate cellular processes and pathways in developmental and pathological contexts through mRNA degradation or modification, transcription regulation, translation regulation, and miRNA sponging, among other functions (Mattick and Makunin, 2006;Bartel, 2009;Cech and Steitz, 2014;Lasda and Parker, 2014). We previously reported that a group of anti-atherogenic miRNAs might suppress proatherogenic gene expression (Virtue et al., 2011Mai et al., 2012). We also found that microRNA-155 (miR155) is among several miRs significantly upregulated in atherogenic apolipoprotein E deficient (ApoE−/−) mouse aortas. MiR155 promotes HAEC activation and mouse atherosclerosis, whereas the deficiency of miR155 in ApoE−/− background leads to high fat-induced adipocyte hypertrophy and establishes a novel metabolically healthy obese (MHO) mouse model (Virtue et al., 2017). Among 109 miRNAs that we examined in all four hyperlipidemia-related diseases (HRDs), namely atherosclerosis, non-alcoholic fatty liver disease (NAFLD), obesity and type II diabetes (T2DM), miR-155 and miR-221 are significantly modulated in these HRDs. MiR-155 is significantly upregulated in atherosclerosis and decreased in other HRDs. MiR-221 is increased in three HRDs but reduced in obesity. These findings led to our new classification of types I and II MHOs, which are regulated by miR-221 and miR-155, respectively. We also found that the proinflammatory adipokine, resistin, is significantly increased in white adipose tissues (WAT) of the MHO mice, revealing our newly proposed, miR-155-suppressed "secondary wave inflammatory state (SWIS), " characteristic of MHO transition to classical obesity (CO) (Johnson et al., 2018). Our and others' results demonstrate that ncRNAs play significant roles in regulating the pathogenesis and progression of cardiovascular diseases and other diseases.
Two types of lncRNAs (i.e., linear and circular forms) have been found in the lncRNA family. In contrast to previous models, recent work has shown that circRNAs are relatively abundant and may be the predominant transcripts of hundreds of genes. The spliceosome machinery produces circRNAs by direct FIGURE 1 | (A) Lysophosphatidylcholine (LPC)-induced top 5% significantly changed circular RNAs (circRNAs) are studied in human aortic endothelial cells (HAECs) in comparison to that of control HAECS. This is a flow chart to outline how we analyzed the 5% top significantly changed circular RNAs from LPC-treated HAECs in comparison to that from control HAECs. (B) The majority (48 out of 51,94.1%) of circRNA expression increases while related mRNA expressions remain constant or decrease. (C) There is no significant differences in exon numbers among circRNA groups. (D) There are no overlaps among the top 10 pathways in the following four groups.
back-splicing of exons or introns within pre-mRNA or by exon skipping. As a result of a direct back-splicing event, an exon is not associated with an adjacent downstream exon, as seen in linear splicing, but with an upstream exon or intron. This back-splicing gives rise to circular RNA molecules harboring exons or introns that are out-of-order from the genomic context ( Figure 1A). These circular molecules are ubiquitously expressed in human and more than 30,000 different circRNAs have been identified. They are regulated in several disease states and contribute to disease development and progression (Gomes et al., 2018). Although the importance of circRNAs in regulating EC function has been reported, the comprehensive analyses of circRNAs in EC activation, EC function, and cardiovascular diseases remain at their infancy; additionally, the characterization of circRNAs in atherosclerosis-relevant HAECs has not been reported (Maass et al., 2017;Huang et al., 2018;Shang et al., 2018).
Regardless of the progress, we identified several important knowledge gaps to be filled in this study. We hypothesized that circRNAs might play an important role in mediating homeostatic status in HAECs activated by proatherogenic LPC. We analyzed our RNA-Seq raw data sets, obtained from control-, and LPC-treated HAECs (Li et al., 2018b), using two wellaccepted circular RNA databases such as CircExplorer2 and CIRI to obtain a list of significantly changed circular RNAs (circRNAs) (Gao et al., 2015;Zhang et al., 2016) 4. We found that this list of 77 circRNAs may play homeostatic roles in activated HAECs via novel mechanisms including exon usage competition, miRNA sponging (Hansen et al., 2013;Memczak et al., 2013) and potential generation of short proteins for protein interaction disruption. We also found that novel mechanisms, such as (a) flanking intron homology, (b) LPC-induced modulation of spliceosome components, and (c) upstream chromatin long-range interaction, play important roles in promoting the generation of circRNAs in LPC-activated HAECs. These novel insights may lead to the identification of new therapeutic targets for treating metabolic cardiovascular diseases, inflammations, and cancers.

Generation of CircRNA Datasets
Raw RNA-Seq data was obtained from a previous study (Li et al., 2018b). Raw reads generated by the Illumina HiSeq control software were assessed using FastQC (Andrews, 2010). Samples with poor sequence quality were discarded. Sequence reads were mapped to the hg38 genome using STAR (Dobin et al., 2013). Fusion junction reads were parsed using CIRCexplorer2 and captured back-splicing junction reads were then annotated with UCSC annotation files (Zhang et al., 2016). For additional detections, sequence reads were mapped to the hg38 genome using the Burrows-Wheeler Alignment tool (Li and Durbin, 2009), circRNAs were identified using CIRI (Gao et al., 2015). Gene annotation was performed using org.Hs.eg.db (Carlson, 2017). A list of circRNAs common to both control and LPC-stimulated HAECs was generated using the VLOOKUP function in Microsoft Excel. To filter out significantly changed circRNAs, ratios of control to LPC-stimulated circRNA read numbers were calculated, and base 2 logarithms taken of these ratios. Mean and standard deviation of logarithm values were calculated using MATLAB and used to generate a 95% confidence interval (mean ± 2 * SD). Using previously obtained mRNA expression data (Li et al., 2018b), an mRNA fold change confidence interval was similarly calculated using the ten housekeeping genes: C1orf43, CHMP2A, GAPDH, EMC7, GPI, PSMB2, PSMB4, RAB7A, SNRPD3, and VPS29; and it was used to separate circRNAs by corresponding mRNA expression change.

Comparison of Flanking Intron Sequences
The generated circRNA data included flanking intron coordinates; a Python script was used to parse these coordinates and feed them into an NIH NCBI BLAST+ command line query (blastn -max_hsps 1 -outfmt "10 pident e-Value bitscore") comparing the sequences of 3 and 5 flanking intron sequences using local copies of hg38 human genome primary assembly chromosome sequences 1 (Camacho et al., 2009). This returned percent identity, the expected value (number of matches of equal strength expected by chance), and bitscore values (normalized score for match strength) for each pair of flanking introns. These were then exported to a text file and imported into a Microsoft Excel spreadsheet for tabulation. Bitscore values were further used for statistical analysis (described below under "Statistical Analysis").

Matching CircRNAs to Database Entries
The genomic length of each significantly changed circRNA was calculated by subtracting its genomic start coordinate from its genomic end coordinate. A MATLAB script was then used to find circRNAs with matching length and gene locus from a local copy of the circRNAdb database and to export all results to a Microsoft Excel spreadsheet (Schliebs et al., 1996;Chen et al., 2016). Notably, genomic coordinates in circRNAdb were given with respect to the older hg37 human genomic annotation. A local copy of the hg37 genomic annotation was therefore obtained for comparison 2 . Then, NIH NCBI BLAST+ queries between original hg38 coordinates and CircRNAdb hg37 coordinates (blastn -max_hsps 1 -outfmt "10 pident") were performed using a MATLAB script (Schliebs et al., 1996;Camacho et al., 2009), taking only the single best match for each pair of sequences, the percent identity values of which confirmed sequence alignment for all length matches found. CircInteractome also uses the hg37 human genome annotation but a different set of IDs; matches in CircInteractome were determined by hand due to incomplete downloadable data sets.

Analysis of Long-Range Interactions
A complete list of long-range chromatin interaction sites in the human genome was obtained from the 4DGenome database 3 as a tabulated text file (Teng et al., 2015). The grep command line utility was used to filter for entries detected using Hi-C methodology and involving at least one circRNA-related gene. The resulting filtered data was imported into Microsoft Excel and raw interaction distances calculated as the differences between gene start coordinates. An AWK script was used to determine whether the circRNA-related gene was downstream or upstream of its partner in each interaction pair and to add this information to the data file. The signs of distance values were then updated, with downstream entries designated as positive and upstream values designated as negative, using a Python script. These updated distance values were separated into the same six groups as their corresponding circRNArelated genes and used in pairwise two-sample Kolmogorov-Smirnov tests (described further below under "Statistical Analysis"). Distance distributions for all upregulated and all downregulated circRNAs were compared by groups overall as well as by only upregulated and downregulated circRNAs with increased, unchanged, or decreased corresponding mRNA expression, respectively.

Characterizing Open Reading Frame Sequences
Open reading frame (ORF) peptide sequences and internal ribosomal entry site (IRES) data were obtained to identify significantly changed circRNAs from circRNAdb (Chen et al., 2016). A MATLAB script was used to call remote NIH NCBI BLAST+ peptide alignments for all ORF sequences against the NCBI non-redundant protein sequences database (nr 4 ) (Schliebs et al., 1996;Camacho et al., 2009), restricted to entries for Homo sapiens (blastp -db nr -remote -entrez_query "Homo sapiens [Organism]" -max_target_seqs 1 -outfmt "10 sacc pident e-Value"). This returned the subject accession, percent identity, and expected value of the single best database match, which was then exported with the gene name and circRNAdb ID to a Microsoft Excel spreadsheet. The obtained accessions were individually and manually searched on the NCBI protein database 5 to determine if they corresponded to canonical mRNA transcripts from the same gene locus as the corresponding circRNA. A modified version of the hg38/hg37 comparison MATLAB script was also used to run NIH NCBI BLAST+ searches for the Kozak consensus sequence gccRccAUGG in genomic sequences of all significantly changed circRNAs (Kozak, 1987;Schliebs et al., 1996;Camacho et al., 2009).

Pathway Analysis
QIAGEN Ingenuity Pathway Analysis (IPA) software, which constructs predicted upstream and downstream causal networks for input datasets from a curated research literature base, was used to elucidate potential downstream pathways for sponged miRNAs (Krämer et al., 2014). Lists of miRNAs were input into IPA and ran through the miRNA target filter to generate a list of potential mRNA targets. The list of mRNA target genes was then run through an IPA core analysis. All canonical downstream pathways returned in the resulting output were exported into a Microsoft Excel spreadsheet. The top ten pathways were extracted for further qualitative consideration.

Graphical Figure Generation
For three-group Venn diagrams, lists of genes were input into an online Venn diagram generator ( 6 Evolutionary Genomics, Ghent University, Gent, Belgium) (Draw Venn Diagram). This tool was used to produce both diagrams and lists of overlapped genes between groups. For six-group Venn diagrams, gene lists were input into the InteractiVenn online tool 7 to generate diagrams (Heberle et al., 2015), while the above Ghent University tool was again used to produce lists of overlapped genes between groups (Draw Venn Diagram). Explanatory and conceptual graphics were produced using Microsoft Paint.

Statistical Analysis
Descriptive summary statistics were reported by group. Data were checked for normality assumption and, if found to be not normally distributed, subsequently transformed using various functions such as the log10 and cubic-root to find the optimal transformation for the underlying chromatin long-range interaction distance distribution data. Chromatin long-range interaction distance density functions were then estimated and plotted by group under the optimal transformation using the non-parametric kernel density approach with a normal weight function (Jones et al., 1996;Hollander and Wolfe, 1999;Silverman, 2018). Pairwise comparisons of median chromatin long-range interaction distance among the six groups were performed with multiple comparison adjustments using the Dwass, Steel, and Critchlow-Fligner method based on the Wilcoxon test for downstream and upstream data separately (Douglas and Michael, 1991;Conover, 1999;Hollander and Wolfe, 1999;Shalabh, 2011). Pairwise comparisons between groups for chromatin long-range interaction distance distribution, median distance location shift, and distance distribution scale were implemented using the Kolmogorov-Smirnov two-sample test, Hodges-Lehmann estimation method and Fligner-Policello test, and Ansari-Bradley test, respectively, again for downstream and upstream data separately (Lehmann, 1963;Douglas and Michael, 1991;Conover, 1999;Hollander and Wolfe, 1999;Shalabh, 2011). SAS version 9.4 was used to perform these analyses and generate density function plots for the chromatin long-range interaction distance data.
For rest of the data, single-factor ANOVA and nonparametric Kruskal-Wallis tests were conducted for all multigroup analyses using the Real Statistics Resource Pack add-in for Microsoft Excel 8 (Zaiontz, 2013). For confidence intervals, MATLAB was used to calculate mean and standard deviations for data sets (Schliebs et al., 1996). Probabilities of ratios of upregulation to downregulation were calculated by summing binomial coefficients and dividing by the appropriate power of 2 in MATLAB (Schliebs et al., 1996), as given by the , where u is the total number of genes considered and n is the large number of upregulated or downregulated genes. Our preliminary data processing detected 3,170 distinct circRNAs in control cells and 3,486 circRNAs in a proatherogenic lipid LPC-treated HAECs as we reported (Li et al., 2016a). To better analyze our data set, we derived a list of 1,093 circRNAs found in both the control and LPC-treated data sets before calculating read number ratio -the quotient of the read number in the LPCtreated data set and the corresponding read number in the control data set -for each circRNA in the list to serve as a measure of expression change ( Figure 1A). We then took the base-2 logarithm of each readnumber ratio (log ratio) and calculated a 95 percent confidence interval (−0.3815 ± 2.3395) for log ratio, taking all the 77 circRNAs with log ratios outside this confidence interval to have significantly changed levels of expression in LPC-treated HAECs (Table 1A). Then, we separated these circRNAs into two groups of increased and decreased circRNA expression. To assess potential relationships between circRNAs and their corresponding mRNAs published, we then further divided each of these groups into three subgroups based on increased, unchanged, or decreased expression of corresponding mRNA transcripts, with all fold change levels within a 95 percent confidence interval (1.01 ± 0.1) determined by the housekeeping genes C1orf43, CHMP2A, GAPDH, EMC7, GPI, PSMB2, PSMB4, RAB7A, SNRPD3, and VPS29 considered to indicate unchanged mRNA expression (Table 1B), thereby obtaining six groups of significantly changed circRNAs classified by changes in both circRNA and corresponding mRNA expression.

LPC
We also matched each circRNA, if possible, with its corresponding ID in the circRNA database CircRNAdb (Chen et al., 2016), to provide for more consistent identification, and for later analysis of the database's open reading frame data. Of note, among 51 circRNAs increased, 48 out of 51 related mRNAs (94.1%) are either unchanged (38 out of 51, 74.5%) or decreased (10 out of 51, 19.5%), while only three increased both mRNA and circRNA expression (5.9%, Figure 1B). Similarly, among 26 circRNAs decreased, 19 out of 26 related mRNAs (73.1%) are either unchanged (15 out of 26, 57.7%) or increased (4 out of 26, 15.4%), while only seven decreased both mRNA and circRNA expression (26.9%, Figure 1B). The results suggest that a dominant inversely correlative relationship between circRNA expression and mRNA expression levels potentially indicates exon usage competition during splicing for either circRNA expression or mRNA expression.  3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22  We also compared the average number of exons per circRNA transcript of each group using single-factor ANOVA to assess any potential impact of circRNA length on expression patterns, but no significant differences were found (p = 0.788) ( Figure 1C). Considering the slightly right-skewed distributions of exon numbers, we also performed a Kruskal-Wallis test of the same data, which also failed to indicate significance (p = 0.46309). These results suggest that transcript length has little to no role in modulating circRNA expression following LPC treatment in HAECs. Taken together, our results suggest that the proatherogenic lipid LPC induces a list of circRNAs in HAECs; and that a dominant inversely correlative relationship between circRNA and mRNA expression levels potentially indicates exon usage competition during splicing.

CircRNA-Related Genes Do Not Share Any Top Pathways With LPC-Induced mRNA Pathways; and Genes Associated With Significantly Changed CircRNAs Exhibit Little Relationship With Potential Effects of Inflammation-Related Cytokine Stimulation
To determine the signaling pathways associated with circRNArelated genes, we performed the IPA. As shown in Table 1C, we found the following: first, six top pathways related to the unchanged genes which related circRNAs were increased, including pyridoxal 5 -phosphate salvage pathway, Fcγ receptor-mediated phagocytosis in macrophage and monocyte, transforming growth factor-β (TGF-β) signaling, salvage pathways of pyrimidine ribonucleotides, epithelial adherens junction signaling, and Wnt/β-catenin signaling; second, ten top pathways related to the unchanged mRNAs which related circRNAs were decreased, including ephrin B signaling, focal adhesion kinase (FAK) signaling, axonal guidance signaling, ephrin receptor signaling, integrin signaling, interleukin-15 (IL-15) production, tRNA charging, semaphoring signaling in neurons, ephrin A signaling, and regulation of cellular mechanics by calpain protease; third, ten top pathways related to the unchanged RNAs which related circRNAs were increased, including pyridoxal 5 -phosphate salvage pathway, Fcγ receptor-mediated phagocytosis, TGF-β signaling, salvage pathways of pyrimidine ribonucleotides, epithelial adherens junction signaling, Wnt/β-catenin signaling, actin cytoskeleton signaling, cell cycle control of chromosomal replication, unfolded protein response, remodeling of epithelial adherens junctions; and fourth, ten top pathways related to the decreased mRNAs which related circRNAs were decreased, including role of osteoblasts, osteoclasts, and chondrocytes in rheumatoid arthritis, neuroinflammation signaling pathway, cardiomyocyte differentiation via bone morphogenetic protein (BMP) receptors, April mediated signaling, B cell activation factor signaling, Wnt/Ca++ pathway, netrin signaling, BMP signaling pathway, regulation of IL-2 expression in activated and anergic T lymphocytes and factors promotion of cardiogenesis in vertebrates. The Vann Diagram analysis found that no signaling pathways are overlapped among the four groups ( Figure 1D).
The results suggest that, although we cannot determine directly circRNAs functional pathways in this LPC-induced HAEC CircRNAs may play roles in weakening the roles of related mRNAs via their roles in exon usage competition. #Of note, there were no top ten pathways identified by the Ingenuity Pathway Analysis in three other groups in Table 1A.
activation, LPC-induced mRNAs-related HAEC activation pathways, that we reported (Li et al., 2018b), have no overlaps with any signaling pathways for circRNAs-related mRNAs. We further hypothesized that LPC-induced circRNA expression could be positively associated with changes in gene expression under other pro-inflammatory stimuli. Thus, to compare the effects of LPC stimulation with that of prototypical inflammation-related responses, we obtained mRNA expression change data for cells stimulated with the pro-inflammatory cytokines tumor necrosis factor alpha (TNFα), interferon gamma (IFN-γ), and interleukin 6 (IL-6) and the anti-inflammatory cytokines TGF-β, IL-4, and IL-10 from the database ArrayExpress (Braun et al., 2013;Teles et al., 2013;Bald et al., 2014;Jin et al., 2014;Kharma et al., 2014;Montoya et al., 2014;Polak et al., 2014Polak et al., , 2017Duncan et al., 2015;Kolesnikov et al., 2015;Tan et al., 2015). We found all circRNA-related genes significantly changed after stimulation in each data set (not shown) and results indicate that cytokines with similar effects did not tend to regulate circRNA-related genes in similar directions, nor did pro-inflammatory and anti-inflammatory cytokines tend to regulate circRNA-related genes in opposing directions (not shown). This in turn suggests a limited or absent canonical role for circRNA-related genes in inflammatory response and fails to indicate general trends for pro-or anti-inflammatory cytokine modulation of circRNArelated genes. However, it neither supports nor disproves the roles of inflammatory signaling in regulating circRNA levels, as circRNA levels are largely independent of the mRNA expression levels used in this analysis and, furthermore, circRNAs may fulfill different physiological functions such as homeostasis than their canonically spliced mRNA counterparts do such as promotion of HAEC activation induced by LPC.

68.8% of the CircRNAs Feature Significant Flanking Intron Homology
To identify the potential mechanism underlying the generation of circRNAs, previous reports found that many circRNA transcripts feature reverse complementary flanking intron sequences, which enable self-pairing that brings 3 and 5 exonic sequences into close proximity, an arrangement that favors head-to-tail alternative splicing and circRNA biogenesis (Figure 2A; Albury et al., 2015). To confirm this pattern in our dataset, we ran NIH NCBI BLAST sequence alignments using NIH NCBI BLAST+ command line tools (blastn -max_hsps 1 -outfmt "10 pident e-Value bitscore") comparing the 3 and 5 flanking intron sequences for each circRNA, as previously described (Albury et al., 2015). We found that 53 out of 77 circRNAs (68.8%) possessed flanking intron sequences with significant reverse complements as determined by expected value (number of matches expected given random sequences) <1 * 10 −20 ( Table 2). This number is lower than that of human 88% circRNA with the Alu repeats in flanking regions (Albury et al., 2015), and that in C. elegans circRNAs (96.3%); and high than that of Zebrafish circRNAs (2.5%) (Ivanov et al., 2015). To test potential association of reverse complementarity of flanking introns with different expression patterns of circRNA and mRNA, we performed a single-factor ANOVA using the bitscores (normalized values reporting strength of alignment between sequences) of circRNAs from each of the six groups sorted previously, assigning a bitscore of zero (indicating no matches better than those resulting from chance) for all circRNAs whose flanking introns did not demonstrate reverse complementarity. This initial analysis indicated significant differences between groups (p = 0.0011), but further inspection revealed the presence of an outstanding data point (bitscore = 4617) in the "circRNA increased, mRNA increased" group, which combined with the small size (n = 3) of the group greatly weakened the generalizability of the test (Figure 2B). We thus performed a non-parametric Kruskal-Wallis test on the same data set, but this failed to reach significance (p = 0.74053). We also performed another single-factor ANOVA with the outstanding point removed, the results of which marginally failed to indicate any significant difference in bitscore between groups (p = 0.07877). Thus, we cannot conclude that there is an association between flanking intron homology and circRNA expression patterns, although that our statistical analysis approaches significance suggests that the area may warrant further study, especially with a larger sample. Our results indeed demonstrate that  We and others reported that downregulation of certain spliceosome components promotes circRNA formation over canonical mRNA splicing (Shen et al., 2017;Wilusz, 2018); and that proinflammatory cytokine TNF-a, inflammation (Liang et al., 2017), autoimmune disease (Xiong et al., 2006b), and tumors modulate alternative splicing (Ng et al., 2004) and spliceosome . To determine whether this potential mechanism plays any roles in shaping LPCmodulated circRNA expression patterns, we compared the list of significantly changed mRNAs from RNA-Seq (determined via housekeeping gene confidence interval, as previously described) with lists of known spliceosome proteins (Deckert et al., 2006;Yang et al., 2006;Behzadnia et al., 2007;Bessonov et al., 2008; Table 3). Six spliceosome components were significantly downregulated after LPC stimulation, while another 11 were significantly upregulated. Given that 51 (66.2%) circRNAs were significantly upregulated while only 26 (33.8%) were significantly downregulated (p = 0.0029, where p represents the probability of a ratio equal to or more extreme than the experimental value, assuming circRNA upregulation and downregulation equally likely) and that among upregulated circRNAs, 10 corresponding mRNAs experienced significant downregulation while only 3 experienced significant upregulation (p = 0.0461, where p represents the probability of a ratio equal to or more extreme than the experimental value, assuming mRNA upregulation and downregulation equally likely) ( Figure 1A and Table 3), this disruption may have been sufficient to shift the splicing patterns from canonical mRNA splicing toward back splicing for circRNA formation (Figure 3), as reported previously (Wilusz, 2018). While the data currently available is insufficient to establish a causal link between LPC-induced spliceosome disruption and a shift from canonical splicing to circRNA formation, it does demonstrate a strong association suitable for more intensive study. Taken together, our results suggest that LPC stimulation disrupts the There are 6 protein downregulation and 11 proteins upregulation.   (Wahl et al., 2009), chromosome conformation capture-on-chip (4C-Seq) (Dekker et al., 2002;Splinter et al., 2011), and chromosome conformation capture carbon copy sequencing (5C-Seq) (Stadhouders et al., 2013) allows determination of the relative frequency of direct physical contact between a pair of linearly separated chromatin segments with 3C-Seq, and genome-wide interactions involving a single anchor region, and interactions involving multiple genomic regions with 4C-Seq and 5C-Seq, respectively ( Figure 4A). Differences in chromatin long-range interaction patterns between genes have previously been hypothesized to influence alternative splicing (Dostie et al., 2006) and the transcription of inflammatory genes such as cytokine (Acuña and Kornblihtt, 2014) and cytokine receptor (Park et al., 2016) and cardiovascular disease causative genes (Deligianni and Spilianakis, 2012). We hypothesize that chromatin long-range interactions differentially regulate the gene promoters to differentiate the mRNA increase from circRNA increase versus decrease, derived from the same gene, under LPC stimulation in HAECs. To test this hypothesis with respect to circRNA biogenesis, we obtained long-range interaction data for all genes related to significantly changed circRNAs from the well-accepted 4DGenome database with a huge collection of 4,433,071 experimentally derived chromatin long-range interactions (Teng et al., 2015). We then calculated the distances between interacting genes with respect to the circRNA-related gene, with distances designated as positive if the circRNA-related gene was localized downstream of its partner and as negative if the circRNA-related gene was localized upstream of its partner. The two-sample Kolmogorov-Smirnov test of the chromatin long-range interaction distance between genes corresponding to downregulated and upregulated circRNAs indicated a significant difference between the two distance distributions (p < 0.001). First, the majority of longrange interaction sites with circRNA-related genes were focused on 10 6 kb upstream and downstream of the circRNA-related gene, suggesting a common feature of the long-range chromatin interaction sites with LPC-induced circRNA-related genes; Second, in the three groups with circRNA related genes such as (a) upregulation; (b) downregulation; or (c) no changes induced by LPC, the group with LPC-induced upregulation of mRNAs had the lowest peak, the group with LPC-induced downregulation of mRNAs had the middle peak, the group with LPC-induced no changed expression of mRNAs had the highest peak, suggesting that high concentration of chromatin long-range interaction sites prevents HAECs from responding to LPC stimulation in changing gene expressions (Figures 4B,C); Third, in both groups of mRNA upregulation and mRNA downregulation induced by LPC, circRNA increase was associated with increased upstream chromatin long-range interaction sites, suggesting that an increase in upstream chromatin long-range interaction sites is associated with circRNA generation presumably by favoring back-splicing in generating circRNAs over canonical splicing in generating mRNAs ( Figure 4C). In contrast, a decrease in upstream chromatin long-range interaction sites is associated with mRNA generation presumably by favoring canonical splicing in generating mRNAs over back-splicing in generating circRNAs ( Figure 4C). Altogether, these results suggest that different long-range interaction mechanisms are likely to govern circRNA upregulation and downregulation. Indeed, upregulated circRNAs tended to interact at somewhat longer distances than downregulated circRNAs did (Figures 4B,C). Since the 4DGenome database contains the experimental data derived from human non-aortic ECs (Teng et al., 2015), the future work will be needed to use circular chromosome conformation capture sequencing (4C-Seq) to examine HAECs stimulated with LPC to map the specific upstream interaction sites for a favorable selection of circRNA over mRNAs. Therefore, our results may suggest that certain ideal long-range interaction distances exist that are more likely to promote circularization; and that favorable distances tend to be greater in magnitude than unfavorable ones.

Six Significantly Changed CircRNAs Sponge miRNAs, Which May Be Associated With Inhibition of Type 2 Diabetes Mellitus, Inhibition of Migration, and Inhibition of Vascular Dysfunction and Immune Responses in Endothelial Cells
A number of circRNAs can sponge miRNAs, preventing miRNAmediated degradation of mRNA transcripts and disruption of canonical signaling and effector pathways (Hansen et al., 2013;Memczak et al., 2013). Thus, by regulating miRNA activity, significantly changed circRNAs can act through canonical pathways, potentially playing a significant role in HAEC activation. We examined the sponging potential of all significantly changed circRNAs using the CircInteractome database (Montefiori et al., 2018), recording two miRNAs with four or more predicted binding sites in a single circRNA transcript, a threshold above which meaningful sponging activity is likely to occur Memczak et al. (2013). Another four significantly changed circRNAs are experimentally shown to sponge miRNAs (Dudekula et al., 2016;Chen et al., 2017;Yan et al., 2017;Wang et al., 2018), for six total circRNAs with miRNA sponging activity including miR125, miR143, miR1272, miR153, miR515-5p, and miR196a-5p (Table 4). In Figure 5A, to demonstrate the proof of principle, we showed how circRNA has_circ_30156 could inhibit type 2 diabetes mellitus (T2DM) by sponging miR143 and preventing miR143 from binding to oxysterol-binding protein-related protein 8 (ORP8) 3 untranslated region (3 UTR).
To elucidate potential downstream pathways for these miRNAs that might influence the LPC-induced inflammatory response, we separately performed IPA for miRNAs sponged by upregulated and downregulated circRNAs, respectively. Notably, the top two pathways for miRNAs sponged by upregulated circRNAs were the pro-inflammatory type 1 T helper cell (Th1), and Th2 activation pathway, and the anti-inflammatory glucocorticoid receptor signaling pathway ( Figure 5B). Similarly, miRNAs sponged by downregulated circRNAs were prominently implicated in the anti-inflammatory transforming growth factor-β (TGF-β) signaling pathway and the pro-inflammatory interleukin-12 (IL-12) signaling and production in macrophages pathway ( Figure 5C). Owing to these numerous conflicting results, the pro-and anti-inflammatory effects of the predicted downstream pathways likely cancel in large part, resulting in little overall pro-or anti-inflammatory effects ( Figure 5D). Specifically, microRNA-125a plays crucial roles both in development and in the adult tissues, including it: (1) contributes to the control of phase transitions in development and/or cell differentiation; (2) regulates the expression of several target proteins that are involved in cell proliferation, apoptosis, and migration; (3) interferes with the expression of the hepatitis B virus surface antigen in liver cells, thus counteracting viral replication . miR-125b plays a significant role in immune response and apoptosis (Potenza and Russo, 2013). The novel long intergenic non-coding RNA UCC promotes colorectal cancer progression by sponging miR-143 (Ribeiro and Sousa, 2014). The function of miR1272 remains unreported. miR153 is tumor suppressive (Huang et al., 2017), and miR153 targeting of KCNQ4 contributes to vascular dysfunction in hypertension (Zhao et al., 2013). miR-515-5p controls cancer cell migration through MARK4 regulation (Carr et al., 2016). Higher expression of miR-196a was both associated with poor overall survival and disease-free survival and recurrence-free survival in different kinds of cancers (Pardo et al., 2016). Taken together, sponging of most of six miRs in ECs may be associated with inhibition of migration, inhibition of vascular dysfunction and immune responses.

74% Significantly Changed CircRNAs Contain Open Reading Frames, the Majority of Which Encode Fragments of Known Proteins
Despite their classification as ncRNAs, numerous circRNAs can be translated into peptides, some of which are novel forms with uncharacterized yet potentially relevant function (Cai et al., 2016;Legnini et al., 2017;Pamudurti et al., 2017). Moreover, the circRNAdb database used to assign circRNA IDs provides ORF peptide sequence data for all circRNAs with ORFs (Chen et al., 2016). To assess the potential of significantly changed circRNAs to encode small proteins/peptides, we performed NIH NCBI BLAST peptide sequence alignments against the NCBI nonredundant protein sequences database (nr), restricted to entries for Homo sapiens, for all significantly changed circRNAs with ORF data in circRNAdb using the NIH NCBI BLAST+ command line utility (blastp -db nr -remote -entrez_query "Homo sapiens [Organism]" -max_target_seqs 1 -outfmt "10 sacc pident e-Value") ( Table 5A). Out of 77 significantly changed circRNAs, 57 (74.0%) contained ORFs; of these, 48 ORF peptide sequences matched substantially (with an expected value < 1 * 10 −4 ) to fragments of known proteins. Notably, 47 of the 48 ORFs  with substantial matches mapped to a protein encoded by mRNA corresponding to the same gene as the circRNA, the sole exception originating from a circRNA corresponding to the non-coding gene PVT1. A majority of significantly changed circRNAs, therefore, have protein-coding potential and could produce fragments of known proteins. Such protein fragments remain both uncharacterized and experimentally undetected, but could theoretically serve physiological functions, including ones different from those of the full canonical protein (Figure 6). Furthermore, using predicted IRES data from circRNAdb, we found that all 57 circRNAs with ORFs contain at least one predicted IRES, and 32 of them (56.1%) contain at least one predicted IRES with an R-value (a metric indicating suitability of a sequence as an IRES) greater than a previously determined cutoff value of 1.54 and a predicted pseudoknot structure, which together indicate a greater probability of ribosomal entry and subsequent translation . We previously reported that this type of IRES drive secondary open reading frames is translated into short proteins (Wu et al., 2009); and that short peptides can disrupt the protein-protein interactions (Xiong et al., 2006a). The Kozak consensus sequence is a sequence which occurs on eukaryotic mRNA and has the consensus (gcc) gccRccAUGG, which plays a major role in the initiation of the translation process (Yang et al., 2005). However, NIH NCBI BLAST searches for Kozak (1996) sequences insignificantly changed circRNAs did not yield any matches. We previously reported that short peptides with high homology to the interaction domain sequences of protein interaction partners could interfere the protein interaction (Yang et al., 2005). To demonstrate the proof of principle, in Figure 6, as suggested by the experimental data in the NIH-NCBI-Gene database 9 , focal adhesion kinase 1 (PTK2) has as many as 225 interaction proteins (Table 5B and Figure 6), thus, the circRNA hsa_circ_23485 encoded short peptide(s) could disrupt 9 https://www.ncbi.nlm.nih.gov/gene/5747 the protein interaction between PTK2 and interaction proteins. Taken together, 74% significantly changed circRNAs contain open reading frames, the majority of which encode fragments of known proteins, which may interfere with protein interactionbased signaling pathways. Of note, in addition to perform experiments in the future to verify the circRNAs, we also notice that various novel computational methods for finding ncRNAs versus coding RNAs (Abbas et al., 2016) and identifying small peptides (Pauli et al., 2015) have been developed, so that we will collaborate in the future with bioinformatics experts to verify the generation of those potential short peptides in HAECs.

DISCUSSION
We previously published experimental data outlining the role of LysoPC (LPC) on HAEC (Li et al., 2018c). In this previous publication we demonstrated that LysoPC stimulation leads to: (1) the upregulation of inflammatory adhesion molecules in HAEC; (2) upregulation of cytokine transcript expression; and (3) promotion of prolonged HAEC activation and expression of innate immune co-stimulation molecule expression. We believe that this RNA-Seq (RNA sequencing) data analyses provide strong experimental evidence for the role of LysoPC stimulation of proinflammatory mRNAs. Given that alternative splicing events lead to the generation of circRNAs from mRNA precursors in other cell types (Wilusz, 2018), there is strong justification to explore the role of LysoPC-induced circRNAs.
Although the importance of circRNAs in regulating EC function has been reported, the comprehensive analyses of circRNAs in EC activation, EC function, and cardiovascular diseases remain at their infancy; and the characterization of circRNAs in atherosclerosis-relevant HAECs has not been reported (Maass et al., 2017;Huang et al., 2018;Shang et al., 2018). To fill in this important knowledge gap, we profiled the circRNAs in HAECs stimulated with we proposed "conditional danger associated molecular  pattern" LPC (Wang X. et al., 2016), and compared the circRNAs significantly modulated in HAECs by LPC with the mRNAs modulated in the same sets of RNA-Seq samples (Li et al., 2018b).
In the current study, we made the following findings: first, LPC modulates the expression of a list of 77 circRNAs in HAECs; and the majority of which, 94.1% of increased circRNAs, compete with related mRNAs for exon usage; second, circRNAs-related genes do not share any top pathways with LPC-induced mRNA pathways, and genes associated with significantly changed circRNAs exhibit little relationship with potential effects of inflammation-related cytokine stimulation; third, 68.8% of the circRNAs feature significant flanking intron homology; fourth, LPC stimulation disrupts the expression levels of 17 spliceosome components, potentially promoting back splicing for circRNA formation over canonical splicing for mRNAs; fifth, shifting chromatin long-range interaction sites from downstream to upstream promotes induction of a list of circRNAs in lysoPC-activated HAECs; and in contrary, shifting chromatin long-range interaction sites from upstream to downstream promotes induction of a list of mRNAs; sixth, six significantly changed circRNAs sponge miRNAs, which may be associated with inhibition of type 2 diabetes mellitus, inhibition of migration, inhibition of vascular dysfunction and immune responses in ECs; and seventh, 74% significantly changed circRNAs contain open reading frames, the majority of which encode fragments of known proteins, suggesting that short proteins generated from the secondary open reading frames potentially driven by internal ribosome entry sites may interfere the protein-protein interactionbased signaling.
As reported previously, LPC stimulation causes significant changes in mRNA expression levels for many genes (Li et al., 2018b). IPA of these changed genes implicates several common upstream regulatory signaling pathways, the most significant of which are the sterol regulatory element binding  transcription factor 2 (SREBF2), SREBF chaperone (SCAP), and SREBF1 pathways (not shown), which jointly serve as key lipid metabolism regulators and contribute to pro-inflammatory signaling (Shao and Espenshade, 2012). This establishes a novel association between upregulated pro-inflammatory lipid metabolism pathways and significant expression changes for 77 circRNAs. Our data are insufficient to conclude any roles of lipid metabolism pathways in circRNA expression modulation. However, given the strength of the association between the significantly changed upstream pathways and the significantly changed circRNAs, more direct future investigation of the effects of lipid metabolism pathway disruption on circRNA levels will be appropriate and could confirm or extend this observation. It is highly unlikely that the extreme ratio of upregulated to downregulated circRNAs or that of upregulated to downregulated corresponding mRNAs among upregulated circRNAs under LPC stimulation would occur by chance, as noted earlier. For this reason, spliceosome depletion, which is experimentally demonstrated to cause such responses in circRNA and mRNA expression (Liang et al., 2017;Wilusz, 2018), continues to be a compelling area for future study despite the inconclusive results presented in this study. For instance, knockdown of downregulated or overexpression of upregulated spliceosome components identified after LPC stimulation, potentially combined with more robust circRNA detection methods such as qPCR verification of back-spliced sequences (Jeck and Sharpless, 2014;Szabo and Salzman, 2016), could clarify the role of each component in LPC-induced disruption of circRNA expression patterns.
Additionally, the dramatic expression changes of many circRNAs under LPC stimulation suggests that they play currently undetermined, possibly significant roles in HAEC activation and the proatherogenic response. Indeed, the circRNA circular antisense non-coding RNA in the INK4 locus (circANRIL) performs anti-atherogenic, anti-proliferative functions in vitro by inhibiting ribosomal biogenesis through assembly factor binding (Holdt et al., 2016), although this circRNA was not detected in our data set. Nevertheless, it is possible that some significantly changed circRNAs fulfill regulatory roles in atherogenesis through similar mechanisms. As circRNAs are better characterized and prediction of RNA-protein binding becomes more sophisticated, it may be possible to detect proteinbinding functions, including those related to atherogenesis or atheroprotection, in our significantly changed circRNAs.
Likewise, many of our significantly changed circRNAs could encode proteins. Although the ORFs detected match the loci FIGURE 6 | One example of potential disruption of protein-protein interactions by circRNA encoded short peptides is demonstrated. of canonically spliced mRNAs, many of the peptides encoded are likely to be fragments of canonical proteins and therefore may serve different functions. Knock-in experiments of these novel predicted peptides could uncover new insights about their potential roles in LPC-induced atherogenesis. Furthermore, as our current findings do not include protein expression data, additional assay experiments under LPC stimulation, such as Western blots and high-performance liquid chromatography-linked mass spectrometer analyses specific to these novel predicted peptides, would provide additional insight on their expression patterns. Given the established coding potential of circRNAs (Legnini et al., 2017;Pamudurti et al., 2017;Yang et al., 2017, such microproteins remain an underexplored pathway for circRNA regulatory roles in atherogenesis. Furthermore, the field of circRNAs is still in its infancy. Software detection methods may give widely differing results, and not all circRNAs predicted from RNA-Seq data can be verified experimentally through methods such as RNase R-catalyzed degradation (Hansen et al., 2016). Simply repeating an LPCinduced HAEC activation experiment, treating with RNase R to remove spurious linear transcripts, and using additional circRNA algorithms to corroborate findings could provide more accurate data on circRNA expression change; using qPCR confirmation for particularly meaningful or significant detections would likewise further improve data reliability. Additionally, recent studies have outlined a number of novel experimental methods for verifying the expression and function (Jeck and Sharpless, 2014;Petkovic and Muller, 2015;Barrett and Salzman, 2016;Werfel et al., 2016;Liu et al., 2017;Zeng et al., 2017;Gomes et al., 2018;Li et al., 2018b). Overall, however, changes in circRNA expression levels following LPC stimulation follow unusual and poorly understood patterns, and further, more intensive profiling under otherwise identical conditions could reveal regulatory pathways not apparent in this study, which uses relatively simple detection methods. To this end, future studies are required using clinical arterial samples from atherosclerotic patients and healthy controls through the collaboration with vascular surgeons using the Temple University Institutional Review Board approved protocol.
Taken together, based on our findings, we propose a novel working model (Figure 7): first, modulation of 77 significantly changed circRNAs in LPC-activated HAECs is temporally associated with LPC-upregulated cholesterol synthesis-SREBP2 pathway and LPC-downregulated TGF-β pathway that we reported recently (Li et al., 2018b); second, the generation of circRNAs in LPC-activated HAECs may be facilitated by three novel molecular mechanisms including: (i) high flanking intron homologies, (ii) LPC-induced modulation of spliceosome components, and (iii) increased upstream chromatin long-range interactions; and third, circRNAs in LPC-activated HAECs may play homeostatic functions by acting on a few novel aspects, including: S (a) competition with related mRNAs for exon usage, (b) miRNA sponge, and (c) secondary open reading frame-generated short peptide disruption of protein-protein interaction. We acknowledge that circRNAs profiling data have brought a huge amount of new findings for detailed experimental verifications in the future. Nevertheless, these novel insights may lead to identification of new therapeutic targets for treating metabolic cardiovascular diseases, inflammations, and cancers.

AUTHOR CONTRIBUTIONS
AL and YS carried out the data gathering, data analysis, and prepared the tables and figures. CD, YL, DY, YZ, XL, CJ, CY, WY, KM, XJ, JS, TR, WH, and HW aided with the analysis of the data. XY supervised the experimental design, data analysis, and manuscript writing. SP aided with the analysis of the data and drafting of the manuscript. All authors read and approved the final manuscript.

FUNDING
This work was supported by NIH grant to XM and HW.