Identification of the Association Between Toll-Like Receptors and T-Cell Activation in Takayasu’s Arteritis

To explore the relationships between Toll-like receptors (TLRs) and the activation and differentiation of T-cells in Takayasu’s arteritis (TAK), using real-time fluorescence quantitative polymerase chain reaction, mRNA abundance of 29 target genes in peripheral blood mononuclear cells (PBMCs) were detected from 27 TAK patients and 10 healthy controls. Compared with the healthy control group, the untreated TAK group and the treated TAK group had an increased mRNA level of TLR2 and TLR4. A sample-to-sample matrix revealed that 80% of healthy controls could be separated from the TAK patients. Correlation analysis showed that the inactive-treated TAK group exhibited a unique pattern of inverse correlations between the TLRs gene clusters (including TLR1/2/4/6/8, BCL6, TIGIT, NR4A1, etc) and the gene cluster associated with T-cell activation and differentiation (including TCR, CD28, T-bet, GATA3, FOXP3, CCL5, etc). The dynamic gene co-expression network indicated the TAK groups had more active communication between TLRs and T-cell activation than healthy controls. BCL6, CCL5, FOXP3, GATA3, CD28, T-bet, TIGIT, IκBα, and NR4A1 were likely to have a close functional relation with TLRs at the inactive stage. The co-expression of TLR4 and TLR6 could serve as a biomarker of disease activity in treated TAK (the area under curve/sensitivity/specificity, 0.919/100%/90.9%). The largest gene co-expression cluster of the inactive-treated TAK group was associated with TLR signaling pathways, while the largest gene co-expression cluster of the active-treated TAK group was associated with the activation and differentiation of T-cells. The miRNA sequencing of the plasma exosomes combining miRDB, DIANA-TarBase, and miRTarBase databases suggested that the miR-548 family miR-584, miR-3613, and miR-335 might play an important role in the cross-talk between TLRs and T-cells at the inactive stage. This study found a novel relation between TLRs and T-cell in the pathogenesis of autoimmune diseases, proposed a new concept of TLR-co-expression signature which might distinguish different disease activity of TAK, and highlighted the miRNA of exosomes in TLR signaling pathway in TAK.


INTRODUCTION
Toll-like receptors (TLRs) are a type of pattern recognition receptor that can initiate multiple immune responses by combining with pathogen-associated molecular patterns (PAMPs) and damageassociated molecular patterns (DAMPs). A growing number of studies have been demonstrated that TLRs take part in the pathogenesis of many kinds of autoimmune diseases (AIDs) through several mechanisms (1)(2)(3). Recently, NI-0101, an antitoll-like receptor 4 monoclonal antibody was used in RA, which was the first clinical trial to target TLRs to treat autoimmune diseases (4). The dearth of knowledge on the role of TLRs in the pathogenesis of Takayasu's arteritis (TAK) is especially notable when compared to the relatively higher number of studies on their role in that of rheumatoid arthritis (RA), systemic lupus erythematosus (SLE), and multiple sclerosis (MS).
TAK is a primary large vessel vasculitis mainly affecting the aorta and its major branches. Patients with onset TAK are typically female, and more than 90% of them are under 30 years old. Annual TAK incidence rates are estimated to be 1.5 cases per million in Japan, and 0.2~2.6 cases per million in Europe and North America (5). The autoimmune inflammation of TAK appears to be dominantly driven by T-cells (6). Studies have shown that there are numerous associations between TLRs and the activation and differentiation of T-cells. For instance, TLRs expressed in innate immune cells, such as DCs and macrophages, can regulate the activation and differentiation of T-cells. On the other hand, TLRs expressed in T-cells can influence T-cells more directly. For instance, TLR1 and TLR9 are highly expressed in CD4 + T cells, and CD8 + T cells have abundant TLR3 and TLR4 expression (7). TLR7-mediated suppression of Th17 cells does not require dendritic cell involvement (8), TLR2 signaling alters the transcriptional program of differentiating and increases the proliferation of Th17 cells (9), TLR8 signaling suppresses glucose uptake and metabolism in Treg cells (10), and TLR2 signaling enhances the movement of Treg cells (11), which all act on T-cell directly. Thus, we hypothesized that TLRs regulate the activation and differentiation of T-cells in TAK. However, there have been only two studies in the literature on TLRs in the pathogenesis of TAK in PubMed database (12,13). Kabeerdoss et al. found that the higher mRNA expression of TLR4 and its ligand S100s in peripheral blood mononuclear cells (PBMCs) of TAK patients compared to healthy controls and that after being stimulated with TLR4 ligand (12), PBMCs from TAK patients had a higher mRNA expression of IL-1b and IL-1R2 compared to that of HC (13). Taken together, it has been currently unknown whether TLRs are related to the disease activity or the activation and differentiation of T-cells in the pathogenesis of TAK.
In organisms, genes form molecular networks, these molecular networks tend to be modular, and similar modules combine to function (14)(15)(16). If a network has a high clustering coefficient, it suggests the presence of local cliques or clusters of connected molecules (17). Genes with high expression similarity or linked by the shortest path are often involved in the same biological pathway or are subjected to shared regulatory pathways (18)(19)(20)(21). The identification of stable and reliable gene co-expression networks is essential to unravel the interactions and functional correlations between genes (22). Analysis of the gene co-expression network is one of basic approaches currently adopted by research on the relations between two clusters of genes (14,15).
Additionally, in this study, we defined 'co-stimulatory molecules' as both negative and positive co-stimulatory molecules (23). After T cell receptor (TCR) activation, the fate of the T cells is controlled by signals from T cell co-stimulatory molecules and cytokines largely.
The aim of this study was to determine whether a relation exists between TLRs and the activation and differentiation of Tcells in TAK and to analyze the key molecules that play an important role in this regulation. This study may provide experimental evidence for targeting TLRs in the treatment of TAK. Figure 1 Summarized the basic workflow of this study.

Gene Function Annotation
Universal Protein Resource (UniProt) SwissProt database was used to annotate gene functions (24).

Functional Enrichment Analysis
Functional enrichment analysis was performed by the Metascape webserver (25) (https://metascape.org/), and the pathway databases consisted of Gene Ontology (GO) (26), KEGG (27), and WikiPathways (28). Enrichment analysis was performed by using a hypergeometric test with a p-value cutoff of 0.01. The enrichment results were visualized using the R Package ggplot2, the R Package Circlize, and Cytoscape V. 3.8.2.

Networks Based on Public Databases
The gene co-expression network was constructed using the COEXPEDIA database (29) (www.coexpedia.org), and the protein-protein interaction (PPI) networks were constructed using the STRING database (30) (https://www.string-db.org/).

Patients
Treated TAK patients fulfilling the 1990 ACR criteria (31) were enrolled. And we assessed the disease activity of TAK by the 1994 NIH criteria (32), which included the following. 1. Systemic features, such as fever, musculoskeletal (no other cause identified). 2. Elevated ESR. 3. New onset or aggravated features of vascular ischemia or inflammation, such as claudication, diminished or absent pulse, bruit, vascular pain (carotodynia), asymmetric blood pressure in either upper or lower limbs (or both). 4. Typical angiographic features.
If a TAK patient had two or more features, he was defined as "active TAK patient"; otherwise, we diagnosed the patient was FIGURE 1 | The basic workflow of this study. (A) mRNA abundance of 29 target genes in PBMCs was measured using RT-qPCR, and these genes were highly enriched in TLR signaling pathways and the function of the activation and differentiation of T-cells by functional enrichment analysis. (B) Sample clustering analysis suggested that TAK patients had a different expression pattern of the selected genes from the healthy controls. Gene expression differential analysis demonstrated the upregulation of TLR signaling pathway in TAK, but no substantial difference in TLR signaling pathway was found between the inactive-treated TAK patients and the active-treated TAK patients. (C) Correlation analysis showed that the inactive-treated TAK group exhibited a unique pattern of inverse correlations between the TLRs gene cluster and the gene cluster associated with T-cell activation and differentiation, while the active-treated TAK group did not. And TLRs and their correlation cluster could distinguish active patients from inactive patients in TAK. (D) Dynamic gene co-expression network was constructed. TLR-co-expression signature of different stages was observed, and the degree of functional association between the other genes with TLRs was assessed. (E) To account for the co-expression in the greatest cluster of the inactive-treated TAK group, functional enrichment analysis, miRNA database prediction, and miRNA sequencing of plasma exosomes was performed, which indicated that miRNAs might play an important role in the cross-talk between TLR and T-cell in TAK patients. at remission stage, and the patient was defined as "inactive TAK patient". A patient that never received any TAK medication was defined as "untreated TAK patient", and a patient that was under treatment was defined as "treated TAK patient". Both the untreated and the treated patients were classified into the inactive group and the active group.
Written informed consent was obtained from all participants and the study was performed in accordance with the Declaration of Helsinki. And this study was approved by the Institutional Review Board of Peking Union Medical College Hospital, Beijing, China (S-478).

Collection and Processing of Human Blood Samples
PBMCs were isolated from patients by density-gradient centrifugation. Total RNA was prepared from the PBMCs using Trizol reagent (15596026, Thermo Fisher Scientific) (33). The RNA samples were diluted in RNase-free water, denatured at 65°C for 10 min. RNA concentration and purity were determined spectrophotometrically, and the RNA integrity was verified by denaturing RNA gel electrophoresis.

Real-Time Fluorescence Quantitative Polymerase Chain Reaction (RT-qPCR)
RNA was reverse transcribed using the PrimeScript ™ RT reagent Kit with gDNA Eraser (RR047A, Takara). Genomic DNA (gDNA) was eliminated at 42°C for 2 min. Reverse transcription was performed using the following conditions: 37°C for 15 min, 85°C for 5 sec. RT-qPCR reactions were performed with the iTaqTM Universal SYBR ® Green Supermix (725124, Bio-Rad) and primers were listed in Supplementary Table 1. The temperature cycle parameters in an Applied Biosystem 7900HT1 were: 95°C for 30 sec and 40 cycles of 95°C for 30 sec, 56°C for 30 sec and 72°C for 40 sec followed by a hold at 72°C for 40 sec. Gene expression was calculated using the 2 −DDCq method. Melting curve analysis was performed from 65 to 95°C.

The Selection of Reference Genes
A list of 9 genes previously reported that were stably expressed in human PBMCs from the literature was compiled, including ACTB, b-glucuronidase, B2M, GAPDH, HPRT1, PGK1, RPL13A, SDHA, and YWHAZ. We Assessed the gene expression stability and selected the most appropriate housekeeping genes for each analysis using the geNorm (34), the NormFinder (35), and the BestKeeper software (36). As a result, B2M and YWHAZ were used as the internal reference genes for the comparison of healthy controls and untreated TAK patients, B2M and SDHA were used for the comparison of healthy controls and treated TAK patients and the comparison of the active-treated and the inactive-treated TAK patients, and HPRT1 and YWHAZ were used for the comparison of the untreated and the treated TAK patients and correlation analysis.

Dynamic Gene Co-Expression Network
We constructed the gene co-expression network of the healthy control group, the untreated TAK group, the treated TAK group, the inactive-treated TAK group, and the active-treated TAK group based on the qPCR dataset, respectively. The Pearson correlation has an advantage in predicting the interaction of molecules and is widely used as a measure of gene co-expression in many public databases (37,38), so it was adopted the main approach in the network analysis and the results were available in the main text. However, Spearman correlation has an advantage in revealing the functional associations, so it was adopted as a complementary approach and the results are available in the Supplementary Information (39). Statistical analysis was performed using IBM SPSS statistic V.23 (Armonk, New York, USA). A p-value of less than 0.05 was considered significant. Cytoscape V. 3.8.2 was used to edit the networks.

Unsupervised Hierarchical Clustering Analysis and the Co-Expression Cluster
Samples were clustered based on the Pearson correlation coefficient for the profile of those 29 genes using the R package ggcorrplot. Hierarchical clustering of samples was carried out using Pearson correlation and Spearman correlation respectively (38). Data were partitioned into five clusters by cutting the clustering tree at the height of 1.0 and the co-expression clusters were defined which consist of similarly expressed genes. The clustered heatmap and hierarchical clustering trees were performed using the R package ggcorrplot and the R package ggplot2 by R V.4.1.0. A p-value of less than 0.05 was considered significant.

Scoring System for the Assessment of Disease Activity
The disease activity score was calculated using the algorithm as published before (40). To assess the disease activity of TAK, highly co-expressed gene pairs (defined as |r| > 0.73, p < 0.01) of the inactive-treated were selected. Regression equations were calculated by linear least squares. To qualify how well the data of an individual sample fit a regression line, we calculated the relative error (RE) as the ratio of the absolute error (AE) and the observed value of gene expression level. The AE is the absolute value of the difference between the predictive value and the observed value.
M-value was introduced which was equal to the RE: M = j (y predictive value − y observation )=y observation j : (40) Qualified the data of each sample to fit every regression line one by one. The M-value of each active or inactive patient was calculated. Then the patients were ranked by M-value. The threshold was set using the Youden index and was adjusted when appropriate. A score of 1 was given if a patient had an Mvalue less than the threshold, 0 if a patient had an M-value greater than the threshold. The total score of a patient was the sum of each score of M-value. The threshold of the total score was set using the Youden index. Detailed protocols are available in ref.
Evaluate the Closeness of Gene-to-Gene Functional Relationships Between TLRs and the Other Genes To predict the functional relationship between TLRs and the other genes through quantification, the following metrics were employed, and its rationale has been published (19)(20)(21)(41)(42)(43). Take BCL6 for an example.
1. The number of TLR-genes that were co-expressed with BCL6. 2. The number of common neighborhoods with TLRs. In the metric, all TLRs were processed as a whole-body corresponding to one gene. 3. The number of TLRs in the co-expression cluster that consisting of BCL6, which cluster was cut at the height of 1.0.
One gene of each metric was counted as one point. Calculated the total score of each gene assessed.

miRNA-Gene Network Prediction
The pipeline was composed of three major steps: 1. miRNAs targeting each gene were obtained from the miRDB database, which is a miRNA target prediction database (44,45). 2. Counted the number of genes of each miRNA listed in (i), and the miRNA which targets more than one gene was selected. 3. Mapping according to the miRNA-gene pair list acquired in (ii). In this step, the miRNAs were pooled for miRNA families if there were more than one miRNA belonging to the same miRNA family. 4. Ranked the miRNAs or the miRNA family according to the number of genes targeted.

Plasma Preparation
Whole blood was collected from controls and TAK patients in ethylenediaminetetraacetic acid (EDTA) tubes and stored at 4°C. Plasma was isolated immediately from the whole blood by centrifugation (4°C, 800 × g for 10 min) in 4 hours. Then, the collected plasma was isolated by further centrifugation (4°C, 3,000 × g for 15 min) to remove cellular debris and large vesicles. Plasma samples were stored at −80°C until analysis.

Exosome Isolation and Quantification for miRNA Sequencing
Plasma samples were thawed at 37°C and filtered through a filtration membrane (0.8 mm). Next, the plasma was diluted 1:1.5 in filtered 0.01 M phosphate-buffered saline (PBS). Exosomes were purified by consecutive steps of size exclusion chromatography (SEC) on size exclusion columns and ultrafiltration tubes (100 kDa cutoff, Amicon). Successful exosome isolation was confirmed as follows.
RNA was isolated using miRNeasy Mini Kit (Qiagen) from exosomes. All the RNA samples met the following RNA quality threshold.

miRNA Sequencing
Small RNA libraries were prepared with the QIAseq miRNA Library Kit (QIAGEN). The quality of the libraries was validated on an Agilent Bioanalyzer 2100 and qPCR. Pair-end sequencing was performed on Illumina HiSeq2500.

Statistical Analysis
Normality was assessed with a Kolmogorov-Smirnov. Normally distributed continuous variables were provided as mean ± standard deviation and non-normally distributed continuous variables as median (interquartile). A Chi-square test was used for reporting associations between two categorical variables. Differences of continuous variables between groups were analyzed by the Mann-Whitney test. The correlation between gene expression levels was represented by the Pearson correlation coefficient. The models were otherwise validated by examining standardized residuals for normal distribution. Statistical analysis was performed using IBM SPSS statistic V.23 (Armonk, New York, USA). A p-value of less than 0.05 was considered significant.

Demographic Data, Clinical Features, and Laboratory Findings of Patients
Twenty-seven TAK patients were included. Based on the definition described above, 7 patients were classified into the untreated TAK group, and 20 patients were classified into the treated TAK group. Among them, 15 patients were assessed for remission, and 12 were assessed as being at the active stage following the NIH criteria described previously. Besides, 10 healthy subjects served as controls.
The demographic data and laboratory findings of patients were listed in Supplementary Table 2, and the individual patient data were seen in Table 1. There was no significant difference between the inactive-treated TAK group and the active-treated TAK group in age (p=0.82), sex (p=0.35), disease duration (p=0.1), erythrocyte sedimentation rate (ESR) (p=0.33), IL-6 (p=0.1), tumor necrosis factor (TNF)-a (p=0.66), and the dose of corticosteroid (p=0.37). The active-treated TAK group had a higher hypersensitive-C reactive protein (hs-CRP) level than the inactive-treated TAK group (p=0.02).

Functional Enrichment Analysis and Interaction Analysis of Selected Genes
The function of the selected genes according to UniProt was summarized in Supplementary Table 3. To gain further insight into the function of the 29 selected genes, we performed a functional enrichment analysis based on GO, KEGG, and WikiPathways databases. The functional enrichment analysis results confirmed that the 29 selected genes were closely related to this research topic. Figure 2A showed these genes were highly enriched in TLR signaling pathways and the function of the activation and differentiation of T-cells. Figure 2B showed the one-to-one correspondence between the genes and some GO terms and KEGG pathways which were closely related to this research topic, including positive/negative regulation of T-cell activation, I-kappaB kinase/NF-kappaB signaling, and Toll-like receptor signaling pathway., and other functions or pathways that are crucial to the pathogenesis of AIDs, such as B-cell mediated immunity, regulation of interleukin production, external side of plasma membrane, NAD+ nucleosidase activity, DNA-binding transcription repressor activity, RNA polymerase II-specific, NOD-like receptor signaling pathway, RIG-I-like receptor signaling pathway, PI3K-Akt signaling pathway, and MAPK signaling pathway. The detailed enrichment results were seen in Supplementary Table 4.
Next, to characterize possible molecular interactions across these selected genes, we built a protein-protein interaction (PPI) map using the STRING database. Figure 2C showed the PPI network clustering-based K-means (k = 3). The 3 clusters were each labeled with a different color, in which the proteins were functionally related. The genes labeled as red coded the first signaling molecule of T-cell activation or co-stimulatory molecules except FOXP3, the genes labeled as blue coded key subset transcription factors, and most genes labeled as green belonged to TLR signaling pathway. Notably, NR4A1 was in the same cluster as RELA, NFKB1 and, NFKBIA. In addition to PPI networks, we also identified gene interaction using gene coexpression network (16). Figure 2D showed the gene coexpression network based on Coexpedia database, indicating the strong functional association among these genes. And these interactions would align with our experimental results later.

Sample Clustering Analysis and Gene Expression Differential Analysis in TLR Signaling Pathway
The heat map displayed the expression levels of the selected genes in individual samples, with red color indicating greater expression level, which showed a distinct pattern of gene expression in TAK patients compared with healthy controls, in untreated patients compared with treated patients, and in inactive patients compared with active patients ( Figure 3A). Next, to identify the power of the selected genes to distinguish different sample groups, we performed the unsupervised hierarchical clustering of different samples based on the gene expression level calculated by Pearson correlation coefficients. A sample-to-sample matrix revealed that 80% (8/10) of healthy controls could be separated from the TAK patients by hierarchical cluster analysis, which suggested that TAK patients had a different expression pattern of the selected genes from the healthy controls ( Figure 3B). Among the selected genes, TLR1, TLR2, TLR4, TLR6, TLR8, CD83, IkBa, p50, p65, TNF, and CCL5 were key genes in TLR signaling pathway. As gene expression differential analysis showed, compared with the HC, the untreated TAK patients had higher mRNA levels of TLR2 (p=0.043), TLR4 (p=0.014), IkBa (p=0.00021), p50 (p=0.00046), and p65 (p=0.00031), and the treated TAK patients had higher mRNA levels of TLR2 (p=0.015), TLR4 (p=0.044), and IkBa (p=0.00016). Compared with the untreated TAK patients, the treated TAK patients had lower The gene co-expression analysis based on the COEXPEDIA database. Functional enrichment analysis indicated these genes were closely related to TLR signaling pathways and the function of the activation and differentiation of T-cells. PPI analysis and gene co-expression analysis indicated that there were strong associations among these genes, and experimental co-expression relationships would be compared with these database-interactions for newly discovered coexpression relationships.
mRNA levels of p50 (p=0.013), and p65 (p=0.021). Compared the active-treated TAK patients with the inactive-treated TAK patients, p50 mRNA level (p=0.038) was higher and no increase in the mRNA level of TLRs was observed. ( Figure 3C). Although there was evidence supporting the upregulation of TLR signaling pathway in TAK, no substantial differences in TLR signaling pathway were found between the inactive-treated TAK patients and the active-treated TAK patients.

TLRs and Their Correlation Cluster Distinguish Active Patients From Inactive Patients in TAK
To more rigorously evaluate the functional relations among these genes, we calculated the pairwise Pearson correlation and the pairwise Spearman correlation between each pair of genes.
Although Pearson correlation is better suited for the establishment regression equations, we adopted Spearman correlation as a complementary analysis as Spearman correlation has some advantages in reflecting the functional associations between genes (39). Due to space limitations, most of the results based on Spearman correlation were provided in the supplementary materials, and Pearson correlation coefficient was presented in the text and adopted in the subsequent mathematical modeling. The most striking results to emerge from the data was that the inactive-treated TAK group exhibited a unique pattern of inverse correlations between the TLRs gene cluster (including TLR1/2/4/6/8, BCL6, TIGIT, NR4A1, IkBa, p50, TNF, CD83, PD-1, PD-L1, and TIM3), and the gene cluster associated with T-cell activation and differentiation (including TCR, CD28, T-bet, GATA3, FOXP3, CCL5, CD3, CD40L, CTLA4, PD-L2), either performed by Pearson correlation ( Figure 4A) or Spearman correlation analysis (Supplementary Figure 1A) analysis, while the active-treated TAK group did not ( Figure 4B and Supplementary Figure 1B). To more clearly demonstrate this feature, Supplementary Figures 1C-F only showed the high correlations (defined as |r| > 0.73, p < 0.01). CCL5 is a potent chemoattractant for blood monocytes, memory T-helper cells, and eosinophils, which is important in recruiting T-cells into inflammatory sites, and also activates the apoptotic cell death pathway in T cells (46). Among these inverse correlations, CCL5 was negatively co-expressed with TLR1 (Pearson's r= -0.675, p= 0.046), TLR2, TLR4 (Pearson's r= -0.878, p= 0.002), TLR6 (Pearson's r= -0.903, p= 0.001), and TLR8 (Pearson's r= -0.818, p= 0.007), and T-bet was negatively co-expressed with TLR4 (Pearson's r= -0.713, p= 0.031), TLR6 (Pearson's r= -0.755, p= 0.019), and TLR8 (Pearson's r= -0.837, showed that the inactive-treated TAK group exhibited a unique pattern of inverse correlations between the TLRs gene cluster and the gene cluster associated with T-cell activation and differentiation, while the active-treated TAK group did not. (C) The scatter plots, the linear regression, and the receiver operating characteristic curve (ROC) of the TLR4-CCL5 pair, the TLR6-CCL5 pair, the TLR8-CCL5 pair, and the TLR8-Tbet pair. (D, E) The ROC analysis was used to evaluate the assessment accuracy for the disease activity of TAK. The total score of each treated-TAK patient (D) and the ROC (E) when assessing the disease activity using the summed scoring system consisting of the above 4 gene pairs. TAK, Takayasu's arteritis. p= 0.005). We took the TLR4-CCL5 pair, the TLR6-CCL5 pair, the TLR8-CCL5 pair, and the TLR8-T-bet pair to establish linear regression models to assess the disease activity, and the values of AUC/sensitivity/specificity were 0.848/72.7%/100%, 0.818/ 63.6%/100%, 0.828/90%/80%, and 0.778/54.5%/100%, respectively ( Figure 4C). We then calculated the total score from the above 4 gene-pair models for each patient to assess the disease activity as mentioned in Method, the AUC, the sensitivity, the specificity, the positive predictive value, and the negative predictive value were 0.939, 90.9%, 88.9%, 90.9%, and 88.9%, respectively, which had greater diagnostic accuracy than did the single models ( Figures 4D, E). The assessment of individual patients was shown in Supplementary Table 5. The above results suggest that TLRs have a potential relationship with the disease activity of TAK, which requires further explanation.

Topological Analysis of Dynamic Gene Co-Expression Network Identifies the Functional Association Between TLRs and T-Cell Activation in TAK
The More Active Communication Between TLRs and T-Cell Activation in TAK Compared to Healthy Control Gene co-expression network analysis is a useful method to link tightly co-expressed gene modules to phenotypic traits (47). To determine gene modules associated with distinct disease stages of TAK, we constructed a dynamic gene co-expression network for each group based on the pairwise Pearson correlations or Spearman correlations for all genes. (Figure 5 and Supplementary Figure 2). Unsupervised hierarchical clustering was used to assess the function of genes. As a result, genes were organized into 5 clusters by cutting the clustering tree at the height of 1.0, which was indicated by red frames. Notably, genes belonging to the same cluster were like to have similar functions and were labeled with the same color. In active-treated TAK, the 5 clusters were as follows. In inactive-treated TAK, the 5 clusters were as follows. Comparing the two results, it could be seen that the 29 genes were divided into two broad categories: (i) TLR signaling pathway and (ii) the activation and differentiation of T-cells.
But compared with the inactive group, CD83, TNF, IkBa, p50, p65 as well as the co-stimulatory molecules (including PD-L1, PD-L2, and LAG3) were more closely related to TLRs in the active group. Besides, to present the main framework of the network clearer, another dynamic network with the observed value of the Pearson correlation was 0.01 was shown in the right column of Figure 5 and Supplementary Figure 2. This coexpression clustering revealed a functional association among the genes, providing insight into gene functions and networks.
To assess the functional communication among these genes at distinct stages, we calculated a number of topological network parameters commonly used to describe the network. Each of these networks had a short characteristic path length and a large clustering coefficient, suggesting that they participate in the biological processes that might be functionally related. Additionally, compared with the healthy controls, the activetreated TAK and the inactive-treated TAK showed shorter characteristic path length (healthy controls vs. active-treated TAK vs. inactive-treated TAK, Pearson correlation 3.265 vs.  Table 2. TLR-Co-Expression Signature: The Regression Equation Relating the TLR6 mRNA Level to the TLR4 mRNA Level Serves as a Biomarker of Active Disease in Treated TAK There was a tight interplay among TLRs. Figure 6A was the PPI network based on the STRING database analysis. TLRs function as a homodimer or heterodimer, such as TLR1/TLR2, TLR2/ TLR6, and TLR4/TLR6. Besides, there is some cross-talk between TLRs, for example, TLR7 and TLR9 ( Figure 6B). Interestingly, the co-expressed-TLR-pairs in different groups were different ( Figure 6C). The inactive-treated TAK group had high coexpression of TLR1 and TLR4 (r=0.804, p=0.009) which was absent in the active-treated TAK group, the untreated TAK group, or the HC group. The high co-expression of TLR4 and TLR6 (r=0.892, p=0.001) exist in the HC group (r=0.847, p=0.002) as well as the inactive-treated TAK group (r=0.892, p=0.001), while was absent in the active-treated TAK group. The different co-expressions might mean different functional relationships between TLRs at different stages.
Studies have found that a heterodimer of TLR4 and TLR6 promote a protracted sterile inflammatory response after being triggered by oxidized low-density lipoprotein (LDL) and bamyloid, which involves the pathogenesis of atherosclerosis and Alzheimer's disease (48)(49)(50). We found that the regression equation relating the TLR6 mRNA level to the TLR4 mRNA level might be a biomarker of active disease in treated TAK. Figure 6D showed the scatter plots illustrating the difference in the coexpression of TLR4 and TLR6 between the inactive-treated TAK and the active-treated TAK. The points of the active-treated TAK group, except the point presented Patient No. 12, were scattered around the fitted line of the scatter points of the inactive-treated TAK group. The regression model was established and the threshold of M-value for assessing the disease activity was set to 0.173 according to the Youden index. Patients with an Mvalue of more than 0.173 got one score and were assessed as being in the active stage. And the AUC, the sensitivity, specificity, positive predictive value, and negative predictive value were 0.919, 100%, 90.9%, 100%, and 90%, respectively ( Figures 6E, F). The details of the assessment were shown in Supplementary Table 6. The results indicated that the inactive stage of TAK could be characterized from the treated TAK by the co-expression of several TLR genes.

Genes Key to the Cross-Talk Between TLRs and the Activation and Differentiation of T-Cell in TAK
To identify genes closely related to TLRs in these genes key to the activation and differentiation of T-cells, we assessed the degree of functional association between the other genes were to TLRs, and the evaluation protocol is described in Method. Specific scores are shown in Table 3 and the results are detailed below. In inactive-treated TAK, most of the genes exhibited functional association with TLRs. The 7 genes with the highest scores were BCL6 (12), CCL5 (11), FOXP3 (7), GATA3 (7), CD28 (6), T-bet (6), and NR4A1 (6) according to Pearson correlation analysis, and the 5 genes with the highest scores were BCL6 (16), TIGIT (14), IkBa (13), NR4A1 (12), and FOXP3 (11) according to Spearman correlation analysis. However, in the active-treated group, the functional association between these genes and TLRs did not seem to be as strong as the inactive-treated group ( Table  3). The 3 genes with the highest scores were BCL6 (7), PD-1 (1), and LAG3 (1) according to Pearson correlation analysis, and the 5 genes with the highest scores were BCL6 (3), CD40 (3), and LAG3 (3) according to Spearman correlation analysis.
NR4A1 is a key transcription factor that drives T cell dysfunction and plays an important role in the apoptosis of self-reactive T cells (53,54). The results suggested that NR4A1 is likely to be functionally related to TLRs in TAK. Compared to the healthy controls, the untreated TAK group (p=0.0001) had an increased mRNA level of NR4A1, and that compared to the untreated TAK group, the treated TAK group (p=0.000005) had a decreased mRNA level of NR4A1 (Supplementary Figure 3). The details of the assessment were provided in Supplementary Table 7.

Different Signaling Pathways at Distinct Stages in TAK
To explore the possible mechanism and signaling pathways, we summarized the co-expression variations across different conditions, which reflected changes in the activated signaling pathways ( Figure 9A). Compared with the HC group, the untreated TAK group had 35 gene co-expression relations uniquely and lost 39 gene co-expression relations. Compared with the untreated TAK group, the FIGURE 5 | Dynamic gene co-expression networks based on Pearson correlation. Left panel, the gene co-expression networks consisting of correlation with a pvalue less than 0.05. Right panel, to more clearly demonstrate this feature, only the high correlations (defined as |r| > 0.73, p < 0.01) were shown. Middle panel, the hierarchical clustering trees using complete method. As a result, genes were organized into 5 clusters by cutting the clustering tree at the height of 1.0, which was indicated by red frames. Genes belonging to the same cluster were like to have similar functions and were labeled with the same color. The negative correlation was indicated by the yellow edge, while the positive correlation was gray in the networks. treated TAK group had 80 gene co-expression relations uniquely and lost 27 gene co-expression relations. Compared with the inactivetreated TAK group, the active-treated TAK group had 47 gene coexpression relations uniquely and lost 54 gene co-expression relations. As the treated TAK group, the inactive-treated TAK group, and the active-treated TAK group had so many gene co-expression relations that the relations with a p-value less than 0.01 were listed only, and the complete list was shown in Supplementary Table 8. We also built a Venn diagram to visualize overlapping co-expressions among the five conditions ( Figure 9B). Besides, the newly discovered co-expression relationships which had never been reported in STRING or Coexpedia database were indicated red in Figure 9A. The results showed the differences in signaling pathways at distinct stages.
Next, to classify the genes based on function, the heatmaps of hierarchical clustering based on the Pearson correlation were conducted ( Figure 9C). The heatmap of the treated TAK group was shown in Supplementary Figure 4. Notably, the five heatmaps of five groups had different clustering structures. The two largest clusters of each group for the subsequent analyses are detailed below.    Additionally, we performed GO enrichment analysis to further explore the signaling pathways associated with the abovementioned clusters using Metascape database, and visualized the enrichment results using ClueGO to interrogate functionally grouped gene ontology. The enrichment results are detailed below. As the results showed, in the active-treated TAK group, genes in Cluster 1 were significantly enriched in the signaling pathways related to the activation and differentiation of T-cells, while genes in Cluster 2 were significantly enriched for the regulation of cytokine production and response to bacterium and lipopolysaccharide ( Figure 10C). In the inactive-treated TAK group, genes in Cluster 1 were enriched for the regulation of cytokine production, the response to bacterium, peptide, and lipopolysaccharide, the TLR signaling pathway, the I-kappaB kinase/NF-kappaB signaling pathway, and the regulation of defense response, while genes in Cluster 2 were enriched for the regulation of T-cell activation and the regulation of IL-10 production ( Figure 10D).

miRNAs Might Play an Important Role in the Cross-Talk Between TLR and T-Cell in TAK Patients
The co-expression in Cluster 1 of the inactive-treated TAK group was not presently understood, while Cluster 1 of the active-treated TAK group was led by the activation and differentiation of T-cells. We surmised that the miRNA network might take part in the expression of these genes, which could account for the coexpression of Cluster 1 of inactive-treated TAK group. We predicted the miRNA that targeted TLR1, TLR2, TLR4, TLR6, TLR8, BCL6, NR4A1, NFKBIA, LAG3, HAVCR2, TIGIT, PDCD1, and CD40 separately using the miRDB database (44,45). Except for LAG3 among these genes, there were multiple miRNAs for each gene. We summarized the miRNA (rather than miRNA family) that might regulate two or more genes (Supplementary Table 9) and visualized the miRNA-gene regulatory network (Supplementary  Next, to test whether the miRNAs in the miRNA-gene network were differentially expressed, we sequenced miRNAs from the plasma exosomes from healthy control and TAK patients. We found that compared with the healthy controls, miR-548ad-5p showed a 25.6-fold upregulation (p=0.0012), miR-3613-5p showed a 2.35-fold upregulation (p=0.0012), and miR-335-5p showed a 2.07-fold upregulation (p=0.0039); while miR-335-3p showed a 3.23-fold downregulation (p=0.042) and miR-584-5p showed a 2.01-fold downregulation (p=0.0092). The network shown in Figure 11 only consisted of the miRNA that might regulate three or more genes and two differentially expressed miRNAs (miR-335-2p and miR-584-3p). In the network, miR-548 was the node with the highest degree (7 genes), followed by miR-5692 (4 genes), miR-4763 (4 genes), and miR-520 (4 genes). Figure 11 showed the results. The miR-548 family was associated with 7 genes (7/13) in Cluster 1 of the inactive-treated TAK group, including TLR1, TLR2, TLR4, TLR6, TLR8, BCL6, NFKBIA, NR4A1, and TIGIT, which suggested that the miR-548 family plays an important role in the co-expression of Cluster 1 of inactive-treated TAK group. miRNA sequencing identified 29 differentially expressed miRNAs, 17 of which were increased and 12 were decreased. To validate whether the TLR signaling pathway might be regulated by these differentially expressed miRNAs, we performed miRNA sequencing analysis following the workflow in Figure 12A. First, to identify whether the differentially expressed miRNAs could be associated with specific functional categories, we performed an unsupervised hierarchical clustering of miRNA expression level based on Spearman correlation coefficient and conducted GO enrichment analysis for each cluster based on DIANA-TarBase, an experimentally validated database (55). Figure 12B showed miRNAs were partitioned into three clusters by cutting the clustering tree at the height of 1.1, and miRNAs belonging to the same cluster might have a close functional association. Figure 12C showed the hierarchical clustering heat map. The functional enrichment analysis was described as follows.
The experimentally validated database analysis indicated that TLR signaling pathway might be one of the major targets of the differentially expressed miRNA-mediated regulation, especially of the miRNAs belonging to Cluster 1.
Next, to validate whether one differentially expressed miRNA could target multiple selected genes, which might be involved in the gene co-expression of the inactive-treated TAK group, we constructed the miRNA-Gene-network based on the interactions of miRNAs and genes in the miRTarBase database, another experimentally validated database (56), and the screening condition of "support type" was set to "Functional MTI (miRNA target-interactions)". Figure 12E showed the result. Within the network, we identified 7 genes (including TLR1, TLR2, TLR4, BCL6, NFKBIA, NR4A1, and RORC) that had been validated to be regulated by the same miRNA, miR-335-5p, 4 genes (including TLR4, BCL6, FOXP3, and NFKB1) by miR-21-5p, and 4 genes (including CD40, CD40LG, NFKB1, and TNF) by miR-34a-5p. The results demonstrated that as follows: 1. In TAK, differentially expressed miRNAs that targeted multiple selected genes do exist. 2. TLRs, BCL6, and FOXP3 might be regulated by common miRNAs in TAK.
To sum up, these results suggested that miRNAs might play an important role in the cross-talk between TLR and T-cell in TAK patients.

DISCUSSION
There have been a number of studies demonstrating that TLRs play an important role in the pathogenesis of many AIDs, such as RA, SLE, and MS, but in TAK, it is currently unclear whether TLRs are associated with the disease activity (1, 2). T-cell is one of the major driving forces for the inflammatory response in TAK progression (6). It has been currently unknown whether TLRs are related to the disease activity or the activation and differentiation of T-cells in TAK. In this work, we selected 29 genes that were functionally enriched for the TLR signaling pathway and the activation and differentiation of T-cells. Twenty-seven TAK patients were enrolled which were grouped into the untreated and the treated group (both were further separated into the inactive and the active group), and 10 adult healthy controls were included. The relative mRNA level data of PBMCs were acquired by RT-qPCR. First, differential gene expression analysis showed that the untreated TAK and the treated TAK had an increased mRNA level of TLR2 and TLR4 compared to healthy controls. A sample-to-sample matrix revealed that 80% of healthy controls could be separated from the TAK patients. These findings suggested that TAK patients had a different expression pattern of the selected genes from the healthy controls. Second, we identified the association between TLRs and the disease activity, as the linear regression models consisting of the TLR4-CCL5 pair, the TLR6-CCL5 pair, the TLR8-CCL5 pair, and the TLR8-T-bet pair could distinguish between active and inactive disease in TAK [the AUC/ sensitivity/specificity, 0.939/90.9%/88.9%]. Third, we identified the The results indicated that TLR signaling pathway might be one of the major targets of the differentially expressed miRNA-mediated regulation, especially of the miRNAs belonging to Cluster 1. (E) Differentially expressed miRNAs. The miRNA-Gene-network based on miRTarBase database, an experimentally validated database, and the screening condition of "support type" was set to "Functional MTI (miRNA target-interactions). Red and green colors of miRNA indicated upregulation and downregulation, respectively. The results demonstrated that differentially expressed miRNAs that targeted multiple selected genes do exist in TAK, and TLRs, BCL6, and FOXP3 might be regulated by common miRNAs in TAK. association between TLRs and the activation and differentiation of T-cells in TAK according to the following evidence: (1) As the dynamic gene co-expression network showed, compared with the healthy control group, the active-treated TAK group, and the inactive-treated TAK group had higher network connectivity, shorter characteristic path length, and a larger clustering coefficient, indicating the more active communication between TLRs and T-cell activation in TAK. (2) The inactive-treated TAK group exhibited a unique pattern of inverse correlations between the TLRs gene clusters (including TLR1/2/4/6/8, BCL6, TIGIT, NR4A1, and so on) and the gene cluster associated with T-cell activation and differentiation (including TCR, CD28, T-bet, GATA3, FOXP3, CCL5, and so on). Fourth, we explored the genes key to the cross-talk between TLRs and the activation and differentiation of T-cell in TAK. In inactive-treated TAK, BCL6, CCL5, FOXP3, GATA3, CD28, T-bet, TIGIT, IkBa, and NR4A1 were likely to have a close functional relation with TLRs. However, in the active-treated group, the association between these genes and TLRs did not seem to be as strong as the inactive-treated group, BCL6, PD-1, LAG3, and CD40 were likely to have a close functional relation with TLRs. Fifth, we analyzed the activated signaling pathways in the inactive-treated and the active-treated TAK group. Last, we predicted the miRNAs that involved the greatest cluster of the inactive-treated TAK group and validated that miRNA might play an important role in the cross-talk between TLR and T-cell in TAK by miRNA sequencing. Besides, we proposed a new concept of the TLR-co-expression signature which might distinguish different disease and treatment stages of TAK, such as the co-expression of TLR4 and TLR6, which serves as a biomarker of the active stage in treated TAK (AUC/sensitivity/specificity, 0.919/100%/90.9%). These findings were schematically summarized in Figure 13.

Elevated Levels of TLR2 and TLR4 in PBMCs From TAK Patients
An increased level of TLRs in PBMCs has been detected in some AIDs. For example, comparing SLE patients with healthy controls, TLR2 expression on monocytes was reduced, and intracellular TLR9 expression of CD19 + B-cells was elevated (57). Comparing Behcet's disease patients with healthy controls, TLR1 and TLR2 were elevated in B-cells, TLR1, 2, and 4 were more highly expressed in both CD4 + and CD8 + T-cells, granulocytes displayed a higher expression of TLR1, 2, 4, and 6, and the expressions of TLR2, 4, and 5 were significantly increased on classical monocytes (58). However, few studies probed the effect of the treatment on the expressions of the TLRs in AIDs. For TAK, it has been unclear whether there are abnormal expressions of the TLRs in TAK, except TLR4. Besides, it also has been unclear whether the abnormal expressions of the TLRs are related to the medication and the disease activity, including TLR4. Our results complement and add to findings from previous TAK studies in that the mRNA expression of TLR4 and its ligand S100s in PBMCs from TAK patients were higher compared to healthy controls (12). We found that both the untreated and the treated TAK patients had increased mRNA levels of TLR2 and TLR4, and there was no statistically significant difference between the two groups, which demonstrated that the increased mRNA levels of TLR2 and TLR4 were independent from the medication. Besides, it was unanticipated that there was no statistically significant difference between the inactive and the active TAK patients in the mRNA expression of TLR2 and TLR4, which suggested that the elevated levels of TLR2 and TLR4 might be not resulted from the disease activity in TAK. However, the treated patients in this study had a relatively low serum level of TNF-a, which might imply relatively inactive TLR-signaling pathways in their peripheral blood.

Novel Association Between TLRs and Disease Activity in AIDs
Many studies have shown that TLRs play a dual role in the disease activity of AIDs. On the one hand, TLRs promote disease progression. TLRs can be activated by the DAMPs which are released from the injured tissue, such as high mobility group box-1 (HGMB1) (59) and self-nucleic acids (60,61). On the other hand, TLRs although inhibit autoimmune inflammation. In SLE, TLRs can promote sterile inflammation activated by DAMPs released from damaged cells, which is closely related to the disease progression and autoantibody production (62). For example, TLR-7, -8, and -9, which are not only expressed in the innate immune cells, such as neutrophils, macrophages and dendritic cells (DCs), but also are expressed in the lymphocytes, such as T-cells, and B-cells (60,63,64), can be activated by the self-nucleic acids (DAMPs released from damaged cells) (60,61), and in recent years, there has been an increasing interest in the mechanism of how TLRs directly regulate the adaptive immune response without the innate immune cells, which deepens our understanding of the role of TLRs in the pathogenesis of AIDs. On the other hand, TLRs have a dual role in the pathogenesis of SLE, and TLRs can protect the body against autoimmune inflammation. Low expression of TLR9 due to single nucleotide polymorphisms would increase SLE susceptibility in humans (65), and deletion of TLR8 and TLR9 would accelerate lupus development in mice (66)(67)(68)(69)(70). A variety of mechanisms for TLR9 to promote and limit AIDs have been discovered. The TLR9-dependent function of medullary thymic epithelial cells is required for the generation of regulatory T cells (Tregs) and the establishment of central tolerance (71). TLR9 can inhibit the production of autoantibodies mediated by TLR7, such as anti-RNA Ab, anti-Sm Ab, anti-RNA Ab, anti-IgG2a RF (68,69,72). The balance between TLR7 and TLR9 is important for the formation of autoreactive B cells (73). And UNC93B1 plays an important role in regulating the balance between TLR7 and TLR9 (74)(75)(76)(77). In the endoplasmic reticulum, TLR7 competes with TLR9 to bind to UNC93B1, and the D34 position of UNC93B makes it biased to bind to TLR9 (78). In RA, TLRs also promote sterile inflammation. For example, TLR4 can be activated by high mobility group box-1 protein (HMGB1) in RA (59). Besides, the activation of TLRs caused by PAMPs also plays an important role in the pathogenesis of RA (79). In the clinic, some cytokines induced by TLRs are immunotherapeutic targets in AIDs, such as TNF-a, interleukin (IL)-6, interferon (IFN)-a, and IL-1b (7). Recently, NI-0101, an anti-toll-like receptor 4 monoclonal antibody was used in RA, which was the first clinical trial to target TLRs to treat autoimmune diseases (4). In this study, we observed a unique pattern of inverse correlations between the TLRs gene clusters and the gene cluster associated with T-cell activation and differentiation the inactive-treated TAK group, which suggests a novel relationship between TLRs and the disease activity in AIDs.
However, these potential mechanisms of these inverse correlations need to be elucidated in depth, and there are several fundamental questions to address. First, the inverse correlations between TLRs and the activation of T-cell should be described in more detailed, and some experimental data for the difference in the T-cell activation status are needed to be collected, including the expression of activation markers on T cells (such as CD69, CD25 and HLA-DR) and some cytokines (such as IFN-g + Th1, IL-4 + Th2, and IL-17 + Th17). Then, association studies cannot be used to infer causality. While there is a possible regulatory relationship between TLRs and Tcells in TAK, it is also possible that TLRs and T-cells are under the control of common upstream regulators. Furthermore, another question is why the inverse correlations exist in inactive TAK patients, but not in active TAK patients.

Activated Signaling Pathways in the Inactive and the Active TAK Patients
T-cell is one of the main forces for the autoimmune inflammation in TAK (6). It has been reported that multiple T-cell subtypes are related to the pathogenesis of TAK, including Th1, Th17, Th9, Tfh, and Th2-like Treg cells (80)(81)(82)(83). In this study, the functional enrichment of the greatest cluster of the active-treated TAK group was centered on the activation and differentiation of T-cells, indicating the activation of T-cells at the active stage, which is consistent with previously reported studies. Notably, the functional enrichment of the greatest cluster of the inactive-treated TAK group was centered on the regulation of cytokine production, the response to bacterium, peptide, and lipopolysaccharide, the TLR signaling pathway, the I-kappaB kinase/NF-kappaB signaling pathway, the regulation of defense response, and so on, which suggested the activated TLRsignaling pathways at the inactive stage.
The previous study has shown that BCL6 regulates nearly a third of the TLR4-regulated transcriptome, and that 90% of the BCL6 cistrome is collapsed following TLR4 activation in the macrophages (84). In B-cells, PELI1, which is activated by TLR3 and TLR4, directly interacted with BCL6, inducing lysine 63mediated BCL6 polyubiquitination, and PELI1 expression levels positively correlated with BCL6 expression (85). Our results showed that BCL6 had a very similar expression pattern to that of TLRs, indicating that there is a stronger and closer association between multiple TLRs and BCL6, such as direct interactions and regulatory relationships, and this association warrants further study. The transcription factor BCL6 is required for driving Tfh cell differentiation and regulates their function. Recent studies have demonstrated the implications of Tfh cells with the disease activity of TAK (82), and our results suggested that TLRs might play a role in the regulation of Tfh cells.

miRNA-548 Family
miRNAs play an important role in regulating gene expression and the TLR-signaling pathways (86,87). miRNAs with high sequence similarity form the miRNA family, co-regulating complex biological processes. Among them, the miR-548 family, with over 80 identified miRNAs, regulates the immune process in some diseases, such as cancer (88,89). A large amount studies have shown that miRNA-548 suppresses tumor proliferate by binding WNT2 (90), murine double minute 2 (91), metastasis tumor-associated protein-2 (92), specificity protein 1 (93), cancerous inhibitor of protein phosphatase 2A (94), HMGB1 (95), and so on. Inhibiting or down-regulating miR-548 promotes the tumor proliferation (96). One study proposed to serve miRNA-548 as a prognosis predictor in primary central nervous system lymphoma (97). Our study showed that miR-548 played an important role in the pathogenesis of autoimmune disease, which might be related to the TLR signaling pathway.
It is worth noting that miRNA-548 family was increased in the serum exosomes from TAK patients than healthy controls, while the expression levels of its potential target genes were not decreased. While we have no clear explanation for this, we speculated that miRNAs might play an important role in the negative feedback regulation of TLR signaling pathways based on previous studies. Studies have shown that miRNAs play an important role as negative feedback to control inflammatory responses maintaining a level, avoiding excess inflammatory responses (98). For example, serum IL-6 levels increase with aging, whereas miR-192 in extracellular vesicles alleviated the prolonged inflammation associated with aging (99). Besides, the expression level of miRNA-541 [can attenuate pro-inflammatory cytokine expression (100)] in circulating extracellular vesicles was negatively correlated with the pro-inflammatory cytokine expression levels and the number of adverse local symptoms after vaccination (101). In these cases, although the expression levels of some targeted genes are still increased than normal, miRNAs in circulating extracellular vesicles have functioned as avoiding excessive inflammation. However, we have not been demonstrated the correlations between miR-548 as the RNA sam p les a nd se ru m exos om es were not fr om the same participants.

TLR-Co-Expression Signature
There are certain functional relationships among TLRs. First, TLRs function as a homodimer or heterodimer. TLR2 binds TLR1 or TLR6 to recognize distinct PAMPs (102,103). The TLR4-TLR6-CD36 complex is activated by atherogenic lipids and amyloid-beta to stimulate sterile inflammation (48,49). Second, there is some cross-talk between TLRs. In SLE, TLR9 can inhibit the production of autoantibodies mediated by TLR7, such as anti-RNA Ab, anti-Sm Ab, anti-RNA Ab, anti-IgG2a RF (68,69,72). The balance between TLR7 and TLR9 is important for the formation of autoreactive B cells (73). And UNC93B1 plays an important role in regulating the balance between TLR7 and TLR9 (74)(75)(76)(77). In the endoplasmic reticulum, TLR7 competes with TLR9 to bind to UNC93B1, and the D34 position of UNC93B makes it biased to bind to TLR9 (78). In the endothelial cells transfected with TLR1in septic, the ability of TLR4 in these cells to respond to LPS was lost (104). TLR10 on peripheral blood monocytes reduces TLR2-induced cytokine production in Parkinson's disease (105), which is also a co-receptor of TLR2. Third, some TLR-genes locate the same gene cluster in the chromosome, and there are multiple associated SNPs of TLRs, such as TLR1, TLR6, and TLR10. For example, the observed multiple associated SNPs at the TLR6-TLR1-TLR10 gene cluster may play a role in prostate cancer risk (106). STRING database analysis showed strong PPIs among TLRs ( Figure 6A). We observed the different co-expressed TLRs clusters of different disease and treatment stages in TAK, which might serve as a signature of transcriptome analysis for individualized treatment decision. For instance, the co-expression TLR4 and TLR6 could distinguish the active-treated TAK patients from the inactive ones.

Limitations
In the correlation analysis results, there were many non-strong correlations (0.01<p<0.05) between the expression levels, which might be related to the sample heterogeneity. For instance, some non-strong correlation in the treated TAK group became strong correlations (p<or=0.01) after further stratifying patients into subgroups according to disease activity. But the sample size is relatively small, so the active or inactive group could not be separated into subgroups.
Besides, there were limitations in the protocol for evaluating the closeness of functional relationships between TLRs and the other genes. Although the genes with a high score were likely to be closely connected in function with TLRs, the genes with a low score were not necessarily intimately linked with TLRs.
Notably, the reference genes had a large influence on the correlation analysis results. Stably expressed reference genes help find out the co-expression relations. Moreover, different reference genes may cause some differences in the results, and combining Pearson correlation and Spearman correlation help partially address the problem.
Last, larger clinical studies will be needed to validate our findings and to calibrate the diagnostic thresholds. In this cohort, the naïve TAK patients will be also classified into the active group and the inactive group as it might exhibit a different coexpression network from the TAK patients under treatment. Some experimental data for the difference in the T cell activation status (activation markers and cytokines) will be also collected to construct more accurate and more detailed models of the association between TLRs and T-cell activation in TAK.

CONCLUSIONS
This study identified the association between TLRs and T-cell activation in TAK, found a potentially novel role of TLRs in the pathogenesis of autoimmune diseases, and highlighted the function of miRNAs in the cross-talk between TLRs and Tcells in TAK, and more investigation is needed to further confirm the role of miRNAs in the cross-talk between TLR and T-cell in TAK patients and to elucidate the mechanisms.

DATA AVAILABILITY STATEMENT
We have deposited the miRNA data in the OMIX, China National Center for Bioinformation/Beijing Institute of Genomics, Chinese Academy of Sciences (https://ngdc.cncb.ac. cn/omix: accession no. OMIX807).

ETHICS STATEMENT
The studies involving human participants were reviewed and approved by the Institutional Review Board of Peking Union Medical College Hospital. The patients/participants provided their written informed consent to participate in this study.

AUTHOR CONTRIBUTIONS
YT performed the RT-qPCR experiments and analyzed the data. BH performed the experiment of miRNA sequencing and differential miRNA expression analysis. JL, YT, and XT designed the study. JL and YT wrote the paper mainly. XT and XZ made the critical revision. All authors have read and approved the paper.