This analysis was conducted in accordance with a predetermined protocol, following the recommendations of a guideline for the reporting of systematic reviews of prognostic factor studies (Riley et al., 2019). The data collection and reporting were performed in accordance with the Preferred Reporting Items for Systematic Reviews and Meta-Analyses (PRISMA) (Moher et al., 2015) and the PRISMA Extension Statement for Reporting of Systemic Review Incorporating Network Meta-Analyses (Hutton et al., 2015).

We searched the PubMed, Web of Science, Medline, Embase, and Cochrane Library databases to identify studies comparing and evaluating the effects of various drug regimens used to regulate fetal tachycardias from database inception to July 2021, with no limits on the study design, country of origin, tachycardia type, administration route (placental, intravenous, intramuscular or other administration routes), or whether combination treatment was applied. Additional studies for inclusion were identified from among the references of relevant reviews (Alsaied et al., 2017; Hill et al., 2017; Carvalho, 2019; Yuan, 2020; Veduta et al., 2021). The search strategy was established previously, and the terms included digoxin; flecainide; sotalol; amiodarone; tachycardia(s); fetus; supraventricular tachycardia or SVT; and atrial flutter or AF. The PubMed search strategy was described as [Anti-Arrhythmia Agents (MeSH Terms)] OR [Amiodarone (MeSH Terms)] OR [Digoxin (MeSH Terms)] OR [Sotalol (MeSH Terms)] OR [Flecainide (MeSH Terms)] AND [Tachycardia (MeSH Terms)] OR [Tachycardia, Supraventricular (MeSH Terms)] OR [Atrial Flutter (MeSH Terms)] OR [Hydrops (MeSH Terms)] AND (fetus).

Two authors (Jiangwei Qin and Zhengrong Deng) independently performed the literature search and performed data extraction using a standardized, pre-established form that distinguished the tachycardia types or whether hydrops existed. The full texts of the identified studies were independently assessed for inclusion by two other authors (Yifei Li and Yimin Hua), based on predetermined criteria. Each identified study was assessed for the following inclusion criteria: 1) the study population included fetuses treated with medication for sustained SVT or AF, with a diagnosis of fetal tachycardia based on the results of echocardiography or cardiac MRI; 2) the interventions were defined as anti-arrhythmia drugs administrated to the fetuses using a transplacental approach, so comparisons were made among first-line regimens: digoxin monotherapy (D), digoxin and flecainide combination therapy (DF), flecainide monotherapy (F), digoxin and sotalol combination therapy (DS), and sotalol monotherapy (S) (some drugs such as verapamil and DFS combination were discarded because they were not sufficient to produce results). As for combination therapy, to ensure that the evaluation of the efficacy of the drug is put first in this meta-analysis, considering that many doctors would promptly add drugs when a single drug is ineffective, our definition of combination drugs would be described as: if a second drug was started before the third day of treatment, it was defined as combination therapy; 3) the cardioversion success rate was assessed as the primary outcome, which was defined as the reversion to a normal sinus rhythm during or after the administration of therapeutic drugs using a transplacental approach, without recurrence or relapse until birth; 4) and the use of an appropriate cohort study design. Cardioversions that were successfully controlled during the fetal stage were considered successful, even if tachycardia recurred postnatally. The intrauterine death ratio was as used as a secondary endpoint in the treatment arms, while postpartum deaths were not included. The exclusion criteria included 1) AET and PJRT, due to inadequate sample size; 2) failure to report the primary outcomes of interest; 3) the initial diagnosis occurring after 37 gestational weeks; 4) the study did not distinguish between specific types of fetal tachycardia among the enrolled cases.

