Signature of circular RNAs in peripheral blood mononuclear cells from patients with active tuberculosis

Abstract The study was to characterize the expression profiles of circular RNAs (circRNAs) in peripheral blood mononuclear cells (PBMCs) from active tuberculosis (TB) patients and to investigate their function. Microarray was applied to detect circRNA expression and reverse transcription‐quantitative polymerase chain reaction was conducted to validate the microarray results. Meanwhile, receiver operating characteristic curve (ROC) curve was calculated to evaluate the predictive power of the selected circRNAs for TB diagnosis. Additionally, circRNA/miRNA interaction was predicted based on miRNA target prediction software, and gene ontology as well as Kyoto Encyclopedia of Genes and Genomes pathway analysis were used to predict their biological function. In total, 171 circRNAs were found to be dysregulated in TB samples. Specifically, circRNA_103017, circRNA_059914 and circRNA_101128 were confirmed to be increased, while circRNA_062400 was decreased in TB samples. ROC analysis revealed that circRNA_103017 had potential value for TB diagnosis, followed by circRNA_059914 and circRNA_101128. Moreover, circRNA_101128 expression in TB samples was negatively correlated with the level of its possible target let‐7a and bioinformatics analysis showed that circRNA_101128 was potentially involved in MAPK and P13K‐Akt pathway possibly via modulation of let‐7a. Taken together, our results indicated that some dysregulated circRNAs were potential biomarkers for the diagnosis of TB and circRNA_101128‐let‐7a interplay may play considerable role in PBMCs response to Mtb infection.


| INTRODUCTION
Tuberculosis (TB) remains a common deadliest and communicable disease caused by Mycobacterium tuberculosis (Mtb), especially in low-and middle-income countries. 1 Early identification and proper therapy can obviously improve active TB patients' outcome. 2 Great advances and many efforts have been made to treat and understand basal mechanisms of TB infection in recent years. However, remission rates remain suboptimal and TB diagnosis is still challenging due to the limitations in the specificity and sensitivity of the current diagnostic tests. 3 It may be because the underlying molecular mechanism of TB remains largely elusive.
Genetic alterations are the key to revealing the pathogenesis of active TB. Circular RNAs (circRNAs), a recently discovered class of endogenous non-coding RNAs, have been shown to be abundantly expressed in human cells and characterized by covalently closed continuous loop without free 3′ or 5′ ends, which confers them high  4 Therefore, this makes circRNAs to be more suitable diagnostic biomarkers than other types of RNA. Meanwhile, emerging evidence has suggested that circRNAs play crucial roles in the regulation of gene expression at both transcriptional and post-transcriptional levels, and their tissue-specific as well as disease-associated expression patterns have been proposed as significant regulators for underlying pathological mechanisms and therapeutic targets for diseases. 5 However, to our knowledge, their expression patterns and functions in active TB are still limited.
This study was designed to determine whether circRNAs in peripheral blood mononuclear cells (PBMCs) could be used as novel biomarkers for active TB diagnosis and whether aberrantly expressed circRNAs may provide new suggestion for further understanding of the potential mechanisms of TB pathogenesis and therapeutic strategy. Age-and gender-matched healthy patients were from the Affiliated Hospital of Weifang Medical University and served as controls in the study. In addition, healthy controls with recent close contact of TB patient or history of TB were excluded from the study. Meanwhile, we carried out interferon-γ release assay to detect whether the controls were not latently infected with Mtb. The basic demographic characteristics of all the participants in the study were summarized in Table S1.

| Participants and PBMC samples preparation
Fasting venous blood samples were obtained at diagnosis prior to treatment initiation. After sample collection, PBMCs were freshly isolated by Ficoll gradient centrifugation. Resulting PBMC samples were aliquoted and immediately kept in liquid nitrogen until further analysis. Three PBMC specimens were choosed randomly from each group (male/female = 2/1) and were used for microarray analysis. All of the PBMC samples were applied for subsequent reverse transcription-quantitative polymerase chain reaction (RT-qPCR) verification.
All study protocols were approved by Weifang Medical University local ethics committee and were carried out in compliance with the Code of Ethics of the World Medical Association (Declaration of Helsinki). Written informed consent was obtained from all the participants prior to their enrolment in the study. For data analysis, Agilent Feature Extraction software (version 11.0.1.1) was applied to analyse the acquired slide images. Then, R software package (R version 3.1.2) was used for quantile normalization of raw data and subsequent data processing. Next, low intensity filtering was carried out and the circRNAs that at least three out of six samples have flags in "P" or "M" ("All Targets Value") were left for further analysis. CircRNAs with both false discovery rate (FDR) < 0.05 and fold-change (FC) value >2 were considered to be differentially expressed. Moreover, scatter plot was performed to show the variation in circRNA expression and volcano plot filtering was used to identify differentially expressed circRNAs with statistical significance. Meanwhile, hierarchical clustering analysis was applied to show the dysregulated circRNAs expression profile among the samples.

| Reverse transcription-quantitative polymerase chain reaction
Total RNA was prepared from each PBMC sample using TRIzol ® reagent (Invitrogen), and then, PrimeSuperScript TM Ⅱ(TaKaRa,)was utilized to synthesize cDNA templates from 200 ng of total RNA. qPCR was performed in a volume of 20 μL containing 2 μL of cDNA and 0.5 μmol/L of each forward and reverse primers (primer sequences were available if required). PCR thermocycling conditions were set at 95°C for 10 minutes for pre-denaturation, followed by 40 cycles at 95°C for 10 seconds and 60°C for 60 seconds. Housekeeping gene GAPDH was applied as an internal reference. All reactions were performed in triplicate and the 2 ÀΔΔCt method was used to calculate relative gene expression levels.

| Receiver operating characteristic analysis
Receiver operating characteristic curve (ROC) was constructed and the area under the curve (AUC) value as well as 95% confidence intervals (CI) were calculated to evaluate the predictive power of the selected circRNAs for active TB diagnosis.

| Annotation for circRNA/miRNA interaction and KEGG pathway analysis
circRNA/miRNA interaction was predicted with Arraystar's homemade miRNA target prediction software TargetScan (http://www.ta rgetscan.org/) and miRanda (http://www.microrna.org/) to annotate miRNA response elements (MREs) potentially targeted by the confirmed dysregulated circRNAs based on their seed-sequence complementarity. Function of hsa_circRNA_101128 was predicted based on the analysis of its targeted miRNA-mRNAs by gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis.

| Statistical analysis
Student's t test and chi-squared test were performed to evaluate the statistical significance between the two groups. Correlation between circRNA and its possible targeted miRNA among individual samples was assessed. A P value of <0.05 was viewed as statistically significant.

| Analysis of circRNA expression profiles in PBMCs
To identify the potential circRNAs related to human PBMCs response to Mtb infection, we characterized the circRNAs expression profiles using microarray assays. As present in FigureS1A and B, scatter plot visualization showed the normalized expression level of circRNAs and volcano plot displayed the differential expression of circRNAs. When the screening criteria were set as FC ＞ 2 and P < 0.05, a total of 171 circRNAs (FC ＞ 2 and P value <0.05 adjusted with FDR) were found to be differentially expressed; of these, 145 were up-regulated while 26 were down-regulated in the active TB group compared with the healthy controls. Particularly, we found that hsa_circRNA_059914 was the highest while hsa_-circRNA_062400 showed the largest down-regulation in the TB group, respectively. Furthermore, when the screening criteria were set as FC ＞ 3 and P < 0.05, 55 circRNAs were identified to be significantly dysregulated (Table S2), of which, 51 were increased while four were decreased in TB samples. At the same time, hierarchical cluster analysis was used to visualize the differential expression patterns of the 55 circRNAs among the samples in the present study ( Figure 1). Our results demonstrate that circRNAs alterations are involved in human PBMCs response to Mtb infection.
Microarray expression data set in this study has been submitted to the Gene Expression Omnibus (GEO, http://www.ncbi.nlm. nih.gov/geo) database and the GEO accession number is GSE117563.

| Confirmation of circRNA expression
To further verify the microarray results, the top three increased (hsa_-circRNA_059914, hsa_circRNA_103017 and hsa_circRNA_101128) and the top one decreased circRNAs (hsa_circRNA_062400) were selected for confirmation by RT-qPCR in a total of 61 samples from 31 active TB patients and 30 healthy individuals. As shown in Figure 2A, hsa_-circRNA_059914, hsa_circRNA_103017 and hsa_circRNA_101128 were significantly higher, while hsa_circRNA_062400 was obviously lower in TB samples as compared to healthy controls (P < 0.05). Our findings confirmed the data of the microarray analysis.

| Diagnostic value of the confirmed circRNAs in TB infection
Recently, many studies have suggested that circRNAs can serve as novel diagnostic markers for diseases. Therefore, ROC curve was calculated to assess whether the four confirmed circRNAs serve as diagnostic biomarkers for TB infection. As shown in Figure 3

| Annotation for circRNA/miRNA interaction
Recent studies have demonstrated that circRNAs play an important role in acting as an endogenous "sponge" to sequester miRNA-mediated regulation of mRNA expression, 6,7 which suggests that cir-cRNAs have crucial functions in the regulation of miRNA expression.
To investigate circRNAs possible functions, we analysed potential miRNAs binding with circRNAs using Arraystar's home-made miRNA target prediction software. In this study, we focussed on four confirmed circRNAs (hsa_circRNA_103017, hsa_circRNA_059914, hsa_-circRNA_101128 and hsa_circRNA_062400) to predict their MREs.
Five MREs with good mirSVR scores for each circRNA were present in Table 1; of these miRNAs, it is worth noting that certain miRNAs have been reported to be dysregulated in TB, such as let-7a and let-7f, 8-10 which were decreased upon Mtb infection. More importantly, these above-mentioned miRNAs were possible targets of hsa_-circRNA_101128, which was up-regulated in this study. The detailed potential interplay between hsa_circRNA_101128 and its target miR-NAs including let-7a and let-7f was shown based on sequence-pairing miRNA target prediction (Figure 4). To further investigate whether there was a relationship between hsa_circRNA_101128 and let-7a as well as let-7f in PBMCs, expression of both miRNAs was FU ET AL. Confirmation of the expression of candidate differentially expressed circRNAs. hsa_circRNA_103017, hsa_circRNA_059914 and hsa_circRNA_101128 were enhanced, while hsa_circRNA_062400 was reduced in the active TB patient group (n = 31) vs healthy patient group (n = 30). (B) Prediction of possible miRNA targets of hsa_circRNA_101128. hsa-let-7a and hsalet-7f were decreased in the case group (n = 31) compared with the controls (n = 30). (C) Expression level of hsa_circRNA_101128 was negatively correlated with that of let-7a (r = −0.56, P ＜ 0.05). GAPDH and U6 were applied as reference genes of circRNAs and miRNAs, respectively. Data were analysed with the 2 ÀΔΔCt method. Student's t test was utilized for statistical analysis.  (Table S3).

| DISCUSSION
As TB is one of the leading causes of death caused by a single pathogen worldwide, it is imperative to identify novel diagnostic biomarkers and therapeutic targets to improve patient's survival. Recent evidence has indicated that circRNA plays important roles in various pathological processes and its dysregulation has also been involved in the occurrence and progression of many diseases, such as viral infection. 11,12 Due to its important function, relationship between circRNA and TB infection has started to be gained attention in the past 2 years and it was found that circRNAs dysregulation was involved in Mtb infection. For example, hsa_circ_0005836 was reported to act as a new possible biomarker for active TB disease 13 ; circRNA_00074, circRNA_09585, circRNA_14623, circRNA_005538, circRNA_09993 and circRNA_13478 were found to be differentially expressed in leukocytes from active TB patients. 14 These findings imply that cir-cRNAs may play important roles in active TB and have the great potential to serve as a novel type of biomarkers for TB diagnosis.
Despite some efforts devoted to uncovering circRNAs profiling during the last 2 years, in reality, the knowledge of circRNAs in active TB still remains in its infancy stage.
In this study, we conducted a comprehensive microarray analysis of global changes in the expression pattern of circRNAs in PBMCs from active TB patients. We found that 171 circRNAs were dysregulated in TB samples, which implicated that these aberrantly expressed circRNAs were potentially associated with TB infection.
Of them, 145 circRNAs were elevated while only 26 circRNAs were reduced in TB samples compared with the controls, which were in line with the report that the circRNA's magnitude in PBMCs from active TB patients was significantly increased. 15 Moreover, our results showed that the expression pattern of circRNAs in PBMCs is similar to that in plasma, 16  F I G U R E 4 Detailed presentation of interactions between hsa_circRNA_101128 and its five potential miRNA targets. The 2D structure displayed the MRE sequence, pairing of miRNA target nucleotides and seed type. The "local AU" showed upstream and downstream of the seed sequence, black and red bars represented G/C and low accessibility as well as A/U and high accessibility of the seed, respectively, and bar height stands for the extent of accessibility. The position column showed the most likely relative MRE position on the linearized circRNA sequence. MRE, microRNA response elements discrepancy may be due to samples from different stage of active TB and small sample size in the confirmation phase. However, there is no doubt that the data suggest the clinical potential for the different cir-cRNAs profile to serve as appropriate biomarkers for diagnosis of active TB with more validation based on larger cohorts.
We further focussed on hsa_circRNA_101128, which was predicted to target let-7a and let-7f based on bioinformatics analysis, and interestingly, these miRNAs have been shown to be associated with active TB. [8][9][10] CircRNAs have been reported to be predominantly in the cytoplasm and can action as miRNA sponges to play essential regulatory roles during the development and progression of many diseases. 22,23 Considering this, we suppose that whether let-7a or let-7f has correlation with hsa_circRNA_101128. The results in the current study displayed that hsa_circRNA_101128 expression was negatively correlated with let-7a level, which indicated that hsa_circRNA_101128 might be involved in the pathogen-