miRNA Expression Profiles and Potential as Biomarkers in Nontuberculous Mycobacterial Pulmonary Disease

Pulmonary disease (PD) due to nontuberculous mycobacteria (NTM) is increasing globally, but specific biomarkers for NTM-PD have not been established. As circulating miRNAs are promising biomarkers for various diseases, we investigated whether miRNAs have potential as NTM-PD biomarkers. Sera from 12 NTM-PD patients due to Mycobacterium avium, M. intracellulare, M. abscessus, or M. massiliense and three healthy controls were initially evaluated via small RNA sequencing. Multiple miRNAs showed significant differences in expression in patients compared to in healthy controls, with some expression differences unique to PD caused by a specific mycobacterial species. Notably, 14 miRNAs exhibited significant expression differences in PD associated with all four mycobacteria. Validation by quantitative reverse-transcription-PCR in an additional 40 patients with NTM-PD and 40 healthy controls confirmed that four differentially expressed miRNAs (hsa-miR-484, hsa-miR-584-5p, hsa-miR-625-3p, and hsa-miR-4732-5p) showed significantly higher serum expressions in NTM-PD patients than in controls. Receiver operating characteristic curve analysis of these four miRNAs supported the discriminative potential for NTM-PD and their combination provided an improved diagnostic value for NTM-PD. Furthermore, bioinformatics analysis revealed their 125 target genes, which were mostly associated with immune responses. Collectively, this study identified four miRNAs as potential biomarkers for NTM-PD and provided insight into NTM-PD pathophysiology.


patients.
In total, 95 participants were recruited in this study, including 52 patients with NTM-PD and 43 healthy volunteers ( Table 1). The median age of the patients was 56 years, and 45 (87%) were female. The median body mass index of the 52 patients was 20.8 kg/m 2 , and 46 (89%) were never smokers. Among the patients with NTM-PD, the most common underlying respiratory disease was bronchiectasis (n = 48, 92%), followed by previous treated TB (n = 17, 33%). The median value of white blood cell count and erythrocyte sedimentation rate was 6040/µL and 36 mm/hr, respectively. Out of the 52 patients with NTM-PD, 48 (92%) had nodular bronchiectatic form, and 4 (8%) had the fibrocavitary form ( Supplementary Fig. S1). The median forced expiratory volume in 1 second (%) on pulmonary function test at diagnosis was 79%. In the 43 healthy controls, the median age was 48 years, and 29 (67%) patients were female.
Illumina-based small RNA sequencing of serum miRNAs. Serum   each), and three healthy individuals were examined by Illumina small RNA sequencing in the discovery phase. Averages of 11,217,363 and 18,556,427 reads of RNAs ranging from 18 to 30 nucleotides were obtained from pooled serum samples of patients with NTM-PD and healthy controls, respectively. The length distribution of clean sequences in the reference genome was determined. Length distribution analysis of serum pools from both patients with NTM-PD and healthy controls revealed that most reads were in the range of 18-24 nt, which is consistent with the common length of miRNAs 28,29 . The three pools of serum samples contained small RNAs of various lengths ( Supplementary Fig. S2). Thereafter, bioinformatics analysis was performed to investigate the small RNA species and sequencing frequencies.
In serum derived from patients with NTM-PD and healthy controls, multiple and heterogeneous small RNA species, including miRNAs, long intergenic noncoding RNA, ribosomal, small nucleolar, and small nuclear RNA were identified (Table 2). In serum samples of patients with NTM-PD and healthy controls, miRNAs accounted for 39.7% and 44.5%, respectively, of the total amount of small RNAs (sequencing reads), and there was no significant difference between the two groups. We analyzed all clean reads using mirDeep software to identify known and novel miRNAs. In total, 467 and 407 known miRNAs were identified among patients with NTM-PD and healthy controls, respectively (Table 2). Novel miRNA precursors were predicted, and an average of 262 miRNAs in patients with NTM-PD and 193 miRNAs in healthy controls were identified (Table 2).
In addition, differentially expressed miRNAs were observed in the serum miRNA profiles of patients with NTM-PD relative to in healthy controls (Fig. 1). We identified 148 differentially expressed miRNAs based on two-fold changes and p-values less than 0.05. Of these 148 miRNAs, 70 miRNAs showed significant differences in expression in patients with M. avium-PD, 56 miRNAs in patients with M. intracellulare-PD, 46 miRNAs in patients with M. abscessus-PD, and 32 miRNAs in patients with M. massiliense-PD compared to in healthy controls (Supplementary Tables S1-4, Fig. 2a-d). However, only 14 of these miRNAs were differentially expressed in all four patient groups, i.e., in PD due to any of the four NTM organisms (M. avium, M. intracellulare, M. abscessus, and M. massiliense) (Fig. 2e,f).
Differential expression of miRNAs via qRT-PCR analysis. As the 14 miRNAs with differential expression in all four patient groups were considered as potential biomarkers, their expression was verified by qRT-PCR analysis. In the selection phase, the 14 candidate miRNAs were measured in a separate set of individual serum samples from 12 patients with NTM-PD and three healthy controls, i.e., the number of the same subjects used for small RNA sequencing. miRNAs with a mean fold-change (NTM-PD/healthy controls) of >2 and p-values < 0.05 were selected for further analysis. The analysis indicated that five miRNAs (hsa-miR-423-5p, hsa-miR-484, hsa-miR-584-5p, hsa-miR-625-3p, and hsa-miR-4732-5p) of the 14 candidates showed significant differences in expression between patients with NTM-PD and healthy controls (Fig. 3).    www.nature.com/scientificreports www.nature.com/scientificreports/ has-miR-4732-5p (p = 0.014), showed significantly higher serum expression levels in patients with NTM-PD than in healthy controls (Fig. 3).

Go analysis and KEGG pathway analysis of target genes. Gene ontology (GO) and the Kyoto
Encyclopedia of Genes and Genomes (KEGG) pathway analysis were applied to further classify the identified genes into the functional groups. GO enrichment analyses were performed to classify 125 target genes based on biological processes, molecular function, and cellular component (Fig. 6a). In the case of biological processes, the genes were classified into forty-six categories, of which the most significant term was regulation of transcription and DNA-template (GO: 0006355). Based on cellular components, the genes were classified into ten categories, of which the most significant term was mitochondrial lumen (GO: 0005759). Based on molecular function, the

Discussion
This study was conducted to compare the expression profiles of circulating serum miRNA from NTM-PD and healthy controls. In the discovery phase, 70 miRNAs in patients with M. avium-PD, 56 miRNAs in patients with M. intracellulare-PD, 46 miRNAs in patients with M. abscessus-PD, and 32 miRNAs in patients with M. massiliense-PD showed significant differences in expression relative to levels in healthy controls. Differential expression of 14 of these miRNAs was common to PD caused by any of the four organisms. In the validation phase, using qRT-PCR analysis of sera from additional patients and controls, we confirmed that four miRNAs (hsa-miR-484, hsa-miR-584-5p, hsa-miR-625-3p, and hsa-miR-4732-5p) of the 14 miRNAs were differentially expressed in patients. We further determined that 125 gene targets were shared by the candidate miRNAs and were associated with NTM-PD on TargetScan analysis. The present study is the first to demonstrate that hsa-miR-484, hsa-miR-584-5p, hsa-miR-625-3p, and hsa-miR-4732-5p are associated with NTM-PD and are potential biomarkers for NTM-PD. This study is valuable because multiple patients with NTM-PD and disease due to several NTM organisms were examined.
In the present study, four serum miRNAs (hsa-miR-484, hsa-miR-584-5p, hsa-miR-625-3p, and hsa-miR-4732-5p) showed increased expression in patients with NTM-PD compared to in healthy controls. Interestingly, these four miRNAs were previously implicated in infection with Mycobacterium bovis bacillus Calmette-Guérin (BCG) and immune diseases. miR-484 is upregulated in exosomes of BCG-infected macrophages and targets mitochondrial cleavage protein 1 to modulate the intermediate metabolic pathway 30 . Various assays have revealed that miR-484 upregulation clearly influences cell migration and cellular proliferative capacity, and miR-484 is also tumorigenic, promoting the progression of non-small cell lung cancer 31 . miR-584 is reportedly upregulated in breast cancer cells, enhancing nuclear factor kappa B signaling, and inhibiting TGF-β signaling 32 . TGF-β has several inhibitory effects on T cells, B cells, macrophages, and other cells, and increased TGF-β levels are associated with protection and/or recovery from autoimmune diseases 33 . Previous studies reported that miRNA-625-3p in the blood facilitated the detection of malignancy in benign lung tumors 34 , and miRNA-625-3p is reportedly upregulated in the urine of pulmonary TB patients and is an excellent diagnostic urinary biomarker 35 . miR-4732-3p is associated with breast cancer and stomach cancer and is specifically expressed in exosomes released from macrophages after infection with BCG 20,36 . We hypothesized that differential expression of these miRNAs in patients with NTM-PD and healthy controls would result in differences of target mRNA expression; thus, additional experimental validation is further required.
In the present study, 125 genes were predicted to be targeted by the four miRNAs, including nine genes related to M. tuberculosis infection. Our results revealed that various genes associated with the innate immune response, cell proliferation, and apoptosis are potential targets of hsa-miR-484, hsa-miR-584-5p, hsa-miR-625-3p, and hsa-miR-4732-5p. Our results further showed that NFAT5 is a target gene for all four miRNAs; NFAT5 plays an important role in modulating HIV-1 replication by M. tuberculosis during co-infection, by directly interacting with the viral promoter 37 . Hsa-miR-484, hsa-miR-584-5p, and hsa-miR-4732-5p also targeted IL-6, a pleiotropic cytokine that plays a central role in host defense 38 . These results suggest that the four miRNAs play an important role in host-pathogen interactions.  www.nature.com/scientificreports www.nature.com/scientificreports/ We performed comprehensive GO enrichment and KEGG pathway analysis of the target genes of the differentially expressed miRNAs. GO analysis indicated that these genes are involved in immune responses and oxidation-reduction processes and in the regulation of tumor necrosis factor-α. KEGG analysis revealed the involvement of these genes in cell growth, migration, and proliferation and in pathways such as Hippo signaling, Wnt signaling, p53 signaling, and TGF-β signaling. Thus, these miRNAs may influence NTM infection, as they affect immune cell development. In addition, most proteins encoded by the target genes are involved in mitogen-activated protein kinase signaling, Wnt signaling, and TGF-β signaling in NTM-PD and may be useful in elucidating the pathogenesis of NTM-PD.
To further evaluate the diagnostic value of hsa-miR-484, hsa-miR-584-5p, hsa-miR-625-3p, and hsa-miR-4732-5p, we performed a more detailed analysis using ROC curves. The present study found that the ROC curve of the combined four differentially expressed miRNAs showed an optimal discriminative ability with an AUC of 0.998 (95% CI = 0.949-1.000), and diagnostic sensitivity and specificity of 95 and 100%, respectively. Our data showed that the diagnostic value of these four miRNAs combination was higher than that of the individual miR-NAs, indicating that the combinational signature of these four miRNA could serve as discriminating markers for NTM-PD.
There are several limitations to the current study. First, the number of samples in this study was small. Second, the present study did not include patient subjects with other lung diseases as controls for NTM-PD. To further validate the four miRNAs in discriminating NTM-PD from other lung diseases, appropriate PD control groups are required to confirm whether our findings are unique to NTM-PD. Finally, the combination sets of miRNAs target prediction were solely based on TargetScan database. Thus, the use of efficient genetic algorithms for variable selection in logistic regression model, rather than a linear model, may have advantages to select the clinically optimal markers among all extracted data sets, more accurately 39 .
In conclusion, to our best knowledge, this is the first study to determine the serum miRNA profiles in NTM-PD through deep sequencing. We identified 14 differentially expressed miRNAs in patients with NTM-PD compared to in healthy individuals, including nine upregulated miRNAs and five downregulated miRNAs. The present study verified the relative expression levels of hsa-miR-484, hsa-miR-584-5p, hsa-miR-625-3p, and hsa-miR-4732-5p, indicating their potential as biomarkers for NTM-PD. GO enrichment and KEGG analysis showed that the target genes of the four differentially expressed miRNAs were mainly involved in the immune response. These results indicate the presence of a serum-specific miRNA signature in patients with NTM-PD, which may modulate host antibacterial pathways in response to NTM infection and can be considered a potential molecular tool for diagnosing NTM-PD. The ROC curves were used to show the discriminative ability of the selected four miRNAs with sensitivity (true-positive rate) and 1-specificity (false-positive rate). The x-axis shows 1-specificity, and the y-axis shows sensitivity. Multivariate logistic regression analysis was performed, and ROC curves were generated to evaluate the ability of the chosen miRNAs to distinguish NTM-PD patients from controls. (a) has-miR-484, (b) has-miR-584-5p, (c) has-miR-625-3p, (d) has-miR-4732-5p, and (e) has-miR-panel. AUC, area under the curve; Panel 1, has-miR-484+ has-miR-584-5p+ has-miR-625-3p+ has-miR-4732-5p; Panel 2, has-miR-484+ has-miR-584-5p+ has-miR-4732-5p; Panel 3, has-miR-484+ has-miR-625-3p+ has-miR-4732-5p; Panel 4, has-miR-584-5p+ has-miR-625-3p+ has-miR-4732-5p.

Scientific RepoRtS |
(2020) 10:3178 | https://doi.org/10.1038/s41598-020-60132-0 www.nature.com/scientificreports www.nature.com/scientificreports/ showing overlap of gene targets of hsa-miR-484, hsa-miR-584-5p, hsa-miR-625-3p, and hsa-miR-4732-5p. 125 genes were screened to be targeted by the four miRNAs. (b) miRNA-gene network for hsa-miR-484, hsa-miR-584-5p, hsa-miR-625-3p, and hsa-miR-4732-5p. The network was built using 125 targeted genes and predicted interactions from the TargetScan databases. Red represents miRNAs, and green represents genes; their relationship is represented by the edges. (c) Protein-protein interaction analysis performed using the Search Tool for the Retrieval of Interacting Genes database revealed that the target genes of hsa-miR-484, hsa-miR-584-5p, hsa-miR-625-3p, and hsa-miR-4732-5p may interact with each other. The disconnected nodules are hidden. Bold lines indicate stronger association. This study was performed in several phases to identify serum miRNAs that may be useful as markers of NTM-PD ( Supplementary Fig. S3). Initially, in the discovery phase, screening was conducted by Illumina small RNA sequencing to select miRNAs with expression patterns that differed between patients with NTM-PD and control subjects. In this step, 15 serum samples were tested to obtain miRNA profiles, with the sera including three samples from patients with NTM for each of the four different pathogens, for a total of 12 patient samples, as well as three samples from control subjects. A panel of miRNAs showing differential expression in patients with NTM-PD and common to PD caused by all four pathogens was subjected to further testing by quantitative reverse-transcription (RT)-polymerase chain reaction (PCR) analysis using sera from additional patients and control subjects.
RnA extraction. Total RNA was isolated from sera using a miRNeasy Mini Kit (Qiagen, Hilden, Germany) following the manufacturer's protocol. Each sample was eluted in RNase-free water. RNA quantity and purity were measured using a NanoDrop ™ ND-2000 (Thermo Fisher Scientific, Waltham, MA, USA), and RNA integrity was assessed using an Agilent Bioanalyzer 2100 (Agilent Technologies, Santa Clara, CA, USA). The RNA samples were stored at −80 °C until further analysis.
Illumina high-throughput sequencing. Small RNA libraries were constructed using the NEXTFlex Small RNA Library Prep Kit (Illumina, San Diego, CA, USA) according to the manufacturer's instructions. Briefly, an RNA sequencing library was prepared by cDNA amplification, end-repair, adenylylation of 3′ ends, adapter ligation, and amplification. The libraries were sequenced using the HiSeq 2500 system (Illumina) as 50-base pair single reads.

Quantification of miRNAs by qRT-PCR. Reverse transcription reactions were performed using a
Taq-Man miRNA Reverse Transcription kit (Applied Biosystems, Foster City, CA, USA) and miRNA-specific stem-loop primers in accordance with the manufacturer's instructions. Each reaction mixture for qRT-PCR contained 2.5 μL of 2 × TaqMan Universal PCR Master Mix without AmpErase UNG, 0.25 μL miRNA-specific primer/probe mix, and 2.25 μL diluted RT product (1:15) in a total volume of 5 μL. Reactions were carried out using the following thermal cycling parameters: 95 °C for 10 min, followed by 40 cycles of 95 °C for 15 s, and 60 °C for 1 min, followed by holding at 4 °C. Raw data were analyzed using SDS Relative Quantification Software version 2.2.3 (Applied Biosystems), generally using the automatic Ct setting for assigning baseline and threshold values for Ct determination. The expression level of each target miRNA was normalized to that of miR-16 (internal control). miRNA-gene network construction. TargetScan (http://genes.mit.edu/targetscan/index.html) was used to predict genes targeted by the obtained miRNAs. Target genes associated with NTM-PD were screened using the NCBI database (http://www.ncbi.nlm.nih.gov/pubmed/) and Cytoscape (Cytoscape Software, Version 2.8.2, Seattle, WA, USA). The miRNA gene network was constructed using Cytoscape to analyze miRNA-mRNA interactions.
GO analysis and KEGG pathway analysis. GO enrichment was performed for analysis of the biological process, cellular component, and molecular function of target genes based on the GO database (http://www. geneontology.org/). KEGG pathway analysis was performed to identify the enriched pathways of target genes based on the KEGG database (http://www.genome.jp/kegg/). In addition, KEGG pathway annotation of target genes was found using data from the Database for Annotation, Visualization and Integrated Discovery (DAVID) (http://david.abcc.ncifcrf.gov/).

Statistical analysis.
Data are presented as the number (%) for categorical variables and median (interquartile range) for continuous variables. The data were statistically analyzed using SPSS software, version 17.0 (SPSS, Inc., Chicago, IL, USA) or the programs built into the software. Briefly, two-tailed unpaired t-tests were performed for one between-group comparison, and one-way analysis of variance was performed for multiple-group comparisons. Receiver operating characteristic (ROC) curve analysis was performed to evaluate the discriminative factor of the miRNAs. The AUC with 95% confidence interval (CI) were calculated to determine the specificity and sensitivity of discriminative markers of NTM-PD. ROC curve and logistic regression model were calculated by using the MedCalc Software (Version 19.1.13, Belgium). A default p < 0.05 was considered statistically significant.

Data availability
All data generated or analysed during this study are included in this published article and its Supplementary Information Files.