Non-treatment and placebo treatments are clinically unfavorable for patients and are ethically disallowed for consistent fetal tachycardia; therefore, every regimen was used as the reference baseline for all network meta-analysis comparisons for the results of each set of two-by-two comparisons can be clearly described. In this analysis, we only evaluated the therapeutic efficacy of first-line transplacental medication administration on fetal tachycardia and assessed the associated safety of each treatment option. We also collected the publication year of the study, the dosage of treatment for existing CHD comorbidities, and the occurrence of both maternal and neonatal adverse events.

The risks of bias and article quality were assessed with the Newcastle–Ottawa Scale (NOS) (Wells et al., 2014) according to 3 aspects: Selection (4 stars), Comparability (2 stars), Outcome (3 stars). Publication bias was assessed using funnel plots in Stata, version 16.

Five therapeutic regimens were evaluated, including monotherapies and combination therapies, namely D monotherapy, DF combination therapy, F monotherapy, DS combination therapy, and S monotherapy, data of which were then fed into the model (see Supplementary Materials). The sotalol and flecainide combination therapy, amiodarone monotherapy, and verapamil monotherapy were excluded from the meta-analysis due to inadequate data availability. Five subgroups were established according to the type of tachycardia, including total (means all the included patients), SVT, AF, hydrops, and non-hydrops groups. We used a Bayesian hierarchical model to establish a network analysis framework and used Markov Chain Monte Carlo (MCMC) simulation procedures to generate posterior distributions by R package (gemtc) (Harrer et al., 2021).

In Bayesian statistics, we calculated credible intervals (Crl) rather than confidence intervals around our estimates, and a 95% credible interval indicating that ‘there is a 95% probability that the true value of the parameter falls within our interval’ rather than ‘of every confidence interval calculated with each data set, 95% of the credible intervals (Crl) contain the true values of the parameters’ in frequentist statistics. The Bayesian hierarchical model uses uninformative priors on the effect in every comparison which doesn’t have a big impact on the posterior results, and the effects of a multi-arm study stem came from a multivariate (normal) distribution (Harrer et al., 2021). The model considered the replaceability of each treatment and used collected direct and indirect data based on uninformative priors, in order to obtain a relative accurate posterior result. The MCMC sampling allows to estimate the posterior distributions of our parameters, and thus generate the results of network meta-analysis (Harrer et al., 2021). Due to the inability to control for variables during observational studies, we did not assume that each study is an estimator of the same true effect size, but there are “study-specific” true effects estimated by each observed effect size and they are part of an overarching distribution of true effect sizes, and the variation between studies was estimated using a random-effects model.

A network plot was generated for the primary endpoint in each group using a random-effects model, using four Markov chains, 5,000 burn-in iterations, and 20,000 simulation iterations, where the parameters were considered enough for sensitive and well-converging network models, details on parameters could be accessed by command: ?mtc.run of R package (gemtc) (Harrer et al., 2021). Convergence assessments were performed using Gelman–Rubin plots. This plot used the term Potential Scale Reduction Factor (PSRF) to compare the variation within each chain to the variation between chains, and how both develop over time. The PRSF would gradually shrink down with increasing numbers of iterations, and should at least be below 1.05 in the end for an enough converged model. We also evaluated the consistency of our network model using the node-split method. Using the network meta-analysis results, we calculated the probability of a treatment being the first, second, third, fourth, or fifth ranking treatment; presented the comparison in one cumulative probability plot; data of each set of two-by-two comparisons would be arranged in a table. We only represented forest plot by comparing with D monotherapy. The same approach was used to assess the secondary outcome of mortality across the whole population.

We integrated the safety and efficacy data in the total tachycardia population by introducing a metric we referred to as the term safety index, which was calculated as the number of deaths divided by the number of successful cardioversions, to simultaneously assess the risks and benefits for each therapeutic regimen. It was derived from a pharmacologic concept——therapeutic index, which evaluates the safety of a drug and calculated as the median lethal dose divided by median effective dose. We hypothesized that a higher safety index value indicates reduced benefits of a regimen, because there were more death events with the same recovery events. Different with the comparison barely on death rate and cardioversion rate, safety index may complement the analysis of the data. For example:

