Molecular characterization of atherosclerosis in HIV positive persons

People living with HIV are at higher risk of atherosclerosis (AS). The pathogenesis of this risk is not fully understood. To assess the regulatory networks involved in AS we sequenced mRNA of the peripheral blood mononuclear cells (PBMCs) and measured cytokine and chemokine levels in the plasma of 13 persons living with HIV and 12 matched HIV-negative persons with and without AS. microRNAs (miRNAs) are known to play a role in HIV infection and may modulate gene regulation to drive AS. Hence, we further assessed miRNA expression in PBMCs of a subset of 12 HIV+ people with and without atherosclerosis. We identified 12 miRNAs differentially expressed between HIV+ AS+ and HIV+ , and validated 5 of those by RT-qPCR. While a few of these miRNAs have been implicated in HIV and atherosclerosis, others are novel. Integrating miRNA measurements with mRNA, we identified 27 target genes including SLC4A7, a critical sodium and bicarbonate transporter, that are potentially dysregulated during atherosclerosis. Additionally, we uncovered that levels of plasma cytokines were associated with transcription factor activity and miRNA expression in PBMCs. For example, BACH2 activity was associated with IL-1β, IL-15, and MIP-1α. IP10 and TNFα levels were associated with miR-124-3p. Finally, integration of all data types into a single network revealed increased importance of miRNAs in network regulation of the HIV+ group in contrast with increased importance of cytokines in the HIV+ AS+ group.

Studies of human immunology highlighting variability across individuals have revealed that variation in functional responses is largely determined by non-heritable factors, most likely by environmental exposures [1][2][3] . Viral infections such as human immunodeficiency virus type-1 (HIV) interfere with cellular targets and cause derangement of the immune response. Combination anti-retroviral therapy (cART) has rendered HIV infection as a chronic manageable disease. However, persons living with HIV remain at increased risk of co-morbidities despite cART; the most frequent of these is cardiovascular disease (CVD) driven by atherosclerosis (AS) [4][5][6] . It is likely that cART and HIV itself play a role in AS 7 . This observation was also reported in a study that assessed carotid artery thickening, a marker of AS, which was up to 24% higher in persons living with HIV compared with uninfected sex-and age-matched persons, and was comparable in cART-naïve persons and with virologic suppression on cART 8 . Moreover, a large retrospective analysis indicates that HIV infection status independently confers an odds ratio for acute myocardial infarction of ~ 1.93 (95% CI, 1.21-2.93) after adjusting for traditional risk factors such as age and hypertension 9 . Other more recent studies also support this observation [9][10][11] . Thus, persons living with HIV are at higher risk of AS even with ongoing cART therapy 8,12,13 . Moreover, elevated AS prevalence in HIV-1 elite controllers indicates that inflammation can be induced by the presence of HIV infection 14 .
AS is triggered by a pro-thrombotic inflammatory environment that attracts platelets and monocytes to injured endothelial cells. As monocytes migrate into the vessel wall, they differentiate into macrophages. plex Human Cytokine/Chemokine panel kit (Millipore Cat# HCYTMAG-60K-PX29) was used to quantitate the cytokine and chemokine levels using Luminex magnetic microbead array technology. The assay was performed according to the manufacturer's instructions. In brief, plasma specimens were thawed on ice and micro-centrifuged for 5 min at 15,000 RCF at 4 °C to remove any particulate matter. 25 µLs of reference standard dilutions, control, and test specimens were incubated overnight at 4 °C with shaking at 800 RPM in duplicate wells on a 96 well plate with 25 µLs of either assay buffer for specimens or plasma matrix for standards and controls and 25 µLs of microbead solution. Wells were washed 3 times with wash buffer on an automated magnetic plate washer (Bio-Rad Bio-Plex ProII) and 25 µLs of biotinylated detection antibodies were added for 1 h at room tempera-  The trimmed means of the fluorescence intensities across the  beads for a given sample and target were exported from the Luminex Xponent software. The geometric mean of  two technical replicates for each sample were used in downstream analyses. The data for the cytokine/chemokine  panel is available in Supplemental Table 2.

Levels of pro-and anti-inflammatory mediators in peripheral blood plasma.
Peripheral blood mononuclear cell isolation and RNA-sequencing. mRNA and micro-RNA were isolated and sequenced from PBMCs at the UR Genomics Research Center using their standard protocols. mRNA libraries were prepared with the Illumina TruSeq Kit v2, using oligo-dT beads for poly-A selection. mRNA was sequenced as single-ended 100 bp reads, yielding an average of 26.7 million reads per sample, for a total of 25 samples, representing four sample groups (normal control, n = 6; AS+, n = 6; HIV+, n = 6; HIV+ AS+, n = 7). Reads were trimmed with Trimmomatic to remove low-quality read tails or entire reads if the remaining post-trim read length was too short. Genome mapping was performed with STAR 2.5.2b to Human genome GRCh38.p7 with annotation from Gencode Genes version 25. Read counting was performed with the feature-Counts from Subread. Count data is available in Supplemental Table 3. Small RNA was size-selected from total RNA and libraries were prepared with the Illumina TruSeq Small RNA kit. Small RNA was sequenced as singleended 50 bp reads, with an average of 17 million reads per sample, for a total of 12 samples representing two sample groups (n = 6 each for HIV+ and HIV+ AS+). Reads were trimmed, mapped, annotated, and counted with the miRge pipeline 55 . Count data is available in Supplemental Table 4. miRNA were mapped to cell types using an available catalogue of miR expression in blood cell types 56 .
RNA-sequencing analysis. Small RNA-seq and mRNA-seq datasets were read into R (version 3.5.1) 57 , and filtered to remove sequences with low expression based on raw counts. For mRNA, any gene with fewer than 10 reads in all samples was excluded, leaving 19,861 genes for annotated gene types which may be polyadenylated (protein coding, lincRNA, and processed pseudogenes). For miRNA, any miRNA with fewer than 10 reads in all samples was excluded, yielding a set of 778 miRNAs used for downstream analyses. Exploratory analysis of each dataset (e.g. clustering, expression heatmaps, PCA) was performed on variance-stabilized transformed (VST) counts, as computed by the DESeq2 package (version 1.22.2) 58 . Functional and pathway enrichment was investigated with gProfileR (version 0.6.7) 59 . Differential expression analysis was performed with DESeq2 58 for both mRNA and miRNA (independently). The results of differential gene expression analysis for mRNA for all genes in the dataset are available in Supplemental Table 5, and differentially expressed miRNAs are available in Supplemental Table 6. For mRNA, groupwise comparisons were performed a corrected (FDR) p-value and foldchange cutoffs are noted in the results. Predicted miRNA target genes were found for differentially expressed miRNAs with the miRNAtap package (version 1.16.0) 60 which incorporates predictions from five algorithms: DIANA, Miranda, Targetscan, PicTar, and miRDB. A given miR-target relationship was recorded if it was supported by at least two of the five prediction algorithms. The predicted targets of differentially expressed miRNAs are available in Supplemental Table 7, while functional enrichment analysis on these targets is available in Supplemental Table 8. miRNA expression across normal human blood cells. Publicly available miRNA expression data across normal human blood cells were downloaded from GEO (GSE100467) 56 , and read into R with the BioConductor package IsoMirs (version 1.10.1) 61 to obtain read counts for each of the 450 samples, representing ten cell populations from 162 unique donors. A matrix of raw counts was assembled by joining the individual sample results with dplyr, and miRNAs with very low counts (< 3) across all samples were excluded, leaving 1062 unique miRNAs with expression in at least one sample. Variance-stabilized (VST) and library-normalized counts were obtained with DESeq2. To obtain a summary of expression of each miRNA for each cell type, miRNA abundances were discretized in 5 levels, with higher-value bins representing higher expression. The number of bins was determined by ceiling 3 √ max(n) , where n is the number of samples for each cell type 62 . Relative levels of miRNA expression in each cell-type was given by the median of the discretized values.

Data integration analysis.
Three types of data integration were performed; (1) miRNA and mRNA (2) cytokine and mRNA and (3) cytokine, mRNA and miRNA datasets were combined. Predicted miRNA targets were further filtered based on their significant negative correlation with mRNA (p < 0.05, r ≤ −0.5 ). Linear discriminant analysis (LDA) was performed on predicted miRNA target genes to identify genes that distinguish between HIV+ AS− and HIV+ AS+ samples. To identify transcriptional changes associated with cytokine levels, a single person transcription factor score (spTFscore) was calculated. Particularly, transcription factor (TF) target genes were predicted using positional weight matrices from JASPAR 63 as described 64,65 . To calculate spTFscore, PCA was performed on TF target genes and m i=2 s 2 i * v i was calculated from i = 2 to m , the PC at which 50% cumulative variance is explained (including PC 1 ), s is the unweighted PC score, and v is the % variance explained by the PC. Since PC 1 was associated with inter-individual variability in overall transcription levels rather than inter-group and transcription factor specific variability, the score was calculated starting with PC 2 . Thus, the spTFscore captures the contribution of each TF measured by variance of their target gene expression in each person. Pearson correlation was measured between all plasma cytokines and spTFscore. TF target gene weights were calculated as contributions to spTFscore by where w i is the loading of the gene in that www.nature.com/scientificreports/ PC, v is the variance of each PC, v 1 is the variance explained by the first PC and m is the PC at which 50% cumulative variance is explained, as with the spTFscore. In this way, the gene TF weights reflect the degree to which each gene contributes to the spTFscore for that TF. Finally, to integrate miRNA, mRNA, and cytokine datasets an integrated multi-omic network was constructed and differential network analysis performed using xMWAS 66 . In order to include all data types, only samples from HIV+ AS− and HIV+ AS+ groups were included. Briefly, xMWAS constructs a matrix for pairwise correlation analysis using sparse partial least squares then ranks and filters the top association scores by p-value and association score. The edge lists are merged to a global network upon which community detection algorithms are utilized to group nodes into communities. Further, centrality scores are calculated for each node in each group separately, then compared to generate differences between groups. These centrality scores are available in Supplemental Table 9. Graphical representations of these networks were constructed in Cytoscape (version 3.7.0) 67 .

Validation of miRNAs.
In order to validate RNA-sequencing results, levels of miRNAs 144-3p, 144-5p, 183-5p, 451a, and 4732-3p were measured in a subset of samples used for RNA-seq (n = 2 HIV+ AS+ , and n = 2, HIV+ AS−). In addition, we also measured the expression of these miRNAs in new samples (not used in RNAseq, n = 3 HIV+ AS+ and n = 3 HIV+ AS−). RNA from new samples was extracted using miRNeasy mini kit (Cat # 217004, Qiagen). RT-qPCR was performed as per the manufacturer's protocol (Taqman microRNA assay kits, Cat # 4427975 and 4440886) using the Biorad CFX Connect real time PCR machine. miRNA U6 was used as an endogenous control. The data is shown as ΔCt values. The results were analyzed by unpaired t-test using Graphpad Prism software.

Results and discussion:
Gene signatures of AS in persons living with HIV. Differential expression analysis identified seventeen genes in each of the two comparisons, HIV vs HIV+ AS+ and HIV vs control group (p-adj < 0.05). There was one common genes among the two contrasts, integrin alpha D (ITGAD/CD11d), which is strongly associated with inflammation 68,69 . Comparison between the HIV+ AS+ and control groups revealed five genes higher in HIV+ AS+ group: CD8A, CD8B, NPDC1, JAKMIP1, and STYK1; with five higher in control: IL17RE, SPTSSB, TPBG, OVOS, COL13A1. The presence of CD8A and CD8B provides strong evidence that there is an increase in CD8 cell activity in HIV+ AS+ group. Previous reports have shown memory CX3CR1+ CD8(+) cells may contribute to atherosclerosis in persons living with HIV 70 . On the other hand, downregulation of the IL-17 receptor (IL17RE) could be in response to increased immune stimulation. However, large variation in gene expression across individuals within groups led to identification of few differentially expressed (DE) genes.
To investigate transcriptional changes in this highly variable data, we developed a single person transcription factor score (spTFscore) that could be compared with cytokine levels to identify transcription factors regulated by cytokines. The correlations between cytokine levels and the spTFscore revealed regulation by different arms of the immune system. Specifically, Th1, Th17 signature cytokines along with TNFα, MIP-1α and GCSF were associated with BACH2. It is known that Bach2-Batf interactions are required to prevent an excessive Th2 response 71 . IP10 was associated with several transcription factors such as IRF1 and CEBPα. Further analysis indicated that predicted targets of these transcription factors have little overlap, suggesting diverse transcriptional programs potentially associated with IP-10. IL-13, IL-8, TNFα, IL-1R and GMCSF were all associated with ZBTB7C, a TF that represses matrix metalloproteinases (MMPs) (Fig. 1). A number of the cytokines in this list are important modulators of atherosclerosis, including IL-15 which stimulates T cells and is known to be expressed in atherosclerotic lesions [72][73][74][75] .
Many target genes of TFs significantly associated with the cytokines were differentially expressed (p-adj < 0.1 and LFC > 1.0) between HIV+ vs HIV+ AS+ group and HIV+ vs healthy (Fig. 2). The genes differentially expressed between HIV+ vs healthy and targets of TFs associated with IP10 include genes implicated in HIV infection e.g. CCR4, CCR6, and KIF19. CCR4 + and CCR6 + CD4 + T cells are highly permissive to HIV replication 76 . Interestingly, IP-10 enhances the recruitment of HIV-target cells and is a biomarker of early disease onset 77 . Eotaxin, a cytokine that was significantly different between HIV+ and control groups was also associated with several TFs, including GRHL12, PBX1, EOMES, and TBR1. Levels of circulating LAG3, an Eomes target gene, have been associated with HDL cholesterol and risk of coronary artery disease 78 . LAG3 here had increased expression in HIV+ group compared to control.
We also evaluated additional genes different at p-adj < 0.1 and LFC > 1.0 and associated with either HIV infection or AS (Supplementary Table 5). ATP1A2 is an Na + /K + pump component involved in smooth muscle contraction, and it is differentially expressed in vasculature between myocardial infarction (MI) and stable angina 79 , and is also differentially expressed in our data between HIV+ and HIV+ AS+ . SLIT2 is a known inhibitor of HIV transmission 80 that was significantly elevated in HIV+ vs control. GPR15 is a co-receptor for HIV and a possible biomarker for smoking status 81,82 ; it had increased expression in HIV+ AS+ than HIV+ group. In this context, it may be a marker for the deleterious interaction between HIV and smoking status, which leads to atherosclerosis. VCAM1, a critical pro-atherosclerotic factor expressed on the endothelium, was upregulated in HIV+ vs control.
In conclusion, transcriptional changes and their relation with cytokine levels reveals several dysregulated interactions contributing to AS in the context of HIV and suggests further risk in HIV+ persons.
Characterization of micro-RNAs in HIV+ persons with carotid plaques. We performed genomewide miRNA sequencing in HIV+ persons with (n = 6) and without (n = 6) carotid plaques. All 12 participants were males of ≥ 50 age to minimize the inter-person variation and maximize the differences across disease/ infection groups in this small study. Moreover, in this focused group confounding variables were balanced, increasing the statistical power. Sequencing yielded at least 15 million short reads per sample, and most reads www.nature.com/scientificreports/ were ~ 22 nt after trimming, the typical length of miRNAs, indicating successful size selection and good read quality. Less than 10% of reads remained unmapped in all cases. miRNAs were identified and quantified from sequencing data with the miRge pipeline 55 . 778 miRNAs were considered after filtering out miRNAs that did not have at least 10 reads in all samples. Differential expression analysis identified 12 microRNAs with significant expression changes, out of which 3 were higher in HIV+ subjects with AS and 9 were higher in HIV+ subjects without AS (Fig. 3A). A subset of these differently expressed miRNAs that have been previously implicated in cardiovascular disease are discussed below. High miR-124-3p indicated susceptibility to atherosclerosis in a population of smokers 83 . In our dataset, miR-124-3p had no or very low expression in HIV+ AS− samples, but was found in 5 of 6 HIV+ AS+ samples (only two of which were identified as current smokers) indicating that the elevated level of miR-124-3p may be indicative of AS risk outside of the context of smoking. Another miRNA, miR-182-5p is found to be upregulated about two times in myocardial infarction vs controls in whole blood samples, however this miRNA was downregulated in HIV+ persons with AS compared to without 84 . miR-34a-5p is expressed in T-cells (CD3+ and CD4+), and is negatively correlated with T-cell counts after HIV infection, possibly due to lysis of T-cells 85,86 . miR-34a-5p has also been shown to be increased in heart failure and decreased in peripheral arterial disease [87][88][89] , in addition to being upregulated in endothelial progenitor cells from coronary artery disease patients 90 . In ApoE-deficient mice fed a "western-style" diet for 16 weeks, the mice developed severe atherosclerotic lesions and exhibited progressive upregulation of miR-34a-5p 91 . In accordance with these observations linking higher levels of miR-34a-5p to AS, we find this miRNA to be elevated in most of our HIV+ AS+ samples compared to HIV+ alone. Activity of the TFs was measured by single person transcription factor (spTF) scores. spTF was calculated as described in the methods. Briefly, TF targets as described in 64 were used to perform PCA and subject specific PC scores were used to calculate pearson correlation (ranging from low to high depicted by white to red color bar) with cytokine levels. The significant correlations (p < 0.05) are indicated with stars. Analysis and graphical output done using in R (ver 3.5.1) 57 .  PRY2  PRY  LCN2  UCN3  OBP2A  ACTC1  SCUBE3  MUC5B  MTUS1  RIMBP2  GSTM3  SLC4A10  ANKRD18B  CCR4  RPS4Y2  RPL10L  RPL39L  NACA2  B3GALT2  SPTSSB  SIAH3  TMSB4Y  RASSF6  RORC  ISM1  ANKRD18A  BAALC  HES4  SEMA6B  ULBP3  KIF19  HBD  HTRA3  FAT4  SFRP1  MSC  COL13A1  RCAN2  LAG3  THNSL2  CCNJL  PXDN  GLB1L2  DZIP1L  STYK1  DAPK2  LIMCH1  PTMS  C8G  RNF223  ACTA1  AHNAK2  CCR6  SPOCK3  APOL4  MTRNR2L3  MTRNR2L10  ATP2A1  PRSS35  MTRNR2L12  MTRNR2L1  GZMH  ENAH  TMEM132B  CDC25A  SYN2  SPTLC3  ETV7  KLHL23  MTRNR2L7  MTRNR2L5  ALPK2  SLC7A8  POTEJ  MTRNR2L8  ATP6V0A4  MTRNR2L4  DNAH2  www.nature.com/scientificreports/ To validate the miRNAs identified by RNA-seq we performed RT-qPCR on 5 miRNAs in ten participants (Fig. 3C). The miRNAs miR-144-3p, miR-144-5p, miR-183-5p, miR-451a and miR-4732-3p were chosen based on the largest fold change between HIV+ AS+ vs HIV+ AS−. The samples were run in pairs as shown in Fig. 3C. All of the 5 miRNAs were lower in HIV+ AS+ compared to HIV+ AS− in RNA-seq experiments. RT-qPCR validated this observation when a subset of the same samples from RNA-seq were used. Among the new independent samples one pair showed an opposite direction, whereas all other samples validated the observations from RNA-seq. The five validated miRNAs are novel and very little is known with respect to their role in HIV infection and atherosclerosis. miR-144-3p 92,93 and miR-183-5p 94 have been linked atherosclerosis, cholesterol metabolism, and HIV infection in different organisms and model systems. In consensus with our observations miR-451a has been found to be reduced in circulation of subjects with coronary heart disease 95 .
To investigate which specific cell-types are known to express the differentially expressed miRNAs from our study, we used a public compendium of miRNA expression across human blood cell types 56 . Specifically, expression of differentially expressed miRNAs was investigated in the public compendium to determine the cell types that normally express the differentially expressed miRNAs (Supplemental Fig. 1). This analysis revealed that only miR-486-5p was ubiquitously expressed, and that the other DE miRNA may exhibit low or no expression in normal CD8+ or CD4+ T cells and CD14+ monocytes (Supplemental Fig. 1). To further characterize the effect of these miRNAs, we performed pathway analysis of their predicted targets (Fig. 3B). Importantly, similar to previous findings, miR-124-3p was enriched in 'Signaling to ERKs' and ' Activation of AMPA Receptors' . miR-183-5p was associated with a number of pathways including the 'Twik related potassium channel (TREK)' , 'Signaling www.nature.com/scientificreports/ by Receptor Tyrosine Kinases' , and 'ROBO receptors bind AKAP5′. Finally, miR-182-5p was associated with 'DAG and IP3 Signaling' , 'Signaling by Rho GTPases' , and 'FCγ receptor (FCGR) dependent phagocytosis' 96,97 . Thus, most of our findings corroborate with previously known roles of miRNAs in atherosclerosis, and reveal novel insights about development of AS in persons living with HIV. However, most previous studies have been performed with different model systems, experimental designs, and protocols, necessitating extended study in people living with HIV.
Role of miRNAs in regulatory networks. miRNAs are important regulators of mRNA turnover and are thought to regulate 10-30% of mRNAs [98][99][100] . miRNAs function post-transcriptionally via the RNA-induced silencing complex (RISC) to inhibit mRNA. Hence, we found predicted targets of differentially expressed miR-NAs that also had significant negative correlation with mRNA expression in our data. A set of 27 of these genes with highest loadings from Linear Discriminant Analysis (LDA) was sufficient to differentiate the HIV+ AS+ and HIV+ AS− groups (Fig. 4A). This list included genes such as SPRED1 which is involved in IL-15 and FGFR1 signaling, GRK6 which is involved in the myometrial relaxation and contraction pathway and calcium regulation in cardiac cells, and RAN which is involved in the export of viral ribonucleoprotein. Also in this set of 27 genes is SLC4A7, a critical sodium and bicarbonate transporter involved in generation of Nitric Oxide (NO) signaling and blood pressure regulation; NO signaling has been linked to AS risk 101 . Importantly, 15/27 of these genes are upregulated in HIV+ AS+ vs HIV but downregulated in AS vs healthy control with no correlation between the two changes, indicating that these genes describe an HIV-specific atherosclerosis state (Fig. 4B).
miRNAs can be transcriptionally regulated by cytokine signaling 102 while cytokine production 103 and signal transduction can be influenced by miRNAs 104,105 . Hence, we asked which cytokines are significantly associated with DE miRNAs. IL-8 was positively associated with most of the miRNAs but negatively associated with miR-124-3p (Fig. 4C). miR-124-3p, which regulates ERK signaling and activation of AMPA receptors, was also negatively associated with MIP-1α but positively associated with TNFα and IP10 (also known as CXCL10). IP10 is induced during antiviral response and was higher in HIV+ and HIV+ AS+ groups (not statistically significant, p = 0.08). Numerous studies have reported abnormally high plasma IP-10 levels in the context of HIV infection, and IP10 is considered an important pro-inflammatory factor in the HIV disease process 106 . IP-10 and TNFα are both associated with miR-124-3p (Fig. 4C) and are highly expressed in HIV+ AS+ group. As discussed before, we found IP10 to be associated with diverse transcriptional programs activating several genes implicated in HIV infection. The role of miR-124-3p in this regulation needs to be further investigated. IL8 has interactions with miRNA in a number of diseases and conditions 107 . Other miRNAs have been found to de-repress expression of IL6 and IL8 in a cell-based model of sepsis 108 . These results suggest novel regulatory relationships between miRNAs and IL-8; however, IL-8 levels were not different in HIV+ persons with and without AS. IL-8 expression (at both the gene and plasma levels) has been reported to be significantly elevated in HIV-infected children that were not yet receiving ART therapy, and decreased after therapy began 109 . Several additional associations were found between miRNAs and cytokines (Fig. 4C). Thus, this work suggests a role for pro-inflammatory cytokines in regulating miRNAs which can in turn regulate gene programs.
Integrated network of cytokines, miRNA and mRNA. Finally, having established several links between cytokines, miRNAs, and mRNAs in a pairwise manner, we set out to investigate such interactions systematically. Hence, we constructed a network, incorporating all three data types-mRNA, miRNA, and cytokine profiles, for HIV+ and HIV+ AS+ groups. The network includes 175 genes, 23 cytokines and 10 miRNAs (Fig. 5, Supplementary Table 1). Differential network analysis using eigenvector centrality identify nodes with higher importance with AS in the context of HIV.
Higher centrality of cytokines and chemokines compared to the genes was expected because a small number of cytokines measured were connected to many miRNA and genes. Three cytokines, IL1β, IL15, and MIP1α, which were also associated with BACH2 were more central in HIV+ AS+ than HIV suggesting their role in regulating miRNA and mRNA expression in persons living with HIV and atherosclerosis. IL13, IL5, TNFβ, and IL1-RA associated with ZBTB7C, were more central in HIV than HIV+ AS+ groups, along with hsa-miR-144-5p, hsa-miR-451a, and hsa-miR-182-5p. hsa-miR-451, one of the miRNA independently validated in this study, has been shown in rats to target the IL6 receptor in arterial smooth muscle cells 110 , and IL-6 knockdown in mice induced release of cytokines including IL6, TNF, CCL5/RANTES, and CCL3/MIP1α 111 . This is corroborated by increased centrality of these downstream cytokines in the HIV+ AS+ group. hsa-miR-182-5p, on the other hand, is a pro-atherosclerotic AKT1-inhibitor implicated in acute coronary disease (atherosclerosis of the heart) 112 via endothelial cell dysfunction. Positive association between pro-AS miRNAs and anti-inflammatory/anti-AS cytokines such as IL5 113 might suggest a compensatory response that delays AS. Thus, our network analysis identifies two critical atherosclerosis-related miRNAS− hsa-miR-182-5p which affects endothelial cells, and hsa-miR-144-5p (also validated by RT-qPCR) which affects macrophages, known to be elevated in high viral titer HIV infection.

Conclusion
Taken together, this study revealed novel miRNAs differentially expressed in AS+ /-people living with HIV, and identified many interesting interactions of miRNAs with cytokines. Centrality measures indicate that miRNAs display importance in regulating molecular networks in people living with HIV, whereas cytokines play a greater role in people living with HIV and AS. Though mRNA measurements were highly variable even with very strict subject selection criteria, novel computational approaches such as single person transcription factor score, and leveraging public domain datasets, provided insights about regulatory factors. Nevertheless, several findings in www.nature.com/scientificreports/  www.nature.com/scientificreports/