Abstract

   Long-term non-progressors (LTNP) and elite controllers (EC) represent
   spontaneous natural models of efficient HIV-1 response in the absence
   of treatment. The main purposes of this work are to describe the
   miRNome of HIV-1 infected patients with different extreme phenotypes
   and identify potentially altered pathways regulated by differentially
   expressed (DE) miRNAs. The miRNomes from peripheral blood mononuclear
   cells (PBMCs) of dual phenotype EC-LTNP or LTNP with detectable viremia
   and HIV-infected patients with typical progression before and after
   treatment, were obtained through miRNA-Seq and compared among them. The
   administration of treatment produces 18 DE miRNAs in typical
   progressors. LTNP condition shows 14 DE miRNA when compared to typical
   progressors, allowing LTNP phenotype differentiation. A set of four
   miRNAs: miR-144-3p, miR-18a-5p, miR-451a, and miR-324 is strongly
   downregulated in LTNP and related to protein regulation as AKT, mTOR,
   ERK or IKK, involved in immune response pathways. Deregulation of 28
   miRNA is observed between EC-LTNP and viremic-LTNP, including
   previously described anti-HIV miRNAs: miR-29a, associated with LTNP
   phenotype, and miR-155, targeting different pre-integration complexes
   such as ADAM10 and TNPO3. A holistic perspective of the changes
   observed in the miRNome of patients with different phenotypes of
   HIV-control and non-progression is provided.

   Keywords: HIV-1 infection, long-term non-progressors, elite
   controllers, miRNA-Seq, miRNome, differential expression, biomarker
   discovery

1. Introduction

   HIV natural infection is characterized by continued viral replication,
   systemic chronic immune activation, and a strong decline of CD4+ T
   cells. In most patients, AIDS is developed within a few years from
   primoinfection, although a small fraction of patients progress to the
   disease in a smaller amount of time. The introduction of antiretroviral
   therapy (ART) allowed the control of AIDS manifestation, but HIV
   infection remains incurable [[36]1]. On that basis, unraveling new
   mechanisms to control the infection is of the essence to develop new
   therapeutic approaches in chasing a functional cure. In the past
   decade, long-term non-progressors (LTNP) and elite controllers (EC)
   have been extensively studied, due to their ability to maintain high
   CD4+ T cell levels and undetectable viral load in plasma, respectively,
   for a long time in the absence of ART [[37]2]. Due to the heterogeneity
   of these populations, it has been difficult to achieve a consensus
   definition for these extreme phenotypes. Even so, it is widely
   considered that both phenotypes represent excellent natural model
   scenarios of immune control of viral infection.

   Long-term non-progression and viremia control are multifactorial
   circumstances based on viral fitness and immuno-genetic factors. A low
   proportion of these patients are associated with attenuated or
   defective viruses that present mutations in viral genes
   [[38]3,[39]4,[40]5]. Replication-competent HIV-1 has been isolated from
   EC, suggesting that host factors are crucial in the control of the
   infection [[41]6]. Among host factors involved in slow progression, the
   presence of “protective” HLA alleles (HLA-B*57 and HLA-B*27), as well
   as several polymorphisms inside and outside the HLA region have been
   described [[42]7,[43]8,[44]9]. The immune response observed in
   ECs/LTNPs has specific features: higher expression of INF-γ by NK
   [[45]10], polyclonal and intensified Gag-specific response in CD4+ T
   cells [[46]11], and a greater production of IL-2, IFN- γ, MIP-1β, TNF-α
   and granzyme B by CD8+ T cells [[47]12]. Pereyra F, et al. [[48]7]
   outline that the relevant markers of HLA altogether with protective
   variants in CCR2 and CCR5 only account for 23% of the observed variance
   of elite control. Hence, the scientific community strives to find
   protective factors that are present in a relevant percentage of
   ECs/LTNPs.

   Some studies have investigated the gene expression profile of these
   individuals [[49]13,[50]14,[51]15]. They have found differences between
   ECs/LTNPs and typical progressors (TP), which suggest that gene
   expression regulation may be associated with the progression of the
   disease. Different mechanisms including epigenetic and
   post-transcriptional pathways are involved in the regulation of gene
   expression. RNA interference (RNAi) is a conserved mechanism in
   eukaryotic cells composed of various small non-coding RNAs, which
   repress mRNA expression [[52]16]. One of the most important families of
   RNAi is micro-RNAs (miRNAs), and the interest in their role in viral
   infection has increased in the last few years [[53]17]. The mature form
   of miRNAs consists of an asymmetric dsRNA of 19–22 nt which associates
   in the cytoplasm with various proteins to build RNA-induced silencing
   complex (RISC), remaining one strand (guide strand) attached to RISC
   and leading the complex to its target sequences. The degree of
   complementarity of the miRNA with the 3′ untranslated region (3′-UTR)
   of the target gene’s mRNA determines if the targeted mRNA is degraded
   or retained from ribosomal interaction [[54]18]. A single miRNA may
   target multiple gene transcripts and a single mRNA can be targeted by
   various miRNAs. Modulation of HIV replication by miRNAs can be exerted
   in two ways: either directly targeting HIV mRNA or targeting the mRNAs
   that encode relevant host cell factors that are involved in HIV
   replication including genes that contribute to immune responses
   [[55]19]. Some findings show that HIV-1 displays several means to
   hinder host mechanisms based on RNA interference, which suggests a
   potential activity of miRNA in viral suppression.

   The study of miRNA expression profiles has provided novel strategies to
   diagnose, to make prognosis, and even to select a concrete treatment
   for multiple infectious diseases [[56]20,[57]21,[58]22,[59]23]. The
   present study aims to identify a miRNA expression pattern associated
   with ECs/LTNPs phenotypes and through the study of their target genes,
   to get a better insight into the molecular mechanisms implied in
   limiting viral replication and disease progression. The analysis of the
   miRNome of ECs/LTNPs has allowed the identification of deregulated
   miRNAs and the characterization of their potential target routes
   involved in the immune control of the HIV-1 replicative cycle.