-In study A, 20 of D group for first-line treatment and 10 of F, 12 of D for cardioversion and 7 of F, and 1 of D for death and 1 of F.

From our point of view, the difference may have relation with existence of the unrecovered but also alive population after first-line therapy, causing some data hidden behind. So safety index does provide another aspect for explaining the results.

The R (gemtc) package (van Valkenhoef et al., 2012) and Stata were used to perform statistical calculations for this network meta-analysis (Salanti, 2012; Shim et al., 2017; Shim et al., 2019; Harrer et al., 2021). We have put the R code and collected data in Supplementary Materials for convenience of replication.

After the database search, 11,075 publications were identified, and 26 observational studies satisfied the established inclusion and exclusion criteria after screening and eligibility procedures were applied (shown in Figure 1). Any discrepancies were resolved between two reviewers by discussion until consensus was reached. Among the 26 included studies, one was a prospective study, and 25 were retrospective studies (van Engelen et al., 1994; Frohn-Mulder et al., 1995; Naumburg et al., 1997; Lisowski et al., 2000; Oudijk et al., 2000; Ebenroth et al., 2001; Jouannic et al., 2002; Krapp et al., 2002; Boldt et al., 2003; Jouannic et al., 2003; Oudijk et al., 2003; D'Alto et al., 2008; Pézard et al., 2008; Lulic Jurjevic et al., 2009; Hahurij et al., 2011; Shah et al., 2012; Uzun et al., 2012; van der Heijden et al., 2013; Ekman-Joelsson et al., 2015; Sridharan et al., 2016; Strizek et al., 2016; Ekiz et al., 2018; Karmegeraj et al., 2018; Miyoshi et al., 2019; O'Leary et al., 2020; Tunca Sahin et al., 2021). Only one randomized controlled trial (RCT; www.fasttherapytrial.com) (Edgar Jaeggi, MD, etc.) was retrieved during this search, but this RCT has not finished; therefore, it was treated as gray literature. Detailed reasons for the exclusion of identified articles are provided in Figure 1. The characteristics of the 26 included studies were summarized in Supplementary Table S1. 21 studies were conducted in Europe with 7 in Netherlands, and 3 studies were established in United States, 1 in Japan, 1 in India. Included gestational age at birth or when outcome recorded was ranged from 23 to 42 weeks but concentrated in 38 weeks. Quantitative data were organized into five groups according to treatment (D, DF, F, DS, and S). The maternal, neonatal and follow-up adverse events data were extracted regardless of the tachycardia type and specific regimens, because this information was not always reported in detail. The data set of extracted results that were used for statistical analyses are presented in the Supplementary Materials.

For the primary outcomes, the data volumes and model stability were assessed for each of the five groups. The network plots represented the data size and the comparison size between two regimens, and the edge thickness corresponded to the numbers of studies included for comparisons. Total, SVT, and hydrops groups were found to form a complete network structure, with relationships identified across all treatment regimens; by contrast, some direct and indirect comparisons were lacking for the AF and non-hydrops groups, including F vs. DS and DF vs. S, suggesting that more attention should be paid when interpreting the results for these two groups (Figures 2A–E). Gelman–Rubin plots indicated that the shrink factor was stably maintained below 1.05 for all five groups (Supplementary Figures S1–S5), indicating enough converged models being applied. To assess inconsistencies, the node-split forest plots (Supplementary Tables S2–S8) revealed the effects of various comparisons when using only direct, only indirect, or all available evidence. Inconsistency was identified in all five groups due to one or more comparisons with p < 0.05, which indicated that the studies included systematically different populations. As shown in Figures 3A–E, 4A–E, the top two best-ranked regimens were consistently identified as DF and F. The odds ratio (OR) and 95% credible interval (95% CrI) when compared with D as the reference therapy was significantly different across all groups for DF: Total, 2.44 (95% CrI: 1.59, 3.52); SVT, 2.77 (95% CrI: 1.59, 4.07); AF, 67.85 (95% CrI: 14.25, 168.68); hydrops, 6.03 (95% CrI: 2.54, 10.68); and non-hydrops, 5.06 (95% CrI: 1.87, 9.88). Other therapies showed significant effects for specific groups when compared with different regimens, details were shown in Table 1. The results indicated there wasn’t one regimen having all statistically significant differences when comparing to the other therapies. Relatively, DF and F administrations demonstrated better efficacy for controlling fetal heart rate, regardless of the tachycardia type or the presence of hydrops, indicating the superiority of the two treatment protocols. DF demonstrated superiority on D in all five subgroups. In total group and AF group, DF had better effectiveness than DS. In AF group, superiority of F regimen could also be observed. Majority of comparisons showed no significant differences.

