A comprehensive analysis of the circRNA–miRNA–mRNA network in osteocyte-like cell associated with Mycobacterium leprae infection

Background Bone formation and loss are the characteristic clinical manifestations of leprosy, but the mechanisms underlying the bone remodeling with Mycobacterium leprae (M. leprae) infection are unclear. Methodology/Principal findings Osteocytes may have a role through regulating the differentiation of osteogenic lineages. To investigate osteocyte-related mechanisms in leprosy, we treated osteocyte-like cell with N-glycosylated muramyl dipeptide (N.g MDP). RNA-seq analysis showed 724 differentially expressed messenger RNAs (mRNAs) and 724 differentially expressed circular RNA (circRNAs). Of these, we filtered through eight osteogenic-related differentially expressed genes, according to the characteristic of competing endogenous RNA, PubMed databases, and bioinformatic analysis, including TargetScan, Gene Ontology, and Kyoto Encyclopedia of Genes and Genomes. Based on these results, we built a circRNA–microRNA (miRNA)–mRNA triple network. Quantitative reverse-transcription polymerase chain reaction and western blots analyses confirmed decreased Clock expression in osteocyte-like cell, while increased in bone mesenchymal stem cells (BMSCs), implicating a crucial factor in osteogenic differentiation. Immunohistochemistry showed obviously increased expression of CLOCK protein in BMSCs and osteoblasts in N.g MDP–treated mice, but decreased expression in osteocytes. Conclusions/Significance This analytical method provided a basis for the relationship between N.g MDP and remodeling in osteocytes, and the circRNA–miRNA–mRNA triple network may offer a new target for leprosy therapeutics.


Introduction
Leprosy, caused by the bacillus Mycobacterium leprae (M. leprae), is a chronic infectious disease, with more than 200,000 cases occurring annually worldwide [1]. Previous studies have shown that these patients are at risk for being misdiagnosed with other diseases, such as nonspecified connective tissue disease and nodular vasculitis, delaying appropriate treatment [2]. The clinical manifestations of leprosy include lesions of the skin and peripheral nerves [3] and deformity and physical disability [4]. Poor oral hygiene also has been noted [5,6]. The mechanism underlying these clinical manifestations remain to be elaborated, however.
Leprosy leads to bone changes, but what causes these changes is controversial. Some researchers have suggested bone destruction as a likely explanation, and indeed, bone damage is observed in 90% of patients [7]. Furthermore, 40% of patients show absorption of terminal phalanges on radiology [8], resulting from the expression of a phosphate-regulating gene with homologies to endopeptidase on the X chromosome (PHEX) that facilitates suppression by M. leprae in osteoblasts [9]. Other authors, however, have reported finding bone formation in patients with leprosy [10]. Toll-like receptor (TLR)1 and TLR4 polymorphisms are associated with protection against leprosy [11,12], and mineralization-related protein is elevated in the condition [13]. Bone islands and osteosclerosis also have been identified by radiology in these patients [14,15]. Thus, the pathways leading to bone alteration in leprosy remain unclear.
Culturing M. leprae in vitro could crucially contribute to understanding of the bone remodeling mechanism, but this species cannot be cultured in artificial medium [16]. M. tuberculosis and M. leprae show phylogenetic proximity [17], and their bacterial products have been used in experimental analyses. Of these, muramyl dipeptide (MDP) is the minimal essential structure for the immunological effects of the cell wall [16] and can occur as N-acetyl MDP and Nglycosylated MDP (N.g MDP). N.g MDP is known to trigger an exceptionally strong immunogenic response [18].
A large class of non-coding RNAs (ncRNAs) has been discovered through large-scale genome and transcriptome studies in recent years [25], including microRNAs (miRNAs) and circular RNAs (circRNAs). The miRNAs consist of~22 nucleotides and have a crucial role in gene silencing [26]. Jorge et al. reported that miR-101, miR-196b, miR-27b, and miR-29c serve as biomarkers in the diagnosis of a subtype of leprosy [27], and Liu et al. found that through miRNA-21 expression, the host cell can ward off infection by M. leprae, preventing bacillus growth, generating natural barriers, and regulating the antibacterial pathway [28]. In addition, miR-342 mediates regulation of Col1a2 expression in bone formation [24]. The circRNAs may have potentially important roles in gene regulation, as well [29]. These ncRNAs are produced through a downstream splice-donor site covalently linked to an upstream splice-acceptor site, in a process called backsplicing [30]. CircRNAs contain many miRNAs response elements and competing binding sites, thereby reducing miRNAs and mRNAs interaction and regulating the expression of target genes at the posttranscriptional level. Thus, it is more likely to have competing endogenous (ce)RNA function and act as miRNA sponges to restrain the expression of miRNA and regulate RNA transcription [30][31][32]. The circRNA Hsa_circ_0074834 has been associated with repair of bone defects and promotion of bone mesenchymal stem cells (BMSCs) osteogenic differentiation by acting as a ceRNA for miR-942-5p [33]. However, the involvement of circRNA in leprosy is unclear.
For a better understanding of gene expression in leprosy, the potential circRNA-miRNA-mRNA triple network needs to be investigated because it represents a possible therapeutic target. In this study, we constructed this triple network, using murine osteocyte-like MLO-Y4 cells treated with N.g MDP to elucidate the activity of this pathway in leprosy.