2. Experimental Section

2.1. Samples

   Peripheral blood mononuclear cells (PBMC) pellets were kindly provided
   by HIV BioBank integrated into the Spanish HIV/AIDS Research Network
   (RIS) [[60]24]. Peripheral blood mononuclear cells (PBMC) were isolated
   on a Ficoll-Hypaque density gradient (Comercial Rafer S.L., Zaragoza,
   Spain) by the staff of the Spanish HIV BioBank following a Standard
   Operating Procedure (SOP) [[61]25]. Samples were processed following
   current procedures and frozen immediately after their reception. In all
   cases patients’ informed consent was obtained and protocols were
   approved by the Institutional Ethical Committee (Instituto de Salud
   Carlos III. CEI PI 10_2011v3). A total of 30 samples were included in
   the study, 16 HIV-infected subjects classified as long term
   non-progressors (LTNP) from the Spanish LTNP-RIS cohort and 14 paired
   samples coming from 7 HIV-positive patients with a typical progression
   (TP) pattern from the Spanish AIDS Research Network cohort (CoRIS), 7
   before receiving antiretroviral therapy (pre-TP) and their counterparts
   after the administration of ART (post-TP). All experiments were
   performed in accordance with relevant guidelines and regulations.

2.2. Small RNA Extraction

   Pellets consisting of 10^7 PBMCs were received and conserved at −80 °C
   until their usage. RNA extraction was carried out with AllPrep
   DNA/RNA/miRNA isolation kit (Qiagen, Hilden, Germany) following the
   manufacturer’s guidelines. Total RNA quality was assayed in a 2100
   Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA) for quality
   control before DNA library preparation, obtaining a mean RNA integrity
   number (RIN) of 9.2, this being the value above 8 in all samples.

2.3. Next-Generation Sequencing

   The cDNA library was constructed with TruSeq Small RNA library
   preparation (Illumina, San Diego, CA, USA) following the producer’s
   recommendations. Quality control of cDNA library was assayed with High
   Sensitivity DNA kit in Bioanalyzer 2100 (Agilent Technologies, Santa
   Clara, CA, USA), checking for mature miRNAs plus adaptors expected to
   be around 150 bp. Samples were pooled together and sequencing was
   performed on an Hiseq 2000 platform (Illumina, San Diego, USA) with
   single-end 50 nt reads. All samples were loaded onto one lane of a flow
   cell.

2.4. Analysis of miRNA

   Raw data obtained from sequencing underwent quality analysis through
   FASTQC software and java tool Trimmomatic v 0.36 [[62]26] was used to
   trim in every read the adaptors and specific primers generated in cDNA
   library preparation, allowing a 2 bp maximum discordance. Reads were
   quality filtered with a 4 bp sliding window and trimmed if the mean of
   the Phred score of the tetranucleotide fell below 15. Furthermore,
   reads under 20 bp length were removed in order to avoid potential
   contamination of the cDNA library.