To assess the secondary outcome of death rates, a network plot (Figure 2F), Gelman–Rubin plot (Supplementary Figure S6), and inconsistency test (Supplementary Table S7) were used to examine intrauterine death. The results demonstrated no significant differences, indicating no differences in intrauterine death outcomes across the examined treatment regimens (Figure 3F and Table 2). As shown in Figure 4F, the two treatments found to be the most efficacious, DF and F, were also the top-ranked regimens for the death rate assessment. Therefore, we attempted to further assess the benefits and risks of each treatment regimen by analyzing the safety index.

It was necessary to mention that only the data of cardioversion must be no zero to assess a safety index, studies with substandard data were removed, resulting in a loss of data volume. For example, data of Jouannic, J, M. et al. in Supplementary Materials as we provided online, the cardioversion for total were 0 in D and DS groups, so we had to remove both data for code running. Figure 2G and Supplementary Figure S7 show acceptable results, but the inconsistency assessment (Supplementary Table S8) indicated an increase in the number of inconsistent comparisons, reducing the credibility of the cumulative probability plot. Table 2 showed the absence of significant results, preventing the determination of absolute ranking among these five regimens (Figures 3G, 4G). To summarize, the DF and F regimens demonstrated relatively superior efficacy for terminating fetal tachycardias, regardless of the tachycardia type or the presence of hydrops, with no significant difference in the associated risks such as death rate.

Supplementary Table S9 indicated that a low risk of bias was determined for all 26 studies included in this network meta-analysis (all studies obtained stars ≥7 (9 stars in total)) The publication bias of primary endpoints in the five groups is shown in Supplementary Figure S8, and no evident bias was observed for the total, SVT, hydrops, or safety index subgroups.

The use of Bayesian model analysis can prevent the second-class errors that ordinary meta-analyses are prone to producing when assessing multiple comparisons. In addition, our model incorporates a larger samples size (799 patients in total group, subgroups data could be accessed in Supplementary Materials) than previous meta-analyses.

Our included samples were limited to cases of SVT or AF, and other rarer manifestations were not included. In addition, Inconsistency identified in all seven data set indicated that the studies included systematically different populations, which were likely due to various causes, including differences in the severity of tachycardia, the gestation period, drug doses, and administration routes other than placental administration, such as muscular injection, umbilical vein injection, or oral administration, which were difficult to unify or control across these 26 observational studies. It should be noted that most of the samples in this paper have normal heart structure, so the conclusions in this paper are more suitable for children with normal heart structure. Whether the included retrospective studies applied similar practical treatment conditions was difficult to determine, and retrospective observation data is unable to meet randomization requirements. Prospective studies and ongoing RCT studies would provide more methodological control. However, Efthimiou et al. (2017) proposed that the inclusion of evidence from non-random studies (observational studies) can improve the accuracy of the results obtained from RCT-based network meta-analyses; therefore, this article can provide supplemental data for any follow-up studies. Due to the complexity of the Bayesian model, which requires similar distributions by default, the results may be influenced by non-similar distributions.