Ethics statement
The study was approved by the Animal Ethics Committee of the Central South University (Approval No. 2021391) and conformed with the ARRIVE guidelines.

Animals
All protocols used in our animal experiments were approved by the Animal Ethics Committee of the Central South University and conformed with the ARRIVE guidelines. The C57BL/6 female mice were supplied by the Experimental Animal Center of Xiang-Ya Second Hospital. After a week of adjustable feeding, mice in the exposure group were infused with 2 μg N.g MDP dissolved in 100 μl of saline solution once a day for 10 days, and the control group was infused with 10 μl saline only. Both groups were euthanized with an overdose of anesthesia. We have previously published the details of femur removal [24]. Investigators were blinded to the group allocation during the experiment and when assessing the outcome.

mRNA/circRNA sequencing
Total RNA was extracted from the samples in Trizol reagent according to the manufacturer's instructions (Invitrogen). The concentration and purity of each RNA sample were determined using the dsDNA HS Assay kit for Qubit (12640ES76, Yeasen). The quality of the library was determined using an Agilent High Sensitivity DNA Kit (5067-4626, Agilent), and integrity and size were quantified with this kit on an Agilent 2100 Bioanalyser (Agilent Technologies, Santa Clara, CA). For the mRNA library, mRNA was purified via two rounds of hybridization to Dynal Oligo beads (N411-03, Vazyme). After depleting the samples of ribosomal RNAs, we applied fragmentation buffer (AM8740; Invitrogen) and synthesized cDNA from the fragments using random primers. Following end repair second-strand digestion and adaptor ligation, the purified fragments were PCR amplified. For the circRNA library, total RNA was subjected to ribosomal RNA depletion using the QIAseq FastSelect RNA Removal Kit (333180, QIAGEN). Remaining RNA samples were treated with RNase R (RNR07250; Epicenter) to remove linear RNAs. After preparation of cDNA, the remaining procedures were similar to those for the mRNA library. The library preparation was performed using the VAHTS mRNA-seq V3 Library Prep Kit for Illumina (NR611, Vazyme).

