Increased UBE2L6 regulated by type 1 interferon as potential marker in TB

Abstract The aim of this study is to identify potential biomarker of tuberculosis (TB) and determine its function. Differentially expressed mRNAs(DEGs) were selected from a blood database GSE101805, and then, 30 key genes were screened using STING, Cytoscape and further functionally enriched. We then found that only 6 of 13 genes related to ubiquitination (the first in the functional enrichment) were increased significantly. ROC analysis showed that UBE2L6, among 6 genes, had the highest diagnostic value, and then, we found that it also had mild value in differential diagnosis. Moreover, our analysis showed that UBE2L6 may be upregulated by type I interferon, which was further confirmed by us. In addition, we also found that UBE2L6 inhibits the apoptosis of Mycobacterium tuberculosis(Mtb)infected macrophages. Subsequently, we discovered that miR‐146a‐5p, which may target UBE2L6, is reduced in peripheral blood mononuclear cells (PBMC) and plasma of TB, and it also had certain diagnostic efficiency(AUC=0.791). In brief, we demonstrated that UBE2L6 as well as miR‐146a‐5p is a potential biomarker for TB and UBE2L6,which may also plays important role in TB by, at least, modulating Mtb‐infected macrophage apoptosis.

At present, TB diagnosis is mostly based on epidemiological characteristics, combined with microbiological tests and chest Xrays. Respiratory symptoms such as cough and expectoration and systemic symptoms such as low fever, night sweats, emaciation and weakness will be noticed; however, some patients have no obvious clinical manifestations. Although there are accepted methods for diagnosis of TB, they have well-known disadvantages, such as insufficient sensitivity (sputum smear) and long turnaround time (culture for 4 to 8 weeks) for the results. 8 In addition, Xpert provide essential information on rifampin resistance, it is less sensitive than TB sputum culture and cannot distinguish Non-tuberculous mycobacterium (NTM) infections. 9 Furthermore, the Xpert assay has limited performance for TB smear negatives and paediatric specimens (around 67% sensitivity). 10 Specific response networks were constructed by mapping disease genome-wide expression. The DEGs screened out can be used as diagnostic and differential diagnostic markers of disease and indicators of disease change.
The purpose of this study was to identify potential biomarker in blood of TB and determine its function in Mtb-infected macrophages. Bioinformatics methods were used for data analysis to identify the most influential genes in the response network of TB, and then, functional enrichment was conducted, followed by validation, diagnostic and differential diagnostic value analysis.
Subsequently, DEGs in the THP-1 cells were enriched and the role of the final hub gene was determined, its function was further studied in RAW264.7 cells. These will provide new thoughts for the diagnosis and treatment of TB. The flowchart of the study was shown in Figure 1.

