RNA sequencing analysis of activated macrophages treated with the anti-HIV ABX464 in intestinal inflammation

RNA-Seq enables the generation of extensive transcriptome information providing the capability to characterize transcripts (including alternative isoforms and polymorphism), to quantify expression and to identify differential regulation in a single experiment. To reveal the capacity of new anti-HIV ABX464 candidate in modulating the expression of genes, datasets were generated and validated using RNA-seq approach. This comprehensive dataset will be useful to deepen the comprehensive understanding of the progression of human immunodeficiency virus (HIV) associated with mucosal damage in the gastrointestinal (GI) tract and subsequent inflammation, providing an opportunity to generate new therapies, diagnoses, and preventive strategies.


Background & Summary
RNA sequencing (RNA-seq) transcriptome analysis using massively parallel next generation sequencing technology provides the capability to understand global changes in gene expression throughout a range of tissue samples and have become the gold standard for in depth characterization of novel transcript isoforms. RNA-seq (Quantification) is used to analyse gene expression of certain biological objects under specific conditions. It is a cost-effective quantification method that produces high reproducibility, high accuracy and wide dynamic range. As described in our related manuscript submitted earlier in Scientific Reports 1 , our aim in this study is to consolidate the protective effect of ABX464 in DSS-induced colitis comparing transcript levels between bone marrow-derived macrophages (BMDMs) ABX464-treated stimulated or not by LPS. To perform this task we sequenced here 24 samples of C57BL6 strains of mice using BGI RNA sequencing (RNA-seq) technology 2,3 .
ABX464 is a first-in-class antiviral drug candidate used to treat patients infected with human immunodeficiency virus. ABX464 is an orally available, well-tolerated small molecule that blocks HIV replication through the modulation of RNA biogenesis. ABX464 demonstrates anti-viral activity in treatment-naïve patients and induces a long-lasting control of the viral load in HIV-infected humanized mice after treatment arrest. This treatment might improve current therapies by hampering the appearance of virus resistance and, more importantly, lead to a long-lasting control of HIV.
Despite the successful control of viremia, many HIV-infected individuals treated with antiretroviral therapy (ART) exhibit residual inflammation associated with non-AIDS-related morbidity and mortality. Several reports have shown that measures of inflammation and immune activation are the best independent predictors of disease progression in HIV-infected individuals. It is also known that elevated inflammation and immune activation predict non-AIDS-associated morbidities and mortality, even in well-suppressed treated HIV-infected individuals. Mucosal damage to the gastrointestinal (GI) tract with resulting microbial translocation significantly contributes to the heightened and persistent chronic inflammation and immune activation characteristic to HIV infection and might contribute to virus persistence during ART. Recently, a dextran sulfate sodium (DSS)-induced nonhuman primate (NHP) colitis model in SIV-infected rhesus macaques a natural host species for SIV that does not manifest GI tract damage or chronic immune activation during infection, was used to demonstrate that DSS-exposure results in colitis with elevated levels of plasma SIV RNA, sCD14, LPS, CRP and mucosal CD4 T cell loss 4,5 . The pathogenic potential of bacterial translocation has also been studied in humanized mice exposed to DSS to provide evidence that HIV chronic infection contributes to systemic inflammation under DSSexposure 6,7 . Given the similarities in colitis induction between DSS-exposed macaques and mice 5 , the antiviral activity of ABX464 in humanized mice model 8 and its pharmacodynamics and pharmacoketics profiles in rodent, NHP 8 and humans 9 , the anti-inflammatory properties of this drug in DSS-treated mice were tested (Chebli et al.). The results presented in this manuscript based on the use of DSS-induced colitis in mice, show for the first time that an anti-HIV drug can specifically act on the immune system to resolve the inflammation.
The principal aim of this study was to provide a comprehensive resource to deposit the raw RNA-seq data sets underpinning the advanced high-resolution transcriptome analysis of ABX464-treated BMDMs stimulated with LPS. To this end the data from three distinct group of BMDM cultures were analyzed for different gene expression levels using RNA-seq. We performed a high-resolution transcriptome analysis of ABX464-treated BMDMs stimulated or unstimulated with LPS. In this work, the RNA-seq libraries were sequenced using the BGISEQ-500 platform which use the DNA nanoball technology. Approximately 1 billion 50-bp single-end RNA-seq reads were generated, with an average of over 45 million sequence reads per sample.
24 intestinal macrophages mRNA profiles of three distinct group of BMDM cultures from mice were generated by deep sequencing, using Illumina protocol. The sequence reads that passed quality filters were analyzed at the transcript isoform level with two methods: HISAT 10 and Bowtie2 (ref. 11) followed by genes expression quantification using a software package called RSEM 12 .
Using an optimized data analysis workflow, about 45 million sequence reads per sample were mapped to the mouse genome (build mm9) and identified 1,383 expressed genes in the LPS-stimulated BMDMs among which 439 and 397 genes increased and decreased in expression > = 1.5 log2-fold, respectively. Altered expression of some genes was confirmed with qPCR, demonstrating the high degree of sensitivity of the RNA-seq method. Hierarchical clustering of differentially expressed genes uncovered several as yet uncharacterized genes that may contribute to inflammatory function. Data analysis with HISAT 10 / Bowtie2 (ref. 11) and RSEM 12 workflows revealed a significant overlap yet provided complementary insights in transcriptome profiling. A comprehensive RNA-seq data resource was provided on the antiinflammatory properties of ABX464 that warrant exploration in both HIV and inflammatory ulcerative colitis (UC) disease. The data generated here will be valuable to the vision research community for characterizing global changes in gene expression and facilitate a better understanding of the complexities during anti-AIDS drugs treatment and then can help to develop better ways to combat HIV.

Macrophages RNA-seq preparation
Cell isolates from the femurs of C57BL6 mice (in Table 1  stimulation time points and treatments) were cultured for 9 days in the presence of GM-CSF (50 ng ml −1 ; Peprotech). Cell cultures were treated for the last 3 days with either ABX464 (5 μM) suspended in DMSO or DMSO alone. After additional stimulation for 6 h with LPS (4 μg ml −1 ; Invivogen) or unstimulation (NS), the supernatants were analyzed for cytokine production using a cytometric bead array technology (CBA; Becton Dickinson). The detection limits were 5 pg ml −1 for IL-6, 52.7 pg ml −1 for MCP-1, 17.5 pg ml −1 for IL-10, 7.3 pg ml −1 for TNF α, and 10.7 pg ml −1 for IL12p70 (Results presented in ref. 1). The cells were lysed in the Qiazol and the RNAs were purified with the miRNeasy kit from Qiagen. Dosage and analysis are performed using Nanodrop/Qubit and 2,100 Bioanalyzer (Agilent Technologies).

RNA-seq library preparation and sequencing
DNase I was used to digestion step, DNase I degrade double-stranded and single-stranded DNA contaminant in RNA samples. RNA molecules were purified using oligo (dT)-attached magnetic beads. The RNA intensity was checked using a 2,100 Bioanalyzer (Agilent Technologies) and only high quality samples with an RNA Integrity Number (RIN) value greater than or equal to 7 were used to construct the sequencing library. Then mRNA were fragmented into small pieces using fragmentation reagent. In cDNA synthesis first-strand cDNA was generated using random hexamer-primed reverse transcription, then was followed by a second-strand cDNA synthesis. Typically, 1 μg of total RNA was used with the TruSeq RNA library preparation kit (Illumina) in accordance with low-throughput protocol, except that SuperScript III reverse transcriptase (Invitrogen) was used to synthesize first strand cDNA. In the step of end repair the synthesized cDNA was subjected to end-repair and then was 3′ adenylated. Adaptors were ligated to the ends of these 3′ adenylated cDNA fragments. We used PCR to amplify the cDNA fragments with adaptors from previous step. PCR products are purified with the SPRI beads, and dissolved in EB solution. In the circularization step the double stranded PCR products were heat denatured and circularized by the splint oligo sequence. The single strand circle DNA (ssCir DNA) were formatted as the final library. After PCR enrichment and purification of adapter-ligated fragments, the concentration of DNA with adapters was determined with quantitative PCR (Applied Biosystems 7,500) using primers QP1 5′-AATGATACGGCGACCACCGA-3′ and QP2 5′-CAAGCAGAAGACGGCATACGAGA-3′. The length of the DNA fragment was measured using a 2,100 Bioanalyzer. The library was amplified with phi29 to make DNA nanoball (DNB) which have more than 300 copies of one molecular. The DNBs were load into the patterned nanoarray and single end 50 bases reads were generated in the way of sequenced by synthesis. Then, RNA-seq libraries were sequenced using the BGISEQ-500 sequencer. Base-calling was performed by the BGISEQ-500 software version 0.3.8.1111.
Data filtering 'Dirty' raw reads were defined as reads which contain the sequence of adaptor, high content of unknown bases and low quality reads. These reads were removed before downstream analysis to decrease data noise. Filtering steps are as follows: The quality of the RNA-seq libraries was first evaluated using the FastQC v0.11.5 software. Then, the reads were subjected to standard quality control (QC) criteria according to the following parameters: (1) trimming and cleaning reads that aligned to primers and/or adaptors, (2) reads with over 50% of low-quality bases (quality value ≤5) in one read, and (3) reads with over 10% unknown bases (N bases). After filtering, the remaining reads are called 'clean reads' and stored as FASTQ format 13 .

RNA-seq data analysis
24 samples of C57BL6 strains of mice were sequenced using RNA-Seq technology 2,3 , averagely generating 45,188,386 raw sequencing reads and then 45,120,825 clean reads after filtering low quality (see Data Filtering in method page). Table 1 briefly summarizes information of sequencing data for each sample.
After filtering, clean reads are mapped to reference using HISAT 10 /Bowtie2 (ref. 11) tools. The average mapping ratio with reference gene is 83.67% and Table 2 lists separate mapping rate for each sample. The average genome mapping ratio is 91.47%. Strict quality control was conducted for each sample to evaluate whether the sequencing data are qualified.
Sequencing data saturation analysis is used to measure whether the depth of sequencing data is sufficient for informatics analysis. With the number of sequenced reads increasing, the number of identified genes is also increased. However, when the number of sequenced reads reaches a certain amount, the growth curve of identified genes flattens, indicating that the number of identified genes tends to reach the saturation. The distribution of the reads on reference gene reflects whether each part of the gene body is evenly sequenced. If the randomness is good, the reads in every position (from 5' terminal to 3' terminal) would be evenly distributed. If the randomness is poor, reads preference to specific gene region will directly affect subsequent bioinformatics analysis. We use the distribution of reads on the reference genes to evaluate the fragmentation randomness, as shown in Fig. 1. Gene's expression level is quantified by a software package called RSEM 12 . The number of identified expressed genes were counted and its proportion to total gene number in database for each sample was calculated. Meanwhile, the distribution of gene number on different expression level for each sample are shown as Fig. 2, to get a    DEGs (Differentially Expressed Genes) screening is aimed to find differential expressed genes 14 between samples and perform further function analysis on them. Poisson distribution method was used to perform this analysis.
All expressed genes of each pairwise are stored in *.GeneDiffExp.xls and screened DEGs are stored in *.GeneDiffExpFilter.xls. They have the same file format and are listed in Figshare repository (Data Citation 1).
For result list of each control-treatment pair above (*.GeneDiffExp.xls), scatter plots were draw for all expressed genes (Fig. 3) to present the distribution of DEGs in screening threshold dimensions. Also, a histogram is shown to demonstrate the number of genes that are significantly up-down in each pairwise comparison (Fig. 4).  Genes with similar expression patterns usually have same functional correlation. We also provide clustering analysis of differentially expressed genes with cluster 15,16 and javaTreeview software 17 according to the provided cluster plans for DEGs. Complete intersection and union DEGs heatmap for each cluster plan were provided (Fig. 5).
Annotation analysis of Gene Ontology (GO) are performed for screened DEGs and then generate a report named GOView.html under the data folder of Figshare (Data Citation 1).
After getting GO annotation for DEGs, WEGO software 18 was used to do GO functional classification to help understand the distribution of gene functions of the specie from the macro level.
Genes usually interact with each other to play roles in certain biological functions. We perform pathway enrichment analysis of DEGs based on KEGG database 19 and generate a report for DEGs in each pairwise respectively, stored in Figshare folder (Data Citation 1). In addition, we generate a scatter plot for the top 20 of KEGG enrichment results (Fig. 6) and a bar plot for the statistics of KEGG terms types (Fig. 7).
Different proteins often form complex proteins through complicated interactions to perform their biological functions. Protein-protein interaction network analysis integrates several famous interaction network databases such as BIND, BioGrid and HPRD and constructs interaction networks for DEGs' coding proteins. The analysis result documents are showed by Medusa software 20 and webpage version is presented by the report network_en.html under Figshare folder (Data Citation 1).

Statistical data analysis
Referring to 'The significance of digital gene expression profiles' 14 , we have developed a strict algorithm to identify differentially expressed genes between two samples. Denote the number of unambiguous clean tags (which means reads in RNA-seq) from gene A as x, given every gene's expression occupies only a We do correction on P-value to differential gene expression test using Bonferroni method 21 . Since DEG analysis generate a large multiplicity problems in which thousands of hypothesis (is gene x differentially expressed between the two groups) are tested simultaneously, correction for false positive (type I errors) and false negative (type II) errors are performed using FDR method 22 . Assume that we have picked out R differentially expressed genes in which S genes really show differential expression and the other V genes are false positive. If we decide that the error ratio 'Q = V/R' must stay below a cutoff (e.g., 5%), the FDR was set to a number no larger than 0.05. We use 'FDR ≤0.001 and the absolute value of Log2Ratio ≥1' as the default threshold to judge the significance of gene expression difference. More stringent criteria with smaller FDR and bigger fold-change value can be used to identify DEGs.

Technical Validation
Quality control-RNA integrity To determine the quality of the RNA-seq libraries, RNA intensity was checked using a 2,100 Bioanalyzer (Agilent Technologies). The integrity of the total RNA was calculated by the RNA Integrity Number (RIN) algorithm, which can be used to estimate RNA quality and degradation level. Higher RIN values indicate better integrity of total RNA, with the highest RIN value being 10. In this study, only total RNAs with a RIN value of greater than or equal to 7 were used to construct the sequencing library, which shows RNA with high integrity was used.

Quality validation and analyses
We applied FastQC v0.11.5 software to determine data quality and analysed several measures, including the sequence quality per base, GC content per sequence, sequence duplication levels, and quality score distribution over all sequences of the fastq files 13 . A representative summary plots are depicted in Figshare repository (Data Citation 3). Here, the quality scores per base were high, with a median quality score above 30 suggesting high quality sequences across all bases. The GC distribution per base over all sequences was examined. The GC composition pattern was similar to theoretical distribution indicating that the samples were contaminant free. In addition, sequence complexity was examined. Approximately 15% of sequences were shown over 10 to 100 times, either from highly expressed transcripts duplications. The quality score distribution over all sequences was analysed to see if a subset of sequences had universally poor quality. The average quality for most sequences was high, with scores above 37, which indicated that a significant proportion of the sequences in a run had overall high. All other fastqc files were shown to have similar quality metrics.