Bioinformatics analysis
TargetScan (http://www.targetscan.org/) software packages were used to predict the potential miRNA target of the mRNAs. The biological processes involving Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) (adjusted P value < 0.05; gene count �2) were obtained from the database for Annotation, Visualization, and Integrated Discovery (known as DAVID) (https://david.ncifcrf.gov/) (gene count �2). The differentially expressed mRNAs and proteins were visualized using the Search Tool for the Retrieval of Interacting Genes (STRING) software (https://string-db.org/). The circRNA-miRNA-mRNA interaction network was built by merging the circRNA database (provided by Guangzhou Forevergen Biosciences), differentially expressed genes, and the predicted potential miRNA. The raw data has been submitted on the GEO database (PRJNA798506).

Immunohistochemistry
Immunohistochemistry was performed as previously described [23]. Briefly, right femurs were fixed in 4% paraformaldehyde for 24-48 h, decalcified in 10% EDTA for 21 d, and embedded in paraffin. Decalcified right femur sections 5 μm thick were deparaffinized in turpentine, 3% H 2 O 2 was used to suppress endogenous peroxidase activity, and the sections were treated with 0.1% trypsin for antigen retrieval. After incubation with a CLOCK polyclonal antibody (AF0323; 1:100, Affinity Biosciences) overnight, sections were incubated with secondary antibodies. Antibodies were detected by staining with a horseradish peroxidase-conjugated rabbit anti-mouse IgG and diaminobenzidine (GTVision III Detection System/Mo&Rb Kit; Gene Tech, Shanghai, China). Specimens were counterstained with hematoxylin.

Statistical analysis
All experiments were carried out in triplicate. Data were analyzed using IBM SPSS Statistics 26. Statistical significance was considered at P < 0.05 using the Student's t-test. The results are presented as mean ± standard deviation.

Differentially expressed mRNAs and circRNAs in N.g MDP-treated osteocyte-like cell
Osteocyte-like cell was observed under a light microscope, and western blots showed DMP-1 expression in osteocyte-like cell, while not be detected in MC3T3-E1 (S1 Fig). With the cutoff criteria of a fold-change >1.5 and P < 0.05, we identified 724 differentially expressed mRNAs and circRNAs between control and N.g MDP-treated samples, with 579 up-regulated and 145 down-regulated genes in differentially expressed mRNAs and 309 up-regulated and 415 downregulated circRNAs (Fig 1). The top-10 dysregulated mRNAs (Table 1) and circRNAs ( Table 2) were summarized based on P values. Bioinformatics analysis of the interactions (Table 3)

De-regulation of Clock expression in N.g MDP-treated osteocyte-like cell in vitro and in femurs of N.g MDP-treated mice in vivo
To verify in vitro and in vivo the accuracy of the filtered results for the differentially expressed mRNAs, osteocyte-like cell was treated with or without (control) 1 μg/ml N.g MDP for 36 h, and quantitative reverse-transcription polymerase chain reaction (qRT-PCR) was used to measure and compare the expression of the identified genes in control and treated samples. We found that the expression of Clock mRNA was decreased approximately 20% (P < 0.05), but found no other significant differences between the two samples ( Fig 4C). Western blots showed an obviously decreasing amount of CLOCK protein (P < 0.05) ( Fig  4D and 4E), along with a tendency for expression of the osteogenic differentiation biomarkers Runx2 and Bglap to decline. The qRT-PCR showed that Clock mRNA was upregulated in BMSCs under osteogenic induction (Fig 4F). To illustrate the expression of CLOCK protein in vivo, mice were infused with 2 μg N.g MDP dissolved in 100 μl of saline solution once a day for

PLOS NEGLECTED TROPICAL DISEASES
circRNA-miRNA-mRNA network of osteocyte in leprosy

PLOS NEGLECTED TROPICAL DISEASES
circRNA-miRNA-mRNA network of osteocyte in leprosy 10 days, and the control group was infused with 100 μl saline only. Immunohistochemistry showed increased CLOCK protein expression in BMSCs and osteoblasts from femurs of N.g MDP-treated animals, but decreased expression in osteocytes (Fig 4G-4I).

Discussion
In this study, we analyzed RNA-seq results of mRNA and circRNA, and integrated these findings with bioinformatics analysis to construct a circRNA-miRNA-mRNA triple network. According to the two different analytical methods we applied, N.g MDP may control the process of osteoblastic differentiation through regulation of Clock expression in osteocyte-like cell, suggesting a potential therapeutic target in leprosy-related bone loss.
In the context of bone metabolism, RNA-seq has influenced every aspect of understanding of genomic function [34,35]. Using a combination of large-scale genome and transcriptome studies with bioinformatic analysis, Xu and colleagues found that targeting the SLIT3 pathway could represent a new approach to treating bone loss [36]. However, no related study has demonstrated an interaction between ncRNA and leprosy. As noted, M. leprae cannot be cultured

PLOS NEGLECTED TROPICAL DISEASES
circRNA-miRNA-mRNA network of osteocyte in leprosy in vitro, representing an obstacle to researching bone remodeling mechanisms in this disease [16,37,38].
MDP is a component of the cell wall and the minimal essential structure for eliciting immunological effects [16]. To further investigate the relationship between MDP and ncRNA, we developed an in vitro MDP-treated osteocyte-like cell model and identified several differentially expressed mRNAs and circRNAs. However, because we could not verify all of these differentially expressed genes in experimental analysis, we used TargetScan and the circRNA database (provided by Guangzhou Forevergen Biosciences Co., Ltd.) to predict the target miRNA and construct a circRNA-miRNA-mRNA triple network. The circRNA-miRNA-mRNA triple network was wildly applied in research, Chen X et al. found that Cir-cRNA_28313/miR-195a/CSF1 could regulate osteoclastic differentiation in OVX-induced bone absorption in mice [39], and Shen WX et al. revealed that CircFOXP1/miR-33a-5p/ FOXP1 could promote osteogenic differentiation and CircFOXP1 can be used as a potential osteoporosis therapeutic target [40]. Based on the predicted results, we screened out the differentially expressed genes using a fold-change cutoff of >2.0 for mRNA [41]. We narrowed the list to six differentially expressed mRNAs using a PubMed search and to two differentially expressed mRNAs using GO and KEGG analysis. We confirmed decreased expression of the Clock gene, one of the two identified with the latter approach, but that expression of Rcan2, Cacnb4, Plag1, Lpar6, Npnt, Lin28b, and Rora did not differ between treated and untreated samples ( Fig 4C). Thus, this circRNA-miRNA-mRNA triple network could provide a new analytical approach to elucidating the relationship between N.g MDP and bone remodeling.
The circadian clock is a key factor in the circadian system, which is associated with cell function, metabolic state, and life expectancy [42,43]. Bone remodeling could be subject to circadian regulation [44]. Previously research has indicated that mice lacking Clock or brain and muscle Arnt -like protein 1 (Bmal1) have a reduced bone mass [45]. A recent study demonstrated that the circadian clock could regulate bone formation through transcriptional control of the 1,25-dihydroxy-vitamin D3 receptor PDIA3 [46]. Decreasing of the bone mineral density and increasing risk of patients developing osteoporosis would happen, if the circadian rhythm has disordered [47]. Bglap and Runx2 also displayed rhythmicity with the expression of Clock gene [48]. BMAL1 could influence numbers of key factors in skeletal development, such as Runx2/Osterix and HIF1α-VEGF signaling pathway [49]. Our results showed decreased expression of Clock, along with expression of the osteogenic biomarkers Runx2 and Bglap (Fig 4D), while the expression of Clock was increased in BMSCs under osteogenic induction. Yuan HP et al. found that intraperitoneal injection of MDP in ApoE −/− mice could reduce alveolar bone loss which expose in Porphyromonas gingivalis [50]. Park OJ et al. revealed that MDP could induce bone formation and osteoblast activation by upregulating the expression of Runx2 and ALP in BMSCs and preosteoblasts, and micro-CT showed increased trabecular bone in the femur of MDP-administered mice [51]. However, the regulatory effect of MDP on osteocyte still not be reported. Immunohistochemical analysis indicated that expression of the CLOCK protein was significantly increased in murine femur BMSCs and osteoblasts, while decreased in osteocytes in exposed mice (Fig 4E). The paradoxical results are possibly because of negative feedback regulation of osteocytes to osteoblasts and BMSCs [52], and indicate that different cells induce different responses in vivo and in vitro. Previous study has generated Irs-1-deficient Irs-1 smla/smla mice, and in BMSCs derived from Irs-1-null mice, the expression of COL1A2 was decreased, then, osteocytes promote osteogenesis through negative feedback regulation [53]. However, the candidate miRNA and circRNA involving in our RNA-seq, miR-322-5p and mmu_circ_0006240, has not been reported. In general, these results suggest a crucial role for Clock in bone remodeling.
In summary, based on the results of transcriptome studies and bioinformatic analysis, we have characterized a circRNA-miRNA-mRNA triple network and screened out the key factors in the N.g MDP-treated osteocyte-like cell. Our approach potentially represents a new analytical method for elucidating the mechanism of bone remodeling in leprosy, and the circRNA-miRNA-mRNA could filter the important factors that might change the presentation of the diseases.