We conducted a systematic review and network meta-analysis based on a registered protocol (PROSPERO ID CRD42021287837) and reported it according to the Preferred Reporting Items for Systematic Reviews and Meta-Analyses (PRISMA) guidelines (Supplementary Material S1). Bayesian network meta-analysis was used because it offers a more straightforward method in generating probabilistic statements and in predicting treatment effects.

We searched PubMed, Embase, Cochrane Central Register of Clinical Trials, Web of Science, China National Knowledge Infrastructure (CNKI), and Wan-fang (WF) data, and included reports until December 2021 in all languages. The main search terms used were “Typhus Scrub,” “Antibiotics,” and “Treatment” (Supplementary Material S2). The references of the included articles to articles that were missing from the first search.

All studies were evaluated and included in the meta-analysis according to the following criteria: (1) Treatment of scrub typhus with antibiotics, including RCTs or quasi-RCTs. (2) Patients were diagnosed with scrub typhus based on clinical symptoms or laboratory tests (7–9).

The exclusion criteria were as follows: (1) RCTs had overlapping patients. (2) The primary outcome was not described or the study reported ambiguous clinical outcomes. (3) The study design was incomplete and/or ambiguous. (4) Combined interventional studies. If several articles were published on the same population intervention, only the study that was complete and of higher quality was included. Two reviewers (DF and ZH) evaluated whether each of the selected RCTs met the intended criteria independently and consulted a third reviewer (TW) in case of a disagreement.

Data (e.g., first author, publication year, and patient characteristics), treatments, and reported outcomes were extracted. Extracted data were assessed by two independent authors (DF and ZH) to avoid potential assessment bias.

The risk of bias was assessed for each trial using the Cochrane Handbook for Systematic Reviews of Interventions. Items were considered as low, high, or uncertain risk of bias (10). We incorporated results into the Grading of Recommendations, Assessment, Development, and Evaluation (GRADE) to assess the credibility of findings from each network meta-analysis. The GRADE system evaluates the quality of evidence at four levels: high, moderate, low, and very low (11). Two investigators (DF and TW) independently assessed the risk of bias of individual studies and GRADE system. Any discrepancies were resolved by consensus and arbitration by the authors (DF, JH, TW, and ZH).

We synthesized evidence to compare different treatments in terms of efficacy and safety, and reported our findings as odds ratios (ORs) for binary outcomes [cure rate (CR) and adverse events (ADs)], and as mean differences (MDs) for continuous outcomes [defervescence time (DT)] along with corresponding 95% credible intervals (CrIs). The primary outcome was CR. Secondary outcomes were DT and ADs as reported by the studies.

Network meta-analyses were performed in a Bayesian framework using a Markov chain Monte Carlo simulation technique in gemtc package (R version 4.0.2). For each outcome measure, the fixed-effects and random-effects consistency models were used. Four Markov Chain Monte Carlo (MCMC) were established for running 50,000 iterations, 20,000 sample burn-ins, and a thinning interval of 1. The Brooks-Gelman-Rubin diagnostic was used to evaluate and visualize model convergence of the iterations. Once convergence was established, posterior distributions for model parameters were obtained as the output of the network meta-analysis estimate (MD/OR and the corresponding 95% CrIs). In the presence of minimally informative priors, credible intervals can be interpreted such as conventional confidence intervals (CIs). Network meta-analysis estimated the overall rankings of treatments by calculating the surface under the cumulative ranking curves (SUCRA) for each, which equals 1 when a treatment is the best and 0 when a treatment is the worst. Transitivity was evaluated using descriptive statistics for the study and population baselines, including sample size, age, and gender (12, 13).

We used nodal analysis to evaluate the inconsistency of this NMA by comparing differences between direct and indirect evidence. Heterogeneity between studies was assessed using the I² and P statistic within a visual forest plot. A p-value of 0.05 was considered statistically significant. Heterogeneity was considered low, moderate, or high for estimated I² under 25%, between 25 and 50%, and over 50%, respectively. The Egger regression test with a funnel plot was used to assess publication bias, and p < 0.05 was considered to indicate significant asymmetry and publication bias. Using the netmeta package in R.

Sensitivity analysis was used to assess the impact of individual or several specific studies on overall heterogeneity. We hypothesized that the inclusion of various study populations and conditions might contribute to heterogeneity and inconsistency. Thus, we assessed the sensitivity of our findings by repeating each network meta-analysis after excluding partial specificity studies (14).

Meta regressions were performed using controlled-placebo adjusted data to examine the relationship between treatment outcome–associated changes and baseline age and sex, publication year and quality of the study, and whether or not it was a randomized study (14, 15). If a multi-arm study was identified based on these meta regressions, the estimates for each group were combined (14). The metafor package in R was used for analysis.