2.5. miRNA Alignment

   After sequence cleaning, reads were mapped to known miRNA and
   precursors of miRNAs (pre-miRNA) with miRanalyzer v 0.3 using the
   dataset of human miRNAs contained in miRBase [[63]27]. Filtered files
   were converted to “read-count” format with groupReads.pl script
   built-in miRanalyzer standalone version. This function removes reads
   with less than 5 hits and also trims sequences up to 30 bp. miRanalyzer
   is a java tool that integrates short read alignment to genome and
   miRNome of reference using Bowtie aligner, allowing the introduction of
   external libraries or databases as a reference [[64]28]. This method
   only considers those read counts which just map to one and no other
   sequence of the reference sequence. Annotation was carried out
   sequentially, first, read counts were aligned to known miRNAs and
   removed from the database after matching. Remaining reads that did not
   align with the miRNA database were mapped against pre-miRNA reference
   sequences of miRBase.

2.6. Differential Expression Analysis

   The software package used to analyze differential expression among
   miRNA from distinct groups was DESeq2 [[65]29], which works in R
   environment. This analytic tool is based on a negative binomial
   distribution model and infers probabilities before (p-value) and after
   False Discovery Rate correction (q-value), to spot differentially
   expressed miRNAs among groups of individuals. The variation in
   expression level is provided as log[2] of fold change (FC). Dual
   comparisons were carried out between the aforementioned groups to
   assess the differences in expression levels of previously mapped miRNAs
   and pre-miRNAs. This approximation allows the evaluation of changes in
   miRNA expression due to the viremia level examining typical progressors
   before and after treatment, miRNA expression variations associated with
   LTNP phenotype compared to TP, and the miRNA profile linked to viremia
   control examining EC-LTNP versus vLTNP. Those miRNA with statistically
   significant differences in their expression are defined as
   differentially expressed (DE). R packages pheatmap, factoExtra,
   gridExtra, and ggplot2 were used for PCA and heatmaps analysis.

2.7. Target Genes Identification and Functional Analysis

   Target genes of DE miRNAs were predicted and those that might be
   involved in HIV infection immunopathogenesis were further analyzed.
   First, miRNAs were linked to targets related to HIV infection through
   the information available in scientific literature, and subsequently,
   the predicted union sites and type of union in 3′-UTR were checked in
   TargetScanHuman v 7.1 [[66]30].

   Likely modifications in signaling pathways were assessed via miRNet,
   which allows making miRNA target identification and refinement, network
   visual exploration, and enrichment of altered pathways. We selected the
   DE miRNAs with the most robust statistical significance (LTNP vs. TP
   comparison), and we studied which genes were targeted specifically in
   two databases. The first one is composed of 212 genes included within
   Kyoto Enciclopedia of Genes and Genomes (KEGG) Human Immunodeficiency
   Virus 1 Interaction—Homo sapiens (human) pathway (hsa05170) [[67]31].
   Once the interacting genes with DE miRNA in this list were identified,
   we applied a filter within the HIV-1 infection KEGG route, highlighting
   those pathways which could be altered in LTNP. Alternatively, we used
   the 4668 genes retrieved from the HIV-1 host proteins interactome
   database available at NCBI [[68]32,[69]33,[70]34,[71]35,[72]36] to
   filter the interacting genes with the DE miRNAs identified in LTNPs
   when compared to TPs. Finally, we conducted a pathway enrichment
   analysis within both networks using an empirical sampling approximation
   [[73]37].

3. Results

3.1. Patient Characteristics

   Patient’s records were provided by the LTNP-RIS and CoRIS cohorts
   alongside the classification of the samples—all the information is
   available in [74]Table 1. Positive-HIV patients are considered LTNP
   when their CD4+ T lymphocyte counts are above 500 cells/mm^3 for at
   least 10 years, in the absence of ART. This group is divided into two
   subgroups (n = 8) according to their capacity to control viral
   replication: elite controllers LTNP (EC-LTNP), defined as those with
   undetectable viremia or viral load < 2000 copies of viral RNA/mL in
   less than 25% of determinations during the follow-up, and viremic LTNP
   (vLTNP) that includes the other individuals with a viral load below
   10,000 copies/mL in all the determinations during the follow-up. The
   whole set of LTNP was followed up to 10 years and diagnoses were
   registered between 1988 and 1999. Typical progression was defined as a
   loss of 50–100 CD4+ T cells/mm^3 per year with a presence of >5000
   copies of viral RNA/mL. After two years of ART, viral load declined
   below the detection level (<20 copies/mL). There are no noteworthy
   differences in age, gender, or provenance between the established
   groups.

Table 1.

   Clinical parameters of the patients included in the analysis.
   Group Patient ID Sex Origin CD4+ T Cell Count
   (Cell/mm^3) Via of Infection HCV Co-Infection
   EC-LTNP EC-LTNP-1 M European 635 IDU Yes
   EC-LTNP-2 M European 555 IDU Yes
   EC-LTNP-3 F European 596 IDU Yes
   EC-LTNP-4 F European 1081 SEX No
   EC-LTNP-5 F European 603 SEX Yes
   EC-LTNP-6 M European 640 IDU Yes
   EC-LTNP-7 M European 514 IDU Yes
   EC-LTNP-8 M European 708 IDU Yes
   vLTNP vLTNP-1 M European 625 SEX No
   vLTNP-2 M European 590 IDU Yes
   vLTNP-3 M European 1049 IDU Yes
   vLTNP-4 M European 593 IDU Yes
   vLTNP-5 F European 492 IDU Yes
   vLTNP-6 M European NA SEX No
   vLTNP-7 M European 889 IDU No
   vLTNP-8 M European 929 IDU Yes
   pre-TP pre-TP-1 M European 290 IDU Yes
   pre-TP-2 M European 332 SEX No
   pre-TP-3 M European 397 SEX No
   pre-TP-4 F European 270 SEX No
   pre-TP-5 M European 624 SEX No
   pre-TP-6 M European 157 SEX No
   pre-TP-7 M European 47 SEX No
   post-TP post-TP-1 M European 720 IDU Yes
   post-TP-2 M European 525 SEX No
   post-TP-3 M European 550 SEX No
   post-TP-4 F European 650 SEX No
   post-TP-5 M European 439 SEX No
   post-TP-6 M European 415 SEX No
   post-TP-7 M European 578 SEX No
   [75]Open in a new tab

   M: Male; F: Female; UDI: Intravenous Drug User; SEX: Sexual
   transmission.

3.2. Next-Generation Sequencing Quality Assessment

   A total amount of 393,350,988 reads were obtained in sequencing, with
   an average of 13,211,670 reads per patient. After adaptor trimming and
   data filtration as referred to in methods, the number of reads
   available for mapping reached 73.44% of initial reads (9,703,085 reads
   per patient). All sequences had an average Phred score per base above
   38 after filtering. There were no statistically significant differences
   between the remaining numbers of reads between the phenotypical groups.

3.3. miRNA Mapping

   Filtered sequences were mapped against different databases to obtain
   annotation of the reads. An average of 4,499,757 reads was mapped to
   known miRNAs and 4,692,188 reads were annotated as pre-miRNAs per
   patient. PBMC miRNome is formed by a mean (±SD) of 355 (±39) miRNAs and
   472 (±48) pre-miRNAs ([76]Figure 1). The hsa-miR-26a was identified as
   an overrepresented sequence in the whole set of patients, accounting up
   to a mean (SD) of 12.57% (±0.035%) of the aligned reads to miRNAs.

Figure 1.

   [77]Figure 1
   [78]Open in a new tab

   Number of identified mature miRNAs and precursor miRNAs (pre-miRNA) per
   patient.

3.4. Global Comparison

   In order to describe miRNA expression alteration, study groups were
   compared in pairs. The analysis of differential expression confirms the
   modification of miRNA expression patterns related to distinct
   phenotypes. [79]Figure 2 shows the lists of DE miRNAs generated in each
   comparison.

Figure 2.

   [80]Figure 2
   [81]Open in a new tab

   Summary of differentially expressed miRNAs (q-value < 0.1). Each
   comparison is represented as a circle split in two colors, green for
   upregulated miRNAs and red for downregulated ones. Overlapping miRNAs
   are depicted inside grey boxes. Solid lines illustrate upregulation in
   the corresponding group, while dotted lines show downregulation.
   Underlined phenotype indicates the reference group in each comparison.
   *Statistical significance of DE miRNA p-value < 0.05. LTNP (Long-term
   non-progressor); TP (Typical progressor); EC-LTNP (Elite controller
   Long-term non-progressor).

   Principal component analysis (PCA) was carried out so that the main
   miRNAs associated with each phenotype could be outlined (see [82]Figure
   3). Each graph was generated by selecting the DE miRNAs found in each
   comparison, in order to avoid background noise. At the representation
   of the whole group of samples, LTNPs cluster all together ([83]Figure
   3A). The most remarkable finding is that 80% of LTNPs present a
   negative value for PC1 and only hsa-miR-3182 contributes with a
   negative loading value. Coefficients from hsa-miR-18a-5p, hsa-3613-5p,
   and hsa-miR-324-5p contribute positively on this component. At PC2,
   hsa-miR30d-5p, hsa-miR-3182, and hsa-miR324-5p play a role as negative
   coefficients while hsa-miR451a, hsa-miR-144-3p, and hsa-miR-144-5p are
   the most relevant positive coefficients. Performing PCA in separate
   groups reveals discrimination of phenotypes according to ART
   administration, [84]Figure 3B, and viral load control in LTNP
   individuals ([85]Figure 3C).

Figure 3.

   [86]Figure 3
   [87]Open in a new tab

   PCA analyses were carried out to identify relationships between
   distinctive miRNA expression patterns and rate of progression or viral
   load. DE miRNAs determined in each comparison were selected to generate
   the graphs (A) LTNP vs. TP, (B) pre-TP vs. post-TP, and (C) vLTNP vs.
   EC-LTNP. [88]Supplementary Figures S1–S3 show the contribution and
   loadings of each miRNA in PCs. (A) All LTNP patients group together in
   negative values of PC1 and PC2. TP patients do not group attending to
   any PC values. (B) TP individuals are only separated by PC1. Patients
   receiving ART (blue) show lower values whereas its value varies before
   treatment (red). (C) EC-LTNPs present high scores on PC1 except for two
   patients which are close to vLNTPs scores.

3.5. Differences in miRNA Expression Associated with Viral Replication

   Attending to paired comparisons individually, we first assessed TP
   groups in order to register miRNA expression profile alterations due to
   the presence of viremia. By comparing typical progressors before ART
   administration (pre-TP) and the same patients after ART (post-TP), 18
   miRNA (q-value < 0.1), and 11 pre-miRNA (p-value < 0.05) were found to
   be dysregulated. Among them, hsa-miR-3614-5p (FC: −1.726),
   hsa-miR-3614-3p (FC: −2.618), hsa-miR-18a-5p (FC: −0.782), and
   hsa-miR-18b-5p (FC: −0.729) are downregulated after treatment, and
   miRNAs hsa-miR-3607-5p (FC: 2.050) and hsa-miR-31-5p (FC: 1.597)
   display upregulation. Although statistical significance after FDR is
   not reached for pre-miRNAs, hsa-pre-miR-1303 and hsa-pre-miR-4454
   suffer a reduction in post-ART patients, with FC of (−1.481) and
   (−2.360) respectively. A complete list of DE miRNAs and pre-miRNAs can
   be consulted in [89]Supplementary Table S1. In [90]Figure 4A, we
   observe that pre-TP patients are highly differentiated from post-TP
   patients based on the expression level of the whole set of DE miRNA and
   pre-miRNA. Some of those miRNAs are present in other comparisons
   ([91]Figure 2).

Figure 4.

   [92]Figure 4
   [93]Open in a new tab

   Heatmaps of differentially expressed miRNA and pre-miRNA for each
   comparison. Expression levels are represented as log[2] of normalized
   counts per patient. Individuals were clustered using an unsupervised
   Euclidean approximation based on miRNA and pre-miRNA counts. The color
   scale on the right indicates the relative expression level of each
   miRNA, green being the higher values and red being the lower scores.
   (A) pre-TP vs. post-TP. Patients are split into their corresponding
   group accurately, showing that this set of miRNA and pre-miRNA is
   associated with the administration of treatment. (B) LTNP vs. TP. LTNP
   individuals are differentiated from TP, finding five mismatches. This
   panel of expression allocates correctly 71% of TP and 68.5% of LTNP.
   (C) EC-LTNP vs. v-LTNP. vLTNPs are associated together except for one
   patient who clusters with EC-LTNP.

   [94]Table 2 presents a full list of experimentally validated targeted
   genes implicated in the HIV replication cycle, predicted union site,
   and type of union based on TargetScan v 7.0 analysis.

Table 2.

   List of genes involved in AIDS immunopathogenesis targeted by DE
   miRNAs.
   Group Gene miRNA Expression Change 3′ UTR Predicted Site Canonical Site
   Types ^1 References ^2