Identification of HIV-1-specific cascaded microRNA-mRNA regulatory relationships by parallel mRNA and microRNA expression profiling with AIDS patients after antiviral treatment.

BACKGROUND: The pathogenesis of human immunodeficiency virus 1 (HIV-1) infection is so complex that have not been clearly defined, despite intensive efforts have been made by many researchers. MicroRNA (miRNA) as regulation factor in various human diseases may influence the course of HIV-1 infection by targeting mRNAs. Thus, studies combining transcription of posttranscriptional miRNA regulation are required.
METHODS: With the purpose of identifying cascaded miRNA-mRNA regulatory relationships related to HIV infection in gene level, the parallel miRNA, and mRNA expression profiles were analyzed to select differential expressed miRNAs and mRNAs. Then, miRNA-mRNA interactions were predicted using 3 data sources and Pearson correlation coefficient was calculated based on the gene expression level for accuracy improvement. Furthermore, the calculation of the regulatory impact factors was conducted to reveal crucial regulators in HIV-1 infection. To give further insight into these transcription factor (TF) regulators, the differentially co-expression analysis was conducted to identify differentially co-expressed links and differential co-expressed genes and the co-expression gene modules were identified using a threshold-based hierarchical clustering method, then modules were combined into a miRNA-TF-mRNA network.
RESULTS: A total of 69,126 differentially co-expressed links and 626 differential co-expressed genes were identified. Functional enrichment analysis indicated that these co-expressed genes were significantly involved in immune response and apoptosis. Moreover, according to regulatory impact factors, 5 most influential TFs and miRNA in HIV-1 infection were identified and miRNA-TF-mRNA regulatory networks were built during the computing process.
CONCLUSIONS: In our study, a set of integrated methods was generated to identify important regulators and miRNA-TF-mRNA interactions. Parallel profiling analysis of the miRNAs and mRNAs expression of HIV/acquired immunodeficiency syndrome (AIDS) patients after antiretroviral therapy indicated that some regulators have wide impact on gene regulation and that these regulatory elements may bear significant implications on the underlying molecular mechanism and pathogenesis of AIDS occurrence.

Introduction

As for the year of 2017, the Joint United Nations Programme on HIV/AIDS (UNAIDS, http://aidsinfo.unaids.org/) estimated that 36.9 million persons were living with human HIV infection worldwide.[ HIV type 1 strain with more infectious and more virulent is the majority of HIV.[ Acquired immunodeficiency syndrome (AIDS), the terminal symptoms caused by human immunodeficiency virus 1 (HIV-1) infection commonly, treated mainly with antiretroviral therapy (ART).[ However, this therapy still has many limitations and side-effects.[ In order to facilitate the development of better treatments, it is extremely imperative to obtain clear understanding of the molecular pathogenesis of gene expression and regulation during viral infections.

MicroRNA (miRNA) and transcription factor (TF)[ both can regulate gene expression by binding to complementary mRNA and manipulating gene transcription, respectively. There is evidence that miRNA may be related to the occurrence of AIDS. For instance, cellular miRNA hsa-miR-29a down-regulates the expression of Nef protein and interferes with HIV-1 replication[; miR-34a enhances HIV-1 replication by targeting PNUTS/PPP1R10.[ Although, previous studies have evaluated the expression of either miRNA or mRNA in cells isolated from subjects with HIV-1.[ Their results are liable to be influenced by individual genetic, viral heterogeneity, and viral load. Especially, the majority of the miRNAs are expressed in a cell- or tissue-specific manner. To date, researches of genome wide mRNA and miRNA profiling with AIDS patients were relatively less.

In the present study isolated RNA and miRNA from peripheral blood mononuclear cells (PBMC) of 3 different groups including HIV-healthy individuals (normal control [CON]), HIV+ ART intervened high viral loads (HVLs) patients and HIV+ ART intervened low viral loads (LVLs) patients.[ The samples were subjected to miRNA and RNA array hybridization. Differentially expressed mRNAs and miRNAs were screened out between HVL, LVL, and control group, than were verified by RT-qPCR. The target genes of the significantly different miRNAs (DEMs) and TFs were predicted. Moreover, functional gene modules were identified by agglomerative hierarchical clustering algorithm and differential co-expression analysis. In addition, the calculation of regulatory impact factor (RIF) was employed to screen regulators (TF and miRNA) which contribute greatly to differential expression of HIV-1 infected patients and healthy individuals.[ Aim to detect novel regulatory relationships, we present a framework to build the complex regulatory network with 3elements: TFs, miRNAs, and their target genes, which reveal the interaction between the 2 different types of regulators in HIV infection course.

Materials and methods

Sample acquisition

The microarray profiles derived from Chengdu University of TCM in Sichuan, China. According to 1994 revised classification system for HIV infection in American,[ all donors in the experimental group were confirmed by Western Blot test for determination of their HIV-1 status. AIDS diagnosis standards were referred to the Guidelines for Diagnosis and Treatment of HIV/AIDS in China (2005).[ The patients were accepted if they conformed to the following criteria: (a) corresponding to the criteria for diagnosis and treatment of AIDS; (b) ranging in age from 20 to 60 years old; (c) agreeing to sign an informed consent form. There were 50 HIV-1 infected AIDS patients and 7 healthy HIV-1 negative donors donated blood samples. In order to explore the relationship between virus load and gene expression level, we screened HVL (n = 14, >50,000 IU/mL) and LVL (n = 16, ≤100 IU/mL) among 50 AIDS patients. Patients’ clinical characteristics were summarized in Table 1. Additional 36 independent blood samples from 33 AIDS patients and 3 healthy donors were applied to verify the accuracy of the differential expression results. All of the patients in our research were in asymptomatic stage and have received ART. PBMCs from morning fasting venous blood samples were isolated by Ficoll gradient separation.[ RNAses Cells were mixed with mirVana RNA lysis buffer (Ambion, Austin, TX). Then the lysates were frozen at –80°C until RNA purification. RNA extraction, hybridization, and microarray images process were performed as previously described.[ Quality control was implemented by measuring the expression of several small RNAs with individual TaqMan miRNA assays (Applied Biosystems, Carlsbad, CA).

Table 1

The clinical characteristics of the 3 groups.

	Unit	CON (n = 7)	LVL (n = 16)	HVL (n = 14)	P value (LVL versus HVL)
HIV-RNA (mean ± SD)	×103 IU/mL	0	0.1	56.6–685 233.46 ± 293.3	.01072
CD3+ & CD4+ (mean ± SD)	cells/μL	955–2860 1684.43 ± 812.3	1–710 178.00 ± 194.4	1–132 34.86 ± 39.7	.01071
PLT (mean ± SD)	×109/L	154–274 200.00 ± 42.5	97–377 259.06 ± 80.6	62–387 204.14 ± 93.9	–
HGB (mean ± SD)	g/L	115–147 130.00 ± 11.7	52–150 107.69 ± 27.5	67–131 102.86 ± 17.3	–
WBC (mean ± SD)	×108/L	45–87 70.00 ± 14.9	23.2–444.6 75.92 ± 100.4	25.1–141 60.73 ± 38.6	–
LYM (mean ± SD)	×107/L	106–336 204.00 ± 94.3	39–231 115.81 ± 67.2	9–117 63.57 ± 32.6	.0113
ALT (mean ± SD)	u/L	12–29 18.57 ± 6.1	6–73 37.25 ± 26.2	7–113 37.64 ± 34.2	–
AST (mean ± SD)	u/L	11–30 18.57 ± 7.3	11–71 32.44 ± 17.2	10–140 41.71 ± 35.9	–

ALT = alanine aminotransferase, AST = aspartate transaminase, CON = normal controls, HGB = hemoglobin, HVL = high viral load, LVL = low viral load, LYM = lymph, PLT = blood platelet, WBC = white blood cell.

Differentially expressed (DE) analysis

MiRNA and mRNA raw data were normalized using package “AgiMicroRna” and “Agi4x44PreProcess” in R v3.0.2. Ultimately, we obtained 499 miRNAs and 14507 genes (see Tables S1 and S2, Supplemental Digital Content, which contain all normalized miRNA and mRNA data). Limma package v3.24[ was utilized to select differentially expressed genes (DEGs) and DEMs between LVL, HVL, and CON groups. Values of fold change, q value (adjusted P value), and cutoff for output were taken into account in analyses.[

Prediction of miRNA targets

We presented a set of integrating methods to identify important regulators and miRNA-TF-mRNA interactions. The bioinformatics pipeline was showed in Figure S1, Supplemental Digital Content (which contains a brief flow chart of analysis strategies). To predict target genes of miRNAs, we used miRBase v21 (http://www.mirbase.org/),[ TargetScan v6.2 (http://www.targetscan.org/),[ and microRNA.org (http://www.microrna.org)[ data sources.

For the purpose of accuracy improvement, Pearson correlation coefficient (PCC) of each miRNAs-mRNA pair was calculated using the expression data, and pairs with PCC < –0.692 (false discovery rate [FDR] < 0.05) were selected. The FDR estimation and PCC threshold calculation methods were displayed in Data S1, Supplemental Digital Content (which contains additional methods) and Figure S2, Supplemental Digital Content (which contains the estimation results of FDR and PCC threshold).

Differential co-expression analysis and prediction of TF targets

DCGL package v2.1.1[ was applied to identify differentially co-expressed genes and links. The genes with lower expression than the median of all were discarded primarily using filtered function in differential co-expression gene and link (DCGL). Then the genes which are significantly more variable than the median gene are retained with P = .25. Differentially co-expressed genes and differentially co-expressed links were identified by “DCp” and “DCe” function and the results were combined by DCsum function. Moreover, TFs and potential targets were predicted by hg19 TfbsConsSites from NCBI and tfbsConsFactors from UCSC Genome Browser (http://genome.ucsc.edu/).[

Regulatory impact factor calculation

Although TFs play an important role in HIV infection, their detection of expression data is limited due to their low and sparse expression.[ RIF is a metric given to each regulator (miRNA or TF) that combines the change in co-expression between the regulator and the DEGs in 2 biological conditions. Hence, RIF is an effective way to screen influential regulator in HIV-1 infection. RIF was calculated according to the following formula:

is the total number of DEGs. and represent the expression of the jth DE gene in conditions 1 and 2; and denote the co-expression correlation of the ith regulator and the jth DE gene in conditions 1 and 2, respectively.[

Identify co-expression gene modules using a threshold-based hierarchical clustering method

To identify the co-expressed gene modules of TFs and miRNAs, threshold-based agglomerative hierarchical clustering analysis was carried out. An N∗N similarity matrix was generated based on PCC calculated above. With threshold T = 0.8, we initialized each gene as a cluster, and found 2 clusters with the shortest distance which means that these 2 clusters have the maximum correlation coefficient. Then, rows and columns of 2 selected clusters were merged into a new cluster and added in matrix. The distance between new merged cluster and each old cluster was recalculated by follow formula: d[(cluster),(picked1, picked2)] = min d[(cluster),(picked1)], d[(cluster),(picked2)]. These steps were repeated until the maximum value of similarity matrix was less than T. The calculation method of threshold and P value was presented in Data S1, Supplemental Digital Content (which contains additional methods) and Figure S3, Supplemental Digital Content (which illustrates threshold and P value calculation of hierarchical clustering method).

Results

Differentially-expressed miRNAs and mRNA in control versus AIDS groups

In order to identify the miRNA and mRNA associated with HIV-1 infection, a parallel genome-wide miRNA and mRNA profiling of 30 AIDS patients and 7 healthy individuals were performed. A total of 1121 upregulated and 1514 downregulated genes, and 58 upregulated and 87 downregulated miRNAs were identified in LVL versus CON (fold-change [FC] > 2 and P < .05). In parallel, 1292 up- and 1534 down-regulated genes, and 50 up- and 78 down-regulated miRNAs were screened from HVL versus CON (Tables 2 and 3). The complete results were contained in Tables S3 and S4, Supplemental Digital Content (which contain DEGs and DEMs between LVL, HVL, and CON). However, the number of both DEGs and DEMs between HVL and LVL was much less than that between DIS (LVL + HVL) and CON (Fig. 1A, C, and Fig. 2A). As presented in Fig. 1C, the heatmap of miRNA between LVL and HVL indicated that there was no significant linear correlation between viral load and the expression level of DEMs. Global mRNA and miRNA expression patterns between control and DIS group were also evaluated by hierarchical clustering (Fig. 1B and Fig. 2B). Patients and control samples were clearly divided into different clusters.

Table 2

Differentially expressed miRNA in 3 groups compared to control group with different statistical significance.

	Differentially expressed miRNA, FC ≥ 2, P < .05			Differentially expressed miRNA, FC < 2, P < .05
	Up-regulated	Down-regulated	Total	Up-regulated	Down-regulated	Total
LVL vs CON	58	87	145	236	202	438
HVL vs CON	50	79	129	229	204	432
LVL vs HVL	1	0	1	17	16	33

Comparative analysis of significantly differentially expressed miRNAs between the different 3 groups.

CON = normal controls, FC = fold change, HVL = high viral load, LVL = low viral load, miRNA = microRNA.

Table 3

Differentially regulated mRNA in low and high viral load groups compared to uninfected control group with different statistical significance.

	Differentially expressed mRNA, FC ≥ 2, P < .05			Differentially expressed mRNA, FC < 2, P < .05
	Up-regulated	Down-regulated	Total	Up-regulated	Down-regulated	Total
LVL vs CON	1121	1514	2635	5922	3464	9386
HVL vs CON	1292	1534	2826	5886	3920	9806
LVL vs HVL	5	19	24	492	742	1234

Comparative analysis of significantly differentially expressed mRNAs between the different 3 groups.

CON = normal controls, FC = fold change, HVL = high viral load, LVL = low viral load.

Figure 1

Expression profile of differentially regulated miRNA in HIV-1 infected versus uninfected samples. (A) The Venn diagram displays the number and overlap of significantly differentially expressed miRNA among the LVL and HVL groups relative to the CON and within the infected groups. (B) Hierarchical clustering using the overlapping miRNA between LVL and HVL groups relative to the CON. MiRNAs in the clustergrams are at a significance cutoff of P < .05 and FC > 2. The clustering is performed on all the samples and the 125 miRNAs expressed in all the samples, using complete linkage and Pearson correlation coefficient (PCC) as a distance metric. Each column represents 1 sample and each row 1 miRNA. The miRNA clustering tree is shown on the left. The color-coded scale (green: expression levels lower than the mean and red: expression level over the mean) for the normalized expression value is indicated at the top of the figure. The clustergrams were generated using gplots package in R software. (C) Hierarchical clustering of miRNA between LVL and HVL. Relative expression levels of miRNAs that were differentially expressed between LVL and HVL at a significance of P < .05 are shown. Columns represent individual patients, and rows represent 1 miRNA. CON = uninfected controls, LVL = low viral load, HIV-1 = human immunodeficiency virus 1, HVL = high viral load, FC = fold change, miRNA = microRNA.

Figure 2

Expression profile of differentially expressed mRNA transcripts in HIV-1 infected versus uninfected subjects. (A) The Venn diagram displays the number and overlap of differentially expressed mRNA among the LVL and HVL groups relative to the CON and within the infected groups. (B) Hierarchical clustering of differentially regulated mRNA between DIS (combined with LVL and HVL) and CON group. The probes in this clustergrams are significantly differentially regulated (P < .05). Red indicates high, green indicates low and black stands for no change in level of expression. Numbers on X-axis represent subject group; Y-axis represents the gene symbol. (C) GO enrichment analysis of the significantly differentially expressed mRNAs between DIS and CON. CON = normal control, GO = gene ontology, HIV-1 = human immunodeficiency virus 1, HVL = high viral load, LVL = low viral load.

For DE mRNA between AIDS patients and controls, we performed Gene Set Enrichment Analysis and the result indicated that most of the core-enriched genes fell into immune response and common disease pathway (Fig. 2C). These results were consistent with our previous research, and further confirmed the significant involvement of these pathways in HIV infection.

To maximize reliability of our results, 7 DEGs and 6 DEMs were randomly selected and measured by qRT-PCR with additional 36 blood samples from 33 AIDS patients and 3 healthy donors (P < .05). Result were consistent with expression profiling and proved the reliability of our DE analysis (see Figure S4, Supplemental Digital Content, which contains qRT-PCR results of 7 DE mRNAs and 6 DE miRNAs in independent samples).

Integrative analysis of putative miRNA-mRNA

The expression of miRNAs and their target mRNAs was negatively correlated. Upregulated miRNAs generally resulted in downregulated target mRNAs and vice versa. Consequently, PCC analysis was implemented to filter miRNA-mRNA pair predicted by database, and then DEGs and DEMs among DIS versus CON were used to screen miRNA-mRNA interactions involved in HIV infection. Nine thousand four hundred fifteen miRNA-mRNA pairs (see Table S5, Supplemental Digital Content, which contains all predicted miRNA-mRNA pairs) were obtained (PCC < –0.692, FDR < 0.05). Subsequently, we investigated the global expression patterns of differential miRNAs and their target genes between DIS and CON. Hierarchical clustering analysis revealed that both miRNAs and mRNAs were divided into 2 sections. In each section the reverse expression pattern of miRNAs and their regulatory mRNAs were shown (Fig. 3A).

Figure 3

The expression patterns and functions of miRNAs and mRNAs. (A) Expression patterns of miRNAs and their regulatory targets using hierarchical clustering analysis. The miRNAs and genes are globally grouped into 2 parts. MiRNAs expression patterns are located above, while mRNAs are below. The expression patterns of miRNAs and targets exhibit the reverse trend. This profile consists of 104 miRNAs and their 1673 target genes. If miRNAs (or mRNAs) are highly expressed in DIS than in CON, the colors are marked in red, otherwise they are marked in green. Enrichment results for pathways analysis of miRNA targets are depicted in B & C. The pathway information was obtained from the KEGG database. (B) Function of up-regulated miRNAs targets. (C) Function of down-regulated miRNAs targets. The orders of biological processes listed in the circle are based on their enriched significance. Red color marks the most enriched functions, whereas purple color corresponds to the lowest enriched ones. CON = normal control, DE = differentially expressed, KEGG = Kyoto Encyclopedia of Genes and Genomes, miRNA = microRNA.

Gene ontology enrichment and Kyoto Encyclopedia of Genes and Genomes pathway analysis for target genes were performed to explore their effect on the progression of HIV infection. Those low-expressed miRNAs, which significantly participate in immune response and apoptosis pathways may have important biological influence in AIDS occurs (Fig. 3B and C). Especially, those miRNAs that involved in calcium signaling pathway and chemokine signaling pathway are crucial for host immune response to virus invasion (adjusted P values < .01). In addition, we found that cell death and apoptosis which associated with infectious complications, substantially regulated by miRNAs. In addition, some miRNAs participated in the cell cycle, cell adhesion, and cell proliferation-related processes, and evidence had shown that the cell adhesion components may relate to HIV infection (Fig. 3C). Analogously, we screened miRNA-mRNA pairs among DEGs and DEMs between LVL and HVL. Gene ontology terms of target genes were mainly relevant to apoptosis and cell death. The target genes specifically participated in the regulation of microtubule cytoskeleton organization, regulation of microtubule-based process, and negative regulation of cytoskeleton organization.

These results suggested that many biological functions which associated with cancers were controlled by DEMs concerned with HIV infection. Novel miRNAs identified in our study significantly control important biological functions, especially cell death and apoptosis process. Interestedly, we also identified several cancer-related pathways as well. This finding is consistent with the fact that HIV virus can trigger cancers or become co-factor of cancers and also compatible with the fact that a number of cancer genes are strongly linked with pro-inflammation response.

Construction of a miRNA-mRNA regulatory network

In general, most miRNAs interact with their target genes in complex cellular networks. Therefore, we integrated dataset of miRNAs and target genes to generate a network that revealed the regulatory relationship from a macroscopic view. The network contained 33 miRNAs and corresponding target genes (Fig. 4). A directed edge from a miRNA to one of its targets states functionally regulated relationship. Important function enrichment pathway of mRNAs were labelled on the network. Metabolism-related sub-network, inflammation and cartilage homeostasis sub-network were interrelated and contribute to cartilage destruction and osteoarthritis development. In the context of virus infection and pathogenesis research it is meaningful to get relevant signaling pathways by mRNA enrichment analysis pathway, such as complement and coagulation cascades calcium signaling pathway and B cell receptor signaling pathway. In addition, several essential cancer-related pathways were discovered as well. MiRNAs could trigger cascade reaction of proinflammatory mediators and cytokines, also miRNAs may interact with chemokines and adhesion molecules during HIV infection.

Figure 4

MiRNA-mRNA correlation networks. MiRNA-mRNA correlation network was rendered in cytoscape. In this network, the miRNAs are demonstrated by diamond shape while mRNAs with ellipse. Upregulated miRNA are in red (A), while downregulated miRNA are in green (B). Upregulated mRNAs are in pink, light red, and pale yellow, while downregulated are in light blue and green, etc. Enrichment results for pathways analysis were added in the network. The pathway information was obtained from the KEGG database. Target genes in a pathway are marked by the same color. KEGG = Kyoto Encyclopedia of Genes and Genomes, miRNA = microRNA.

Construction of TF-mRNA regulatory network

Using DCGL package in R, a total of 69,126 differential co-expression links and 626 differential co-expression genes were obtained (cutoff = 0.5, P = .01) among DIS (HLV + LVL) and CON. The function of differential co-expression genes was predominantly associated with defense response, including immune response, inflammatory response, and wounding response. Meanwhile, on the basis of NCBI and UCSC database, 104 TFs and 166,499 TF-targets interactions were identified (see Table S6, Supplemental Digital Content, which contains putative TF-target pairs). We mapped differential co-expression links to TF-target pairs and allocated the corresponding co-expressed gene to every TF-targeted gene. In terms of each TF, RIF scores were calculated according to the aforementioned formula (1) (Table 4). Subsequently, threshold-based hierarchical clustering method was applied to up-regulated and down-regulated TF-targeted genes separately, which ensured absolute consistency on gene expression pattern. One hundred fifty-six co-expression gene modules (threshold = 0.8, P = 1 × 10–6) were built and partially listed in Table 5. Complete co-expression gene modules were presented in Table S7, Supplemental Digital Content (which contains complete TF-mRNA co-expression gene modules) and the method of threshold calculation for cluster was presented in Data S1, Supplemental Digital Content (which contains additional methods).

Table 4

Co-expression gene modules of some TFs.

TF	Count	Gene module
TCF3	29	NAT14,ATOH7,GPR142,TIAM1,CADM3,GRPR,RORC,RPL3L,DLX2,ACVR1B,COL27A1,CEND1,IRF2BP1,CCR10,EIF4G1,NTN3,EFNA1,DUSP4,CLPS,TMEM145,CLCNKA,PCDHB13,LAMC3,CASKIN2,KRT8,NUMBL,SLC34A3,GFRA3,FAM212A
POU2F1	29	ROBO3,KRT31,CADM3,HCRT,MTSS1L,FAM43B,ADRB3,SPRY3,UBL4B,IRF2BPL,RAPGEF3,FAM181A,PGF,GPR156,KCNQ3,DLX2,GPC4,MYOG,FAM219A,CNTNAP1,CASKIN2,NKX2-5,IRAK1,ZDHHC8,NAT10,GHRHR,GFRA3,B3GALNT1,PACRG
SP1	26	FOXH1,KCNQ3,OLFML2B,SEC14L4,HRC,FAM109B,EIF4G1,TEAD3,ATP6V1G2,RLTPR,NTN3,HRH3,DUSP4,EFNA1,KLF14,MEF2D,IRF2BPL,RTN4RL1,ASB16-AS1,OBSCN,ROBO3,RAPGEF3,LINC01126,CTIF,INHBC,ACTR5
ZEB1	23	COL27A1,KRT14,FAM43B,GPR119,ELF5,CLPS,MTSS1L,MEG3,HMX2,UCKL1,KRTAP2-1,IRF2BPL,CLCNKA,HCRT,KCNQ3,GPR156,BREA2,CASKIN2,LRRN2,KRT8,NUMBL,KLHDC7A,BLVRB
STAT1	22	EFNA1,KLF14,MEF2D,RBMS2,P2RX3,DUSP4,EIF1B-AS1,EMILIN1,PLCXD2,IRF2BP1,NKX2-8,CADM3,RAPGEF3,CTIF,LINC00626,FAM109B,ASB16-AS1,HMX2,NYX,LALBA,FAM219A,SLC34A3
MEF2A	19	MEF2D,KLF14,ZDHHC8,KRTAP4-12,ATP6V1G2,IRAK1,IRF2BPL,CLCNKA,CD244,COL27A1,GRPR,CTIF,LINC01126,HRC,GPC4,FOXH1,KCNQ3,ADORA2A-AS1,LOC100507351
STAT5A	17	KRT14,IRF2BP1,CASKIN2,NKX2-8,MTSS1L,PCDHB13,EIF4G1,NAT14,CIRBP-AS1,INHBC,HMX2,TBATA,CD244,OLFML2B,KCNQ3,GPR156,LALBA
PAX5	17	ATP6V1G2,CLCNKA,OXER1,OBSCN,RTN4RL1,NAT14,CIRBP-AS1,GPC4,FAM181A,FOXH1,BPIFB2,ENHO,EIF4G1,SEC24B-AS1,NKX2-5,EMILIN1,SHISA4
STAT3	16	LALBA,BREA2,LINC00626,KCNQ3,MEF2D,PLCXD2,FAM219A,SLC34A3,SRRD,RORC,NAT14,ADRB3,CRABP1,CD244,KIFC2,IRF2BP1
CEBPA	16	EIF1B-AS1,EMILIN1,CASKIN2,LHB,UBL4A,FAM219A,BREA2,UCKL1,NTN3,ROBO3,MTSS1L,RBMS2,UBL4B,SLC16A11,P2RX3,FAM43B

TF = transcription factors.

Table 5

Top 5 regulators ranked by calculated RIF.

Regulators	Rank	RIF	Description
hsa-miR-23a-3p	1	42.35695674	Developmental process involved in reproduction
hsa-miR-425-5p	2	33.91205464	Extracellular matrix structural protein
hsa-miR-532-5p	3	19.36202123	Cadherin signaling pathway
hsa-miR-1288	4	17.57881395	Apoptosis signaling pathway; Cadherin signaling pathway; FAS signaling pathway
hsa-miR-4685-5p	5	11.39337726	Cadherin signaling pathway; Apoptosis signaling pathway
4303 (FOXO4)	1	23.00992193	Forkhead box protein O4; nucleic acid binding transcription factor activity
6777 (STAT5B)	2	20.48599863	Signal transducer and activator of transcription 5B
2623 (GATA1)	3	19.9980792	Erythroid transcription factor
1050 (CEBPA)	4	18.20150998	CCAAT/enhancer32--binding protein alpha
27022 (FOXD3)	5	15.96317746	Forkhead box protein D3

RIF = regulatory impact factor.

Revelation of integrative regulatory relationships between miRNAs, TFs, and mRNAs

Based on miRNA-mRNA network and TF-mRNA network constructed above, we put forward a “miRNA-TF-mRNA” composite regulatory network (Fig. 5). The framework presents some crucial function module and relevant TF which could give a hint for further understanding the complex mechanism related to HIV infection.

Figure 5

MiRNA-TF-mRNA composite regulatory network. The network was generated by the procedure described in the Results. This network consists of 3 sections. The upper part shows the relationships between miRNA and mRNA. The middle part of the network reveals the regulations between miRNA and TF. Naturally, the network below shows the relationships between TF and mRNA. The larger circles marked with different colors denote various metabolic pathways. The miRNAs are demonstrated by diamond shape while mRNAs with ellipse. The TF is denoted as octagon. An edge represents a regulation from regulators to one of its targets. miRNA = microRNA, TF = transcription factor.

Discussion

Previous bioinformatics analysis shows that over 200 human miRNAs have antiviral effects and potentially interact with virus mRNA.[ A lot of endeavors have been made to explore the transcriptional and posttranscriptional regulatory networks and to detect key transcriptional factors and miRNAs.[ In this study, we designed an integrative strategy to explore miRNA and mRNA profiles. Differential expressed miRNAs which might alter susceptibility to HIV infection were selected primarily. It has been demonstrated that HIV seems to down-regulate some of cellular miRNAs with antiviral effects and our results were accordant with this conclusion. For instance, miRNA-28, miRNA-125b, miRNA-150, and miRNA-223 which are deemed to target the 3’ UTR of HIV-1 transcripts and inhibit viral replication, were all DEMs in our results.[ Sun et al[ have shown that the miR-223 expression level was significantly enriched in HIV-1-infected CD4+ and CD8−, and conversely miR-29a/b, miR-155, and miR-21 levels were significantly reduced. These results were consistent with ours except the miR-155 which was not detected in our research possibly because of the sample difference. Additionally, we observed that HVL patients exhibit an expression profile similar to LVL patients, and the limited difference possibly is a consequence of the fact that our analyses were performed on total PBMC rather than enriched cell subsets, which could dilute the effects caused by HIV. In brief, although the mechanism is not completely revealed, there is evidence suggesting that miRNA expression is actually being regulated by influence HIV-1 virus. Comes to regulatory module construction, at least 3 factors could contribute to the strict selection of regulatory modules in our predictions. Firstly, we demanded that all mRNA in the same TF regulated module must be highly correlated to each other. Secondly, those highly correlated mRNAs were DE between DIS and CON samples. Thirdly, a method similar to hypergeometric distribution was utilized to evaluate the significant of a candidate module in our study, and apparently modules discerned in this manner would have better P values.

To focus on the most meaningful information, we calculated RIF for each miRNA and TF. Among DIS (HLV + LVL) against CON samples, 5 out of the top 10 TF and miRNA have confirmed to be related to in HIV infection (Table 4). The top 1, AFX1, encodes a member of the O class of winged helix/forkhead TF family. TFs of this class involve in the regulation of the insulin signaling pathway, negative regulation of the cell cycle and increased of proteasome activity in embryonic stem cells.[ Moreover, with the highest score among miRNA, hsa-miR-23a-3p were found to be related to chronic inflammatory disease and cancers like psoriasis, lung cancer, and colon cancer in previous literature, which implied that it may play a role in immune reaction after HIV-1 infection.[ In our study, upregulated hsa-miR-23a-3p has 2 downregulated target gene: CXCL12 (6387), FKBP4 (2288). The chemokine CXCL12 specifically binds to the signaling receptor CXCR4, which is one of the co-receptors for HIV-1 virus infection of immune cells. CXCL12 could compete with the virus to bind to the CXCR4 receptors and act as potent blocker of HIV infection.[ FKBP4 encodes a cis-trans peptidyl-prolyl isomerase that binds to the immunosuppressants FK506, making FK506 no immunosuppressive activity. Also, peptidyl-prolyl isomerase interferes with IRF-4 DNA binding and transcription activity by inducing conformational changes, and IRF4 deficiency leads to CD4+ T cells dysfunction.[Figure 6 shows the RIF score of TFs and miRNAs. It indicated that there is no significant difference in TFs and miRNAs, and also roughly implied that gene expression is almost equivalently effected by transcriptional and posttranscriptional regulators.

Figure 6

Scatter diagrams of regulatory impact factors (RIF). Scatter plot of the RIF in contrastive groups: (A) DIS versus CON; (B) HVL versus CON; (C) LVL versus HVL; (D) LVL versus CON. CON = normal control, HVL = high viral load, LVL = low viral load.

Aim to identify meaningful regulators involved in the occurrence of AIDS, we have developed a set of methods to identify HIV-associated miRNA-mRNA and TF-mRNA regulatory modules from mRNA and miRNA microarray data. Based on the above regulatory modules, a miRNA-TF co-regulatory network that marked important regulators and genes was constructed. MiRNA/mRNA expression profile data revealed that some regulators have a wide impact on gene regulation in AIDS patients. However, this study has certain limitations. First of all, we ignored the influence of ART duration on the regulatory relationship, and did not record the treatment times in clinical characteristics. Furthermore, with all PBMCs as the research object, we cannot detect the heterogeneity of the regulatory network at the single-cell level. There is limited understanding of gene transcriptional regulation, and also regulatory networks are not static. Although specific and complete interactive relationship remains uncertain, advanced studies for these mechanisms will hold clue for explicit understanding of complex genomic interactions, and also will make preparations for the new therapeutics to treat AIDS with human. Consequently, our study provided some new insights into the host responses and AIDS pathogenesis, which may be applied to the study of miRNA regulation in other complicated human diseases.

Author contributions

XL and QF conceived and designed the study. XC and YFL performed the main experiments and analysis. GG and JY participated in data analysis. YFL, CC, and GF collected blood samples with AIDS patients. YSL and XZ equivalently to network graph construction. MY managed the samples. XC and YFL wrote the paper. All authors discussed the results and revised the manuscript.

Conceptualization: Jian Yang, Quansheng Feng, Xiao Li.

Data curation: Fangyuan Shen, Yefang Liu.

Formal analysis: Xiaoqiang Chai.

Funding acquisition: Quansheng Feng, Xiao Li.

Methodology: Xiaoqiang Chai, Jian Yang, Quansheng Feng, Xiao Li.

Software: Fangyuan Shen, Lanchun Wang.

Validation: Lanchun Wang, Xiaoqiang Chai.

Writing – original draft: Fangyuan Shen, Xiaoqiang Chai.