The online database search yielded 4,411 results, and 4 were found from a secondary source. After the removal of duplicate entries, a total of 60 articles were further considered for full-text assessment. Finally, 16 studies enrolling a total of 1,582 patients were found to be eligible for further analysis (Figure 1) (7, 16–30). Patients were enrolled to receive any of the following 7 treatments: azithromycin, doxycycline, chloramphenicol, tetracycline, rifampin, moxifloxacin, or telithromycin. Comparative network plots are presented in Figure 2. The main characteristics of included studies are reported in Table 1. The main mechanisms of action of the included drugs are reported in Table 2. The baseline characteristics of studies included in the network meta-analysis have been detailed in Supplementary Materials S3–S5.

The network meta-analysis included all studies for CR (Figure 2A), 13 studies for DT (Figure 2B), 14 studies for risk of GS-AD (Figure 2C), and 7 studies for risk of ABC-AD (Figure 2D).

With respect to CR (Figure 3A), rifampin was found to have the highest benefit vs. moxifloxacin [OR 43.056, 95%CrI (3.01, 864.9)]. Benefits were also observed with the use of tetracycline [13.83 (1.488, 181.78)], chloramphenicol [18.616 (2.258, 198.05)], azithromycin [15.02 (2.258, 135.7)], and doxycycline [9.896 (2.306, 60.13)] vs. moxifloxacin.

With respect to ABC-AD (Figure 3D), azithromycin was found to have the highest benefit vs. doxycycline [0.27 (0.083, 0.844)], whereas with respect to DT (Figure 3B) and GS-AD (Figure 3C), none of the treatments showed significant differences among each other.

The Bayesian ranking results were almost in line with those from pooled analyses using ORs and MDs. The meta-analysis revealed that rifampin was most likely to be ranked first for CR (cumulative probability 82%) and DT (23%). Moxifloxacin was most likely to be ranked last for CR (3%) and DT (85%). Doxycycline was most likely to cause GS-AD (72%), followed by tetracycline (66%). Doxycycline had the highest |probability (87%) of being ranked last in causing ABC-AD, followed by chloramphenicol (65%) and azithromycin (2%). The probability profile of each class of drug is shown in Figure 4 and Supplementary Material S8.

The Brooks-Gelman-Rubin diagnostic method demonstrated better stability and replicability for each MCMC chain, with positive confidence in the results (Supplementary Material S12). Results of the evaluation of Inconsistency for all comparisons are presented in Supplementary Material S9. We used a significance level of P > 0.05 for all cases. This finding indicates positive transitivity and consistency among comparisons. Heterogeneity was assessed using the I² and P statistic statistics. Among the included trials, the heterogeneity in the CR, AD, and DT groups was determined to be as follows: CR: I² = 1.28%, ABC-AD: I² = 0%, GS-AD: I² = 60.48%, DT: I² = 95.7%. Egger regression tests were performed with P-values of CR: 0.14, ABC-AD: 0.63, GS-AD: 0.52, and DT: 0.41, as shown in Supplementary Materials S10, S11.

Overall, in combination with previous studies risk of bias assessments, these trials were considered moderate-high risk for bias. Supplementary Material S7 summarizes the detailed risk of bias assessments. The certainty of evidence was evaluated using the GRADE system (Supplementary Material S13). A small number were rated as “Moderate” (4/21 in CR and 3/3 in ABC-AD), the others were rated as “Low” or “Very Low”.

The sensitivity analysis for 3 outcomes (sensitivity analysis and regression analysis were not performed for the ABC-AD group because of the small number of studies) were assessed by repeating analyses after the sequential exclusion of the following studies: 5 studies that did not use randomization, 3 studies that were published prior to 2000, 2 studies of scrub typhus in adolescents, 3 studies that had fewer than 30 patients per arm, and 7 studies that included intravenous interventions. The inclusion of these studies did not have a significant impact on the results. The assessment of heterogeneity (I²) was similar across studies, except that I² was reduced noticeably in the paired meta-analysis and network meta-analysis when publications prior to 2000 or a group of fewer than 30 patients were excluded in the GS-AD group (Supplementary Material S14).

In meta regressions analysis, greater drug therapy–induced GS-ADs were associated with higher proportion of male rate [Coefficient = 0.42 percentage⁻¹ (0.004–0.081), z = 2.13, p = 0.033] and higher mean age of the study [Coefficient = 0.03 age⁻¹ (0.003–0.057), z = 2.17, p = 0.0301]. We did not find evidence that the incidence of GS-AD was related to study blinding, or paper publication year. The increase in DT was associated with study blinding [Coefficient = 8.71 (4.80–12.62), z = 4.368, p ≤ 0.001] and closer paper publication year [Coefficient = −0.34 year⁻¹ (−0.62 to −0.06), z = −2.35, p = 0.018]. We did not find evidence that DT was related to the proportion of male patients or the mean age of patients in the study. Moreover, we did not find evidence that CR was related to these baseline differences (Supplementary Material S15).