| Data processing and identification of DEGs
The microarray data in this study were retrieved and downloaded from the NCBI gene expression omnibus(GEO) database(http:// www.ncbi.nlm.nih.gov/geo/) using keywords 'TB', 'Mtb' and 'Human'. After comparison, the plasma source data set GSE101805 was selected, its microarray platform is GPL16956 whose chip probes include LNCRNAs and mRNAs. The corresponding table of probe and gene was obtained by BLAST reannotation. The original data had been normalized. The GEO2R online tool was used to analyse the data, and DEGs were selected according to the following criteria:p-value<0.05,|logFC|≥2.

| Screening of hub genes and function prediction
STRING software was used to obtain the protein interaction network (PPI) of DEGs, and then, import the analysis results into Cytoscape V.3.7. Cytohubba and MCODE plug-ins were used to screen key F I G U R E 1 The flowchart of the study subnetworks and key protein modules. The hub genes were obtained by crossing the first 30 Cytohubba genes and the cluster with the highest score of MCODE. DAVID and FUNRICH software were used for Gene Ontology(GO) and Kyoto Encyclopedia of Genes and Genomes(KEGG) analysis of the hub genes.

| Functional enrichment of DEGs in cells
In order to explore the expression and function of the key hub gene in Mtb-infected macrophages, we used the GSE17477 gene data set derived from Mtb-infected THP-1 cells for subsequent analysis. The GEO2R online tool is used to process the data and filter the DEGs according to the same criteria as above. Importing them into DAVID and FUNRICH software for GO and KEGG analysis.

| miRNAs prediction of key hub gene and validation
miRWalk, StarBase and TarBase were used to predict the miRNAs which may target key hub gene, and the online software Bioxinren was used to draw the Venn diagram. We selected a PBMC database GSE131708, which included 4 patients with TB meningitis and 4 healthy individuals, to verify the reliability of the predicted miR-NAs in TB. The raw data have been normalized.

| Clinical samples collection and Quantitative reverse transcription-polymerase chain reaction (qRT-PCR)
A total of 30 plasma samples (  According to the density of 6 × 10 5 / well, the cells were inoculated in six-well plate. When cells adhered to the wall and grew exuberantly, a small amount of medium was discarded, bacterial liquid was added in proportion, and grouping was marked. After 4h, fresh medium was replaced for further culture, which was marked as infection 0. Twenty-four hours after transfection, Mtb infection was carried out at a ratio of 1:10. After 24 h of infection, the supernatant was discarded and 500μl PBS was added, followed by 5μl Hoechst and 5μl PI staining solution. After mixing, the cells were incubated at 4 degrees for 20 min. Finally, fluorescence microscope was used to observe.

| Data analysis
Statistical analysis was performed using IBM SPSS Statistics 26 software. Statistical significance of differences between and among groups was assessed using the t test. Plot uses mean with SD. Box diagram was drawn using originPro9. 1. Significant differences are indicated as follows: * p < 0.05,* * p < 0.01,* * *p < 0.001.

| DEGs and hub genes were screened
Through GEO2R analysis of the mRNA data set of healthy people and TB in GSE101805, we identified 2,019 DEGs, including 1272 upregulated genes and 747 downregulated genes ( Figure 2A). The expression level of DEGs is visualized in the form of clustering heat map ( Figure 2B).
In view of the sample type and research purpose, only the upregulated DEGs were selected for the follow-up study. We uploaded these genes to STRING software to build the protein interaction (PPI) network, which consists of 1249 nodes and 4837 edges. To find the genes that play central roles, we imported the built PPI network into Cytoscape software, whose Cytohubba and MCODE plug-ins were used to mine important nodes in the network. The first 30 hub genes found on Cytohubba plug-ins ( Figure 2C). The gene cluster with the highest score selected by the MCODE plug-in ( Figure 2D).
By combining the two algorithms, a total of 30 hub genes were screened out.

| Functional enrichment of hub genes
In order to determine the function of the 30 hub genes mentioned above, we used DAVID and FUNRICH software to conduct GO and KEGG analysis. Results of DAVID showed that these hub genes were mainly connected with protein ubiquitination (p = 3.39E −21 ),
We then performed validation and diagnostic value prediction.
Three peripheral blood gene data sets GSE34608, GSE42826 and GSE42830 were used for verification. GSE42826 and GSE42830 are the detection and verification sets of the GSE42834 series, they use the same platform and data processing methods, so the data of the two data sets are combined for analysis here. The genes with p ≥ 0.
05 and logFC<0 were excluded by GEO2R analysis, then GSE34608 database remained UBE2E1, UBE2L6, NEDD4L, RNF138, UBE2H, FBXO6 six genes, and the other two databases have UBE2L6, FBXO6, RBCK1 three genes left. We further performed ROC analysis of the above genes using the same databases ( Figure 2E,F)

| Bioinformatics analysis of UBE2L6 in Mtbinfected THP-1 cells
After the diagnostic value of UBE2L6 was clarified in TB, we wanted to further explore the mechanism of UBE2L6 in Mtb-infected macrophages. Therefore, we used the mRNA gene database GSE17477 derived from Mtb-infected THP-1 cells for analysis and study. GEO2R analysis results showed that p-value of UBE2L6 was 7.04E −07 , and the logFC value was 2.3163, indicating that UBE2L6 was significantly upregulated in Mtb-infected THP-1 cells. Next, we used DAVID and FUNRICH software to conduct GO analysis ( Figure 3A-C) and KEGG pathway prediction ( Figure 3E) on 240 DEGs. The results showed that UBE2L6 was mainly located in the cytoplasm and was important for regulating cytokine production (p = 1.33E −12 ), type I interferon production (p = 5.18E −06 ), binding to ISG15 protein (p = 0.002) and protein modification (p = 0.0077). KEGG analysis showed that it was mainly involved in the immune system (p = 2.28E −11 ) and the induction of IFN-alpha/beta pathway mediated by RIG-I/MDA5 (p = 7.95E −07 ). In short, UBE2L6 was not only involved in protein ubiquitination, but also correlated with ISGylation of proteins and other innate immune pathways. UBE2L6 and the predicted genes related to it were imported into STRING software to draw PPI network ( Figure 3D), which consisted of 36 nodes and 309 edges.

| UBE2L6 induced by type I interferon in RAW264.7 cells during Mtb infection and then inhibits apoptosis
On the basis of the above bioinformatics analysis, we want to further explore the expression of UBE2L6 in RAW264.7 cells and its relationship with type I interferon, and explore its effect on cells. Our results showed that UBE2L6 was significantly increased (p = 0.037) in Figure 4B). Neutralizing IFNAR with antibodies to inhibit type I interferon pathways, and then infecting cells, we found that UBE2L6 was down ( Figure 4B). This suggests that UBE2L6 can be induced by type I interferon pathway during Mtb infection. Next, PI/Hoechst Double Staining was used to test the effects of UBE2L6 on cells. siRNA with the most significant interference effect (p = 0.024) was screened and act on cells( Figure 4A), and then, the cells were infected for 24 h, and further stained for observation. We found that compared with NC siRNA, apoptosis was significantly increased after UBE2L6 inhibition(p = 0.009), and there was no statistical difference in cell necrosis( Figure 4C). This suggests that UBE2L6 may promote the survival of intracellular Mtb.

| Prediction and validation of miRNAs that may regulate UBE2L6
miRNA is non-coding regulatory small molecule that can play an important role by suppressing gene expression. We tried to find miRNA that interact with UBE2L6, so we used StarBase, TarBase and miRWalk to predict it, and the intersection of the results was the final predicted miRNAs ( Figure 5A), they were namely miR-145-5p, miR-212-5p, miR-1-3p, miR-130a-3p and miR-146a-5p. We selected GSE131708 gene data set to verify the reliability of above miRNAs in TB. The statistical results ( Figure 5B) indicated that compared with HC, miR-146a-5p was significantly reduced in PBMC of TB.
What about the expression of miR-146a-5p in plasma of TB patients? We collected plasma from 15 TB patients and 15 HC, and then, miR-146a-5p was quantitatively analysed by qRT-PCR.
The result indicated that the expression of miR-146a-5p was Specificity meets TPP's minimum standards, while sensitivity does not, which indicated that miRNA-146a-5p had mild clinical diagnostic value in TB.
We then examined miR-146a-5p in Mtb-infected RAW264.7 cells and exosomes from it. First, the extracted exosomes were identified, and the results showed that a clear and typical goblet like vesicle structure of exosomes was observed by TEM, and the size was relatively uniform ( Figure 5E). NTA showed that the diameter of exosomes was concentrated in the range of 80-150nm, which was in line with the standard diameter of exosomes ( Figure 5F). Western blot (WB) detected the exosome signature protein CD63( Figure 5G).
In contrast to the above results, we found that miR-146a-5p increased significantly in Mtb-infected RAW264.7 cells and exosomes from it ( Figure 5H). It means that miR-146a-5p may indirectly affect the expression of UBE2l6 by influencing other molecules. The above results indicated that the expression of miR-146a-5p was significantly different in different cells during Mtb infection.

| DISCUSS ION
The End TB Strategy of the World Health Organization (WHO) aims to reduce the annual incidence of TB to less than 10 cases per 100,000 people by 2035. 14 Therefore, it is urgent to accelerate the research on the mechanism of TB so as to carry out appropriate treatment and improve diagnostic methods to shorten the testing time, such as metabolic regulation mechanism in the host cell of Mtb and the molecular mechanism that affects the survival of bacteria In addition, Vrieling F found that oxidized low density lipoprotein supports the survival of Mtb in macrophages by inducing lysosomal dysfunction. 19 In this study, we found that selected UBE2L6 is associated with protein ubiquitination and type I interferon in TB. Protein ubiquitination is an important cellular process targeting abnormal or short-lived protein degradation and involves at least three types of enzymes: E1 ubiquitin activating enzyme, E2 ubiquitin conjugating enzyme and E3 ubiquitin ligase. 20 Studies have revealed the unique mechanism of host xenophagocytosis triggered by the direct binding of ubiquitins to pathogen surface proteins. 21  including UBA7 and UBE2L6. [23][24][25] Other studies have shown that UBE2L6 is involved in the lipolysis process in nasopharyngeal carcinoma and associated with its poor prognosis. 26 The limitation of this study is that the changes of UBE2L6 in latent TB and different stages of active TB are not taken into account, which will be needed for further in-depth research in the future.
UBE2L6 is both an E2 enzyme of ubiquitin and an E2 enzyme of ISG15. UBE2L6 has been confirmed by conjugated ISG15, and ISG15 has also been identified as a potential cancer serological marker. 27,28 Previous studies indicated that the expression of ISG15 and USP18 genes could be improved by type I interferon signalling pathway, and USP18, in turn, could play an inhibitory role in IFN signal transduction through STAT2. IFNα and ISG15 could also increase the expression of UBA7, UBE2L6 and HERC5 by activating STAT1 protein, and thus induced ISGylation of target proteins. [29][30][31] Our bioinformatics analysis showed that ISG15 and USP18 were indeed upregulation in Mtb-infected THP-1 cells. In addition, some studies found that protein ISGylation was an antagonistic system of ubiquitylation. 32,33 Another study stated clearly that ISG15 and ubiquitin could form a mixed chain to regulate the stability of proteins and the expression of exogenous ISG15 increased the overall ubiquitylation pattern. 34 38 It has also been found that miR-146a attenuates the production of TNFα in BCG-induced THP-1 cells and promotes the survival of BCG by inhibiting the production of NO in macrophages. 39,40 Thus, miR-146a plays an important role during Mtb infection. Our verification results also showed that miR-146a-5p was declined in PBMC and plasma of TB, and it is negatively correlated with UBE2L6, suggesting that miR-146a-5p may participate in protein metabolism and innate immunity by regulating UBE2L6, but the specific mechanism needs to be further investigated. However, they are expressed in the same direction in the Mtb-infected RAW264.7 cell line, suggesting that miR-146a-5p is not involved in the regulation of UBE2L6 in RAW264.7 cells or there are other unknown and complex regulatory mechanisms, which need to be further investigated. It is regrettable that we have not experimented to verify the regulatory relationship between UBE2L6 and miR-146a-5p.
In conclusion, we found that UBE2L6 was significantly up- for further study of the molecular mechanism of TB.

ACK N OWLED G EM ENTS
This work was supported by the Program for the Natural Science

CO N FLI C T O F I NTE R E S T
All authors: no reported conflicts.

AUTH O R CO NTR I B UTI O N S
Jiao Gao designed the study and wrote the paper. Chonghui Li: col-

DATA AVA I L A B I L I T Y S TAT E M E N T
The data sets used for the current study are available from the corresponding author upon reasonable request.