In the last decade well beyond half a million bacterial genomes have been sequenced worldwide, about 7% of these from Mycobacterium tuberculosis (MTB) strains (Blackwell et al., 2021). One of the reasons behind these extensive sequencing efforts is the use of whole genome sequencing in molecular epidemiology. Molecular epidemiology studies of MTB (and of other pathogens) use microbial genome sequences sampled from different patients to investigate epidemiological dynamics such as transmission, relapses, and the acquisition and spread of antibiotic resistance (Hatherell et al., 2016; Guthrie and Gardy, 2017; Nikolayevskyy et al., 2019). One of the most popular approaches to analyze this data is to cluster strains in groups based on their genetic distance. The identification of clustered MTB strains is commonly interpreted as evidence for recent local transmission, while patients infected with non-clustered strains (singletons) are thought to be novel introductions (i.e. patients that got infected somewhere else). Similarly, when studying the epidemiology of antibiotic resistance, clustered strains are considered as cases of transmission of resistance, while singletons are thought to be more likely to represent instances of resistance acquisition (Hatherell et al., 2016).

In MTB studies, clusters are often defined based on a single nucleotide polymorphism (SNP) threshold: strains with fewer SNPs than a given threshold are grouped together in the same cluster. However, this has been criticized, because it is not clear which value should be used (Stimson et al., 2019). A review of more than 30 publications concluded that a threshold of six SNPs could be used to identify cases of direct transmission (Nikolayevskyy et al., 2019), but other studies proposed different thresholds for different settings (see Table 1 in Hatherell et al., 2016 for some examples). Alternative approaches have been proposed to overcome the limits of SNP thresholds. For example, the method proposed by Stimson et al., 2019 groups strains based on a probabilistic model that takes into account the clock rate and the time of sampling. Other studies sidestepped the choice of one specific value by performing clustering multiple times with different thresholds, and comparing the results (Holt et al., 2018; Meehan et al., 2018; Yang et al., 2018; López et al., 2020; Shuaib et al., 2020; Cox et al., 2021, Liu et al., 2021; Walter et al., 2022; Yang et al., 2022). Despite their limitations, clustering methods are considered useful to study transmission dynamics in MTB. Many studies tested the association of clustered strains with host and bacterial sub-populations, such as different age groups, HIV-positive patients, bacterial lineages and others (Guerra-Assunção et al., 2015; Asare et al., 2020; Sobkowiak et al., 2020; Cox et al., 2021, Gygli et al., 2021; Merker et al., 2021; Yang et al., 2022). In these analyses, a positive association is interpreted as evidence for increased transmission of a certain sub-population. Further studies used clustering rates (the percentage of clustered strains), to characterize the extent of transmission in bacterial sub-populations (Holt et al., 2018; Shuaib et al., 2020). For example, Holt et al., 2018 found that in Vietnam, MTB Lineage 2 (L2) had higher clustering rates compared to MTB Lineage 4 (L4) and MTB Lineage 1 (L1), and interpreted these results as evidence for more frequent transmission of L2 strains, compared to L4 and L1 strains. Finally, in recent studies the distribution of terminal branch lengths (TBL) was used as a proxy for transmissibility in different MTB lineages (Holt et al., 2018; Freschi et al., 2021; Walter et al., 2022), partially complementing classical clustering analyses. All these approaches are based on the assumption that increased transmission results in increased clustering rates and shorter terminal branches. However, it is known that other factors could influence the results of clustering, for example it was posited that higher rates of molecular evolution, and low sampling rates should lead to lower clustering rates (Stimson et al., 2019, Menardo et al., 2019).

There is a consensus that epidemiological dynamics have an influence on the shape of MTB phylogenies, and therefore on clustering and TBL, though how they do so was never explored with quantitative studies. Here, I used simulations to explore what factors influence clustering results and TBL in MTB epidemics. I simulated the molecular evolution of MTB strains under different epidemiological conditions. I then inferred phylogenetic trees and computed clustering rates and TBL from the simulated data. With this approach I investigated the molecular clock rate, sampling proportion, sampling period, transmission rate, basic reproductive number (R0), length of the latency period, and length of the infectious period. I found that all these factors affect the results of clustering and TBL, except for the length of the infectious period. These results are in contradiction with the standard interpretation of MTB epidemiological studies, namely that sub-populations associated with clustering and shorter terminal branches are necessarily transmitting more.

To investigate the expected patterns of genetic diversity under different epidemiological conditions, I assembled a pipeline to simulate the evolution of MTB genomes in different epidemiological settings (Figure 1). The details are reported in the Materials and methods section. Briefly, the pipeline simulates a transmission tree using a birth-death model in which transmission events occur at rate λ and originate in the infectious compartment I, leading to a new exposed individual in compartment E. Exposed individuals become infectious at rate ψ, by moving from compartment E to compartment I. The time spent in the compartment E represents the latency period. Individuals are removed from compartment I by death or self-cure, occurring at rate σ, or sampling, occurring at rate ε. The sampling rate ε corresponds to the rate at which infectious individuals are diagnosed and treated. With the onset of treatment, it is assumed that patients stop being infectious immediately. The time spent in compartment I represents the infectious period.

Figure 1

Download asset Open asset

A similar epidemiological model was used previously in a phylodynamic analysis of two MTB outbreaks (Kühnert et al., 2016; Kühnert et al., 2018). In a second step, the pipeline simulates the evolution of genome sequences along the tree given a clock rate π, and it computes the clustering rates under different SNP thresholds, and the terminal branch lengths. This pipeline allows to test how clustering rates and TBL distributions change under different sampling proportions, sampling periods, molecular clock rates, transmission rates, basic reproductive numbers (R0 = λ/(σ+ε)), lengths of the latency period (determined by ψ), and lengths of the infectious period (determined by σ+ε). Moreover, with this approach it is possible to investigate how the different factors impact the choice of a sensible SNP threshold. To do this, I defined the ‘95% sensitivity SNP threshold’ as the minimum threshold for which at least 95% of strains are clustered in at least 95% of the simulations. A lower SNP threshold would lead to lower sensitivity (i.e. simulated samples, which are the result of recent transmission, would not be clustered), while a larger threshold would lead to low specificity, although this cannot be quantified with this analysis (see Materials and methods for additional information).

Latency

Next, I tested whether differences in the duration of the latent period could result in different clustering rates and TBL. Here, latency is defined as the period in which an individual is infected but not yet infectious. Typically, the shift to infectiousness in TB patients is considered to occur with the onset of symptoms. However, in recent years, the importance of sub-clinical TB has been reconsidered. It is possible that a considerable part of TB transmission occurs from asymptomatic patients, although this has not been yet quantified (Kendall et al., 2021). Given the uncertainty about the length of the latent period I tested three different rates of progression to infectiousness ψ=0.5, 1, 2, corresponding to a median latent period of ~16.6, 8.3, and 4.2 months, respectively. These values represent the range of duration of asymptomatic infection estimated in different countries (Ku et al., 2021). All other parameters were constant in all simulations (π=8 × 10^–8; σ = ε=0.5; λ=1, sampling period = 10 years). I found that longer latency resulted in lower clustering rates and longer TBL (Figure 2, Table 1). Correspondingly, the 95% sensitivity threshold was 10, 6, and 5 SNPs, and the mean of the TBL distribution was 0.87, 0.56, and 0.41, respectively for long, mid, and short latency.

Figure 2

Download asset Open asset

Table 1

Scenario	λ	ε	R0	σ	ψ	π	95% SNP threshold	100% SNP threshold	95% CI TBL	MeanTBL
Short latency	1	0.5	1	0.5	2	8 × 10–8	5	15	0–2	0.41
Mid latency	1	0.5	1	0.5	1	8 × 10–8	6	17	0–3	0.56
Long latency	1	0.5	1	0.5	0.5	8 × 10–8	10	20	0–5	0.87

Transmission, infectious period, and R0

I tested how different transmission and sampling rates impact the results of clustering and the TBL distributions. In these analyses, I set the death rate σ to zero, so that all individuals are sampled, and the length of the infectious period is therefore determined by the sampling rate ε. I found that clustering results and TBL depend on the transmission rate, with higher values leading to higher clustering rates and shorter TBL. Conversely, the sampling rate had no major effect (Appendix 5). I tested whether the threshold on the minimum tree size could bias these analyses (the difference from Appendix 4 is that R0 can be different from one). I found that different thresholds on the minimum number of tips sampled in the simulated tree can influence clustering and TBL, but for the settings used in this study this effect was negligible, and the results were robust to different thresholds (Appendix 5).

Sampling and transmission rates are particularly interesting parameters because together they determine R0, and whether the epidemic will grow or shrink. R0 is the ratio between the rate at which infectious individuals are created (numerator) and the rate at which infectious individuals are removed (denominator). If the numerator is larger than the denominator the epidemic will grow (R0 >1). Conversely, if the denominator is larger than the numerator the epidemic will shrink (R0 <1). Specifically, with the epidemiological model used in this study the numerator is the transmission rate (λ), while the denominator is determined by the sum of the sampling and death rates (σ+ε), so that R0 = λ / (σ +ε). Because only the numerator affects the results of clustering rates and TBL, changes in these two metrics will correlate with R0 only when the denominator (σ +ε) is fixed. To show this I tested three different sampling rates: ε=0.5, 1, and 2, corresponding to a median infectious period of ~16.6, 8.3, and 4.2 months, respectively. These values cover well the possible length of the symptomatic period estimated in different countries (Ku et al., 2021). For each value of ε, I considered three different transmission rates leading to three scenarios: a shrinking epidemic with R0=0.9, a stable epidemic with R0=1, and a growing epidemic with R0=1.1. All other parameters were constant in all simulations (Table 2). When considering scenarios with different sampling rates, and therefore infectious periods, I found no correlation between R0 and clustering rates, nor between R0 and TBL distributions. However, when the duration of the infection period, and all other parameters (except the transmission rate) were fixed, larger values of R0 correlated with higher clustering rates and shorter TBL (Table 2, Figure 3).

Table 2

Scenario(Median infectious period (months) - R0)	λ	ε	R0	σ	ψ	π	95% SNP threshold	100% SNP threshold	95% CI TBL	Mean TBL
Long infectious period, shrinking (17–0.9)	0.45	0.5	0.9	0	1	8 × 10–8	7	16	0–3	0.62
Long infectious period, stable (17 - 1)	0.5	0.5	1	0	1	8 × 10–8	6	16	0–3	0.59
Long infectious period, growing (17–1.1)	0.55	0.5	1.1	0	1	8 × 10–8	6	16	0–3	0.57
Medium infectious period, shrinking (8–0.9)	0.9	1	0.9	0	1	8 × 10–8	5	16	0–3	0.41
Medium infectious period, stable (8 - 1)	1	1	1	0	1	8 × 10–8	5	15	0–2	0.38
Medium infectious period, growing (8–1.1)	1.1	1	1.1	0	1	8 × 10–8	4	16	0–2	0.37
Short infectious period, shrinking (4–0.9)	1.8	2	0.9	0	1	8 × 10–8	5	14	0–2	0.31
Short infectious period, stable (4 - 1)	2	2	1	0	1	8 × 10–8	4	15	0–2	0.29
Short infectious period, growing (4–1.1)	2.2	2	1.1	0	1	8 × 10–8	3	17	0–2	0.27

Figure 3

Download asset Open asset

An example

The results reported above have important implications on the interpretation of molecular epidemiology studies of MTB. To illustrate this in practice let’s consider an example in which two different strains of MTB are causing an epidemic in the same area. The two strains have different epidemiological and biological characteristics. We can think about it like two lineages of MTB, but it could also be two different clones belonging to the same lineage. MTB type 1 is expanding, with a R0 of 1.1, and it is characterized by a slow disease progression, with a median latency period of about one year. Type 1 populations are expected to increase by ~150% every 10 years. Conversely type 2 is shrinking, with a R0 of 0.9, and it is characterized by a shorter median latency period, approximately 5 months. Type 2 populations are expected to shrink by ~50% every 10 years. In addition, type 1 and 2 have moderate differences in their rate of molecular evolution, with a clock rate of respectively 1 × 10^–7 and 7 × 10^–8 nucleotide changes per site per year. In all other aspects the two types are identical. Obviously, under this scenario, type 1 is much more concerning for public health compared to type 2, at least in the long term. I repeated the same analysis presented above for the two types (Table 3).

Table 3

Scenario	λ	ε	R0	σ	ψ	π	95% SNP threshold	100% SNP threshold	95% CI TBL	Mean TBL
Type 1	1.1	0.5	1.1	0.5	0.7	1 × 10–7	8	21	0–4	0.82
Type 2	0.9	0.5	0.9	0.5	1.7	7 × 10–8	5	13	0–2	0.40

Type 2 showed shorter terminal branches: the mean of the TBL distribution was 0.84 and 0.4 SNPs for type 1 and type 2, respectively. Moreover, type 2 had higher clustering rates (Figure 4), and the 95% sensitivity threshold was 8 and 5 SNPs for type 1 and 2, respectively. A typical interpretation of this data would be that type 2 is transmitting more than type 1, because of the shorter terminal branches, and the higher clustering rates. Alternatively, if we picked a specific threshold, we would find that type 2 strains are associated with clustering (for lower thresholds), or that there are no differences in clustering among the two types (for higher thresholds). In any case, a classic molecular epidemiology analysis would conclude that type 2 is transmitting more than type 1, or that there are no differences in transmission between the two types. However, the shorter TBL and larger clustering rates of type 2 are due to its shorter latency and lower clock rate, not to increased transmission. Type 2 has a lower transmission rate compared to type 1, and it is bound to extinction, while type 1 is growing exponentially. This example shows the pitfalls of TBL and clustering analyses to study transmission in MTB epidemics. Admittedly, the simulation parameters for this example were picked to highlight the potential problems. However, the epidemiological characteristics of circulating MTB strains are normally not known, and the differences in latency and clock rates used here are possible. The length of the latent period and the clock rates are well within the range of values estimated with different data sets (Menardo et al., 2019; Ku et al., 2021). Overall, these results show that the standard interpretation of clustering results and TBL distributions can potentially be misleading in some epidemiological settings.

Figure 4

Download asset Open asset

The current manuscript is a computational study, so no data have been generated for this manuscript. Modelling code and simulation results are available at https://github.com/fmenardo/sim_cluster_MTB, (copy archived at swh:1:rev:aa2e0bb7629c46e64a099247d225466615b55b07).

1. 1. Ali S
   2. Khan MT
   3. Anwar Sheed K
   4. Khan MM
   5. Hasan F

   (2019) Spoligotyping analysis of Mycobacterium tuberculosis in Khyber Pakhtunkhwa area, Pakistan

   Infection and Drug Resistance 12:1363–1369.

   https://doi.org/10.2147/IDR.S198314

   • PubMed
   • Google Scholar
2. 1. Asare P
   2. Otchere ID
   3. Bedeley E
   4. Brites D
   5. Loiseau C
   6. Baddoo NA
   7. Asante-Poku A
   8. Osei-Wusu S
   9. Prah DA
   10. Borrell S
   11. Reinhard M
   12. Forson A
   13. Koram KA
   14. Gagneux S
   15. Yeboah-Manu D

   (2020) Whole Genome Sequencing and Spatial Analysis Identifies Recent Tuberculosis Transmission Hotspots in Ghana

   Frontiers in Medicine 7:161.

   https://doi.org/10.3389/fmed.2020.00161

   • PubMed
   • Google Scholar
3. 1. Balaban M
   2. Moshiri N
   3. Mai U
   4. Jia X
   5. Mirarab S

   (2019) TreeCluster: Clustering biological sequences using phylogenetic trees

   PLOS ONE 14:e0221068.

   https://doi.org/10.1371/journal.pone.0221068

   • PubMed
   • Google Scholar
4. 1. Blackwell GA
   2. Hunt M
   3. Malone KM
   4. Lima L
   5. Horesh G
   6. Alako BTF
   7. Thomson NR
   8. Iqbal Z

   (2021) Exploring bacterial diversity via a curated and searchable snapshot of archived DNA sequences

   PLOS Biology 19:e3001421.

   https://doi.org/10.1371/journal.pbio.3001421

   • PubMed
   • Google Scholar
5. 1. Bottai D
   2. Frigui W
   3. Sayes F
   4. Di Luca M
   5. Spadoni D
   6. Pawlik A
   7. Zoppo M
   8. Orgeur M
   9. Khanna V
   10. Hardy D
   11. Mangenot S
   12. Barbe V
   13. Medigue C
   14. Ma L
   15. Bouchier C
   16. Tavanti A
   17. Larrouy-Maumus G
   18. Brosch R

   (2020) TbD1 deletion as a driver of the evolutionary success of modern epidemic Mycobacterium tuberculosis lineages

   Nature Communications 11:1–14.

   https://doi.org/10.1038/s41467-020-14508-5

   • PubMed
   • Google Scholar
6. 1. Bouckaert R
   2. Vaughan TG
   3. Barido-Sottani J
   4. Duchêne S
   5. Fourment M
   6. Gavryushkina A
   7. Heled J
   8. Jones G
   9. Kühnert D
   10. De Maio N
   11. Matschiner M
   12. Mendes FK
   13. Müller NF
   14. Ogilvie HA
   15. du Plessis L
   16. Popinga A
   17. Rambaut A
   18. Rasmussen D
   19. Siveroni I
   20. Suchard MA
   21. Wu C-H
   22. Xie D
   23. Zhang C
   24. Stadler T
   25. Drummond AJ

   (2019) BEAST 2.5: An advanced software platform for Bayesian evolutionary analysis

   PLOS Computational Biology 15:e1006650.

   https://doi.org/10.1371/journal.pcbi.1006650

   • PubMed
   • Google Scholar
7. 1. Brites D
   2. Loiseau C
   3. Menardo F
   4. Borrell S
   5. Boniotti MB
   6. Warren R
   7. Dippenaar A
   8. Parsons SDC
   9. Beisel C
   10. Behr MA
   11. Fyfe JA
   12. Coscolla M
   13. Gagneux S

   (2018) A New Phylogenetic Framework for the Animal-Adapted Mycobacterium tuberculosis Complex

   Frontiers in Microbiology 9:2820.

   https://doi.org/10.3389/fmicb.2018.02820

   • PubMed
   • Google Scholar
8. 1. Chen Y-Y
   2. Chang J-R
   3. Huang W-F
   4. Hsu C-H
   5. Cheng H-Y
   6. Sun J-R
   7. Kuo S-C
   8. Su I-J
   9. Lin M-S
   10. Chen W
   11. Dou H-Y

   (2017) Genetic diversity of the Mycobacterium tuberculosis East African-Indian family in three tropical Asian countries

   Journal of Microbiology, Immunology, and Infection = Wei Mian Yu Gan Ran Za Zhi 50:886–892.

   https://doi.org/10.1016/j.jmii.2015.10.012

   • PubMed
   • Google Scholar
9. 1. Chihota VN
   2. Niehaus A
   3. Streicher EM
   4. Wang X
   5. Sampson SL
   6. Mason P
   7. Källenius G
   8. Mfinanga SG
   9. Pillay M
   10. Klopper M
   11. Kasongo W
   12. Behr MA
   13. Gey van Pittius NC
   14. van Helden PD
   15. Couvin D
   16. Rastogi N
   17. Warren RM

   (2018) Geospatial distribution of Mycobacterium tuberculosis genotypes in Africa

   PLOS ONE 13:e0200632.

   https://doi.org/10.1371/journal.pone.0200632

   • PubMed
   • Google Scholar
10. 1. Cock PJA
    2. Antao T
    3. Chang JT
    4. Chapman BA
    5. Cox CJ
    6. Dalke A
    7. Friedberg I
    8. Hamelryck T
    9. Kauff F
    10. Wilczynski B
    11. de Hoon MJL

    (2009) Biopython: freely available Python tools for computational molecular biology and bioinformatics

    Bioinformatics (Oxford, England) 25:1422–1423.

    https://doi.org/10.1093/bioinformatics/btp163

    • PubMed
    • Google Scholar
11. 1. Cox H
    2. Salaam-Dreyer Z
    3. Goig GA
    4. Nicol MP
    5. Menardo F
    6. Dippenaar A
    7. Mohr-Holland E
    8. Daniels J
    9. Cudahy PGT
    10. Borrell S
    11. Reinhard M
    12. Doetsch A
    13. Beisel C
    14. Reuter A
    15. Furin J
    16. Gagneux S
    17. Warren RM

    (2021) Potential contribution of HIV during first-line tuberculosis treatment to subsequent rifampicin-monoresistant tuberculosis and acquired tuberculosis drug resistance in South Africa: a retrospective molecular epidemiology study

    The Lancet. Microbe 2:e584–e593.

    https://doi.org/10.1016/S2666-5247(21)00144-0

    • PubMed
    • Google Scholar
12. 1. de Jong BC
    2. Hill PC
    3. Aiken A
    4. Awine T
    5. Antonio M
    6. Adetifa IM
    7. Jackson-Sillah DJ
    8. Fox A
    9. Deriemer K
    10. Gagneux S
    11. Borgdorff MW
    12. McAdam KPWJ
    13. Corrah T
    14. Small PM
    15. Adegbola RA

    (2008) Progression to active tuberculosis, but not transmission, varies by Mycobacterium tuberculosis lineage in The Gambia

    The Journal of Infectious Diseases 198:1037–1043.

    https://doi.org/10.1086/591504

    • PubMed
    • Google Scholar
13. 1. Didelot X
    2. Fraser C
    3. Gardy J
    4. Colijn C

    (2017) Genomic Infectious Disease Epidemiology in Partially Sampled and Ongoing Outbreaks

    Molecular Biology and Evolution 34:997–1007.

    https://doi.org/10.1093/molbev/msw275

    • PubMed
    • Google Scholar
14. 1. Eldholm V
    2. Pettersson JHO
    3. Brynildsrud OB
    4. Kitchen A
    5. Rasmussen EM
    6. Lillebaek T
    7. Rønning JO
    8. Crudu V
    9. Mengshoel AT
    10. Debech N
    11. Alfsnes K
    12. Bohlin J
    13. Pepperell CS
    14. Balloux F

    (2016) Armed conflict and population displacement as drivers of the evolution and dispersal of Mycobacterium tuberculosis

    PNAS 113:13881–13886.

    https://doi.org/10.1073/pnas.1611283113

    • PubMed
    • Google Scholar
15. 1. Firdessa R
    2. Berg S
    3. Hailu E
    4. Schelling E
    5. Gumi B
    6. Erenso G
    7. Gadisa E
    8. Kiros T
    9. Habtamu M
    10. Hussein J
    11. Zinsstag J
    12. Robertson BD
    13. Ameni G
    14. Lohan AJ
    15. Loftus B
    16. Comas I
    17. Gagneux S
    18. Tschopp R
    19. Yamuah L
    20. Hewinson G
    21. Gordon SV
    22. Young DB
    23. Aseffa A

    (2013) Mycobacterial lineages causing pulmonary and extrapulmonary tuberculosis, Ethiopia

    Emerging Infectious Diseases 19:460–463.

    https://doi.org/10.3201/eid1903.120256

    • PubMed
    • Google Scholar
16. 1. Freschi L
    2. Vargas R Jr
    3. Husain A
    4. Kamal SMM
    5. Skrahina A
    6. Tahseen S
    7. Ismail N
    8. Barbova A
    9. Niemann S
    10. Cirillo DM
    11. Dean AS
    12. Zignol M
    13. Farhat MR

    (2021) Population structure, biogeography and transmissibility of Mycobacterium tuberculosis

    Nature Communications 12:1–11.

    https://doi.org/10.1038/s41467-021-26248-1

    • PubMed
    • Google Scholar
17. 1. Guerra-Assunção JA
    2. Crampin AC
    3. Houben RMGJ
    4. Mzembe T
    5. Mallard K
    6. Coll F
    7. Khan P
    8. Banda L
    9. Chiwaya A
    10. Pereira RPA
    11. McNerney R
    12. Fine PEM
    13. Parkhill J
    14. Clark TG
    15. Glynn JR

    (2015) Large-scale whole genome sequencing of M. tuberculosis provides insights into transmission in a high prevalence area

    eLife 4:e05166.

    https://doi.org/10.7554/eLife.05166

    • PubMed
    • Google Scholar
18. 1. Guthrie JL
    2. Gardy JL

    (2017) A brief primer on genomic epidemiology: lessons learned from Mycobacterium tuberculosis

    Annals of the New York Academy of Sciences 1388:59–77.

    https://doi.org/10.1111/nyas.13273

    • PubMed
    • Google Scholar
19. 1. Gygli SM
    2. Loiseau C
    3. Jugheli L
    4. Adamia N
    5. Trauner A
    6. Reinhard M
    7. Ross A
    8. Borrell S
    9. Aspindzelashvili R
    10. Maghradze N
    11. Reither K
    12. Beisel C
    13. Tukvadze N
    14. Avaliani Z
    15. Gagneux S

    (2021) Prisons as ecological drivers of fitness-compensated multidrug-resistant Mycobacterium tuberculosis

    Nature Medicine 27:1171–1177.

    https://doi.org/10.1038/s41591-021-01358-x

    • PubMed
    • Google Scholar
20. 1. Hanekom M
    2. Gey van Pittius NC
    3. McEvoy C
    4. Victor TC
    5. Van Helden PD
    6. Warren RM

    (2011) Mycobacterium tuberculosis Beijing genotype: a template for success

    Tuberculosis (Edinburgh, Scotland) 91:510–523.

    https://doi.org/10.1016/j.tube.2011.07.005

    • PubMed
    • Google Scholar
21. 1. Hang NTL
    2. Hijikata M
    3. Maeda S
    4. Thuong PH
    5. Ohashi J
    6. Van Huan H
    7. Hoang NP
    8. Miyabayashi A
    9. Cuong VC
    10. Seto S
    11. Van Hung N
    12. Keicho N

    (2019) Whole genome sequencing, analyses of drug resistance-conferring mutations, and correlation with transmission of Mycobacterium tuberculosis carrying katG-S315T in Hanoi, Vietnam

    Scientific Reports 9:1–14.

    https://doi.org/10.1038/s41598-019-51812-7

    • PubMed
    • Google Scholar
22. 1. Hatherell HA
    2. Colijn C
    3. Stagg HR
    4. Jackson C
    5. Winter JR
    6. Abubakar I

    (2016) Interpreting whole genome sequencing for investigating tuberculosis transmission: a systematic review

    BMC Medicine 14:21.

    https://doi.org/10.1186/s12916-016-0566-x

    • PubMed
    • Google Scholar
23. 1. Holt KE
    2. McAdam P
    3. Thai PVK
    4. Thuong NTT
    5. Ha DTM
    6. Lan NN
    7. Lan NH
    8. Nhu NTQ
    9. Hai HT
    10. Ha VTN
    11. Thwaites G
    12. Edwards DJ
    13. Nath AP
    14. Pham K
    15. Ascher DB
    16. Farrar J
    17. Khor CC
    18. Teo YY
    19. Inouye M
    20. Caws M
    21. Dunstan SJ

    (2018) Frequent transmission of the Mycobacterium tuberculosis Beijing lineage and positive selection for the EsxW Beijing variant in Vietnam

    Nature Genetics 50:849–856.

    https://doi.org/10.1038/s41588-018-0117-9

    • PubMed
    • Google Scholar
24. 1. Huerta-Cepas J
    2. Serra F
    3. Bork P

    (2016) ETE 3: Reconstruction, Analysis, and Visualization of Phylogenomic Data

    Molecular Biology and Evolution 33:1635–1638.

    https://doi.org/10.1093/molbev/msw046

    • PubMed
    • Google Scholar
25. 1. Kendall EA
    2. Shrestha S
    3. Dowdy DW

    (2021) The Epidemiological Importance of Subclinical Tuberculosis: A Critical Reappraisal

    American Journal of Respiratory and Critical Care Medicine 203:168–174.

    https://doi.org/10.1164/rccm.202006-2394PP

    • PubMed
    • Google Scholar
26. 1. Kozlov AM
    2. Darriba D
    3. Flouri T
    4. Morel B
    5. Stamatakis A

    (2019) RAxML-NG: a fast, scalable and user-friendly tool for maximum likelihood phylogenetic inference

    Bioinformatics (Oxford, England) 35:4453–4455.

    https://doi.org/10.1093/bioinformatics/btz305

    • PubMed
    • Google Scholar
27. 1. Ku C-C
    2. MacPherson P
    3. Khundi M
    4. Nzawa Soko RH
    5. Feasey HRA
    6. Nliwasa M
    7. Horton KC
    8. Corbett EL
    9. Dodd PJ

    (2021) Durations of asymptomatic, symptomatic, and care-seeking phases of tuberculosis disease with a Bayesian analysis of prevalence survey and notification data

    BMC Medicine 19:1–13.

    https://doi.org/10.1186/s12916-021-02128-9

    • PubMed
    • Google Scholar
28. 1. Kühnert D.
    2. Stadler T
    3. Vaughan TG
    4. Drummond AJ

    (2016) Phylodynamics with Migration: A Computational Framework to Quantify Population Structure from Genomic Data

    Molecular Biology and Evolution 33:2102–2116.

    https://doi.org/10.1093/molbev/msw064

    • PubMed
    • Google Scholar
29. 1. Kühnert D
    2. Coscolla M
    3. Brites D
    4. Stucki D
    5. Metcalfe J
    6. Fenner L
    7. Gagneux S
    8. Stadler T

    (2018) Tuberculosis outbreak investigation using phylodynamic analysis

    Epidemics 25:47–53.

    https://doi.org/10.1016/j.epidem.2018.05.004

    • PubMed
    • Google Scholar
30. 1. Liu Q
    2. Liu H
    3. Shi L
    4. Gan M
    5. Zhao X
    6. Lyu LD
    7. Takiff HE
    8. Wan K
    9. Gao Q

    (2021) Local adaptation of Mycobacterium tuberculosis on the Tibetan Plateau

    PNAS 118:e2017831118.

    https://doi.org/10.1073/pnas.2017831118

    • PubMed
    • Google Scholar
31. 1. López MG
    2. Dogba JB
    3. Torres-Puente M
    4. Goig GA
    5. Moreno-Molina M
    6. Villamayor LM
    7. Cadmus S
    8. Comas I

    (2020) Tuberculosis in Liberia: high multidrug-resistance burden, transmission and diversity modelled by multiple importation events

    Microbial Genomics 6:325.

    https://doi.org/10.1099/mgen.0.000325

    • PubMed
    • Google Scholar
32. 1. Louca S
    2. McLaughlin A
    3. MacPherson A
    4. Joy JB
    5. Pennell MW

    (2021) Fundamental Identifiability Limits in Molecular Epidemiology

    Molecular Biology and Evolution 38:4010–4024.

    https://doi.org/10.1093/molbev/msab149

    • PubMed
    • Google Scholar
33. 1. Meehan CJ
    2. Moris P
    3. Kohl TA
    4. Pečerska J
    5. Akter S
    6. Merker M
    7. Utpatel C
    8. Beckert P
    9. Gehre F
    10. Lempens P
    11. Stadler T
    12. Kaswa MK
    13. Kühnert D
    14. Niemann S
    15. de Jong BC

    (2018) The relationship between transmission time and clustering methods in Mycobacterium tuberculosis epidemiology

    EBioMedicine 37:410–416.

    https://doi.org/10.1016/j.ebiom.2018.10.013

    • PubMed
    • Google Scholar
34. 1. Menardo F.
    2. Duchêne S
    3. Brites D
    4. Gagneux S

    (2019) The molecular clock of Mycobacterium tuberculosis

    PLOS Pathogens 15:e1008067.

    https://doi.org/10.1371/journal.ppat.1008067

    • PubMed
    • Google Scholar
35. 1. Menardo F
    2. Rutaihwa LK
    3. Zwyer M
    4. Borrell S
    5. Comas I
    6. Conceição EC
    7. Coscolla M
    8. Cox H
    9. Joloba M
    10. Dou HY
    11. Feldmann J
    12. Fenner L
    13. Fyfe J
    14. Gao Q
    15. García de Viedma D
    16. Garcia-Basteiro AL
    17. Gygli SM
    18. Hella J
    19. Hiza H
    20. Jugheli L
    21. Kamwela L
    22. Kato-Maeda M
    23. Liu Q
    24. Ley SD
    25. Loiseau C
    26. Mahasirimongkol S
    27. Malla B
    28. Palittapongarnpim P
    29. Rakotosamimanana N
    30. Rasolofo V
    31. Reinhard M
    32. Reither K
    33. Sasamalo M
    34. Silva Duarte R
    35. Sola C
    36. Suffys P
    37. Batista Lima KV
    38. Yeboah-Manu D
    39. Beisel C
    40. Brites D
    41. Gagneux S

    (2021) Local adaptation in populations of Mycobacterium tuberculosis endemic to the Indian Ocean Rim

    F1000Research 10:60.

    https://doi.org/10.12688/f1000research.28318.2

    • PubMed
    • Google Scholar
36. Software

    1. Menardo F

    (2022) Pipeline to simulate the evolution of genome sequences of MTB under different conditions, version swh:1:rev:aa2e0bb7629c46e64a099247d225466615b55b07

    Software Heritage.

    https://archive.softwareheritage.org/swh:1:dir:3d3b0f76e5e4baa8e2500cbcba392dc0c33311ba;origin=https://github.com/fmenardo/sim_cluster_MTB;visit=swh:1:snp:bb5bad061d6cb101966bc93d3b3b0c13ee82aa47;anchor=swh:1:rev:aa2e0bb7629c46e64a099247d225466615b55b07
37. 1. Merker M
    2. Barbier M
    3. Cox H
    4. Rasigade JP
    5. Feuerriegel S
    6. Kohl TA
    7. Diel R
    8. Borrell S
    9. Gagneux S
    10. Nikolayevskyy V
    11. Andres S
    12. Nübel U
    13. Supply P
    14. Wirth T
    15. Niemann S

    (2018) Compensatory evolution drives multidrug-resistant tuberculosis in Central Asia

    eLife 7:e38200.

    https://doi.org/10.7554/eLife.38200

    • PubMed
    • Google Scholar
38. 1. Merker M
    2. Egbe NF
    3. Ngangue YR
    4. Vuchas C
    5. Kohl TA
    6. Dreyer V
    7. Kuaban C
    8. Noeske J
    9. Niemann S
    10. Sander MS

    (2021) Transmission patterns of rifampicin resistant Mycobacterium tuberculosis complex strains in Cameroon: a genomic epidemiological study

    BMC Infectious Diseases 21:1–10.

    https://doi.org/10.1186/s12879-021-06593-8

    • PubMed
    • Google Scholar
39. 1. Mulenga C
    2. Shamputa IC
    3. Mwakazanga D
    4. Kapata N
    5. Portaels F
    6. Rigouts L

    (2010) Diversity of Mycobacterium tuberculosis genotypes circulating in Ndola, Zambia

    BMC Infectious Diseases 10:1–9.

    https://doi.org/10.1186/1471-2334-10-177

    • PubMed
    • Google Scholar
40. 1. Netikul T
    2. Palittapongarnpim P
    3. Thawornwattana Y
    4. Plitphonganphim S

    (2021) Estimation of the global burden of Mycobacterium tuberculosis lineage 1

    Infection, Genetics and Evolution 91:104802.

    https://doi.org/10.1016/j.meegid.2021.104802

    • PubMed
    • Google Scholar
41. 1. Nikolayevskyy V
    2. Niemann S
    3. Anthony R
    4. van Soolingen D
    5. Tagliani E
    6. Ködmön C
    7. van der Werf MJ
    8. Cirillo DM

    (2019) Role and value of whole genome sequencing in studying tuberculosis transmission

    Clinical Microbiology and Infection 25:1377–1382.

    https://doi.org/10.1016/j.cmi.2019.03.022

    • PubMed
    • Google Scholar
42. 1. Page AJ
    2. Taylor B
    3. Delaney AJ
    4. Soares J
    5. Seemann T
    6. Keane JA
    7. Harris SR

    (2016) SNP-sites: rapid efficient extraction of SNPs from multi-FASTA alignments

    Microbial Genomics 2:e000056.

    https://doi.org/10.1099/mgen.0.000056

    • PubMed
    • Google Scholar
43. 1. Peters JS
    2. Ismail N
    3. Dippenaar A
    4. Ma S
    5. Sherman DR
    6. Warren RM
    7. Kana BD

    (2020) Genetic Diversity in Mycobacterium tuberculosis Clinical Isolates and Resulting Outcomes of Tuberculosis Infection and Disease

    Annual Review of Genetics 54:511–537.

    https://doi.org/10.1146/annurev-genet-022820-085940

    • PubMed
    • Google Scholar
44. 1. Poon AFY

    (2016) Impacts and shortcomings of genetic clustering methods for infectious disease outbreaks

    Virus Evolution 2:vew031.

    https://doi.org/10.1093/ve/vew031

    • PubMed
    • Google Scholar
45. 1. Rambaut A
    2. Grass NC

    (1997) Seq-Gen: an application for the Monte Carlo simulation of DNA sequence evolution along phylogenetic trees

    Bioinformatics 13:235–238.

    https://doi.org/10.1093/bioinformatics/13.3.235

    • Google Scholar
46. 1. Rutaihwa LK
    2. Menardo F
    3. Stucki D
    4. Gygli SM
    5. Ley SD
    6. Malla B
    7. Feldmann J
    8. Borrell S
    9. Beisel C
    10. Middelkoop K
    11. Carter EJ
    12. Diero L
    13. Ballif M
    14. Jugheli L
    15. Reither K
    16. Fenner L
    17. Brites D
    18. Gagneux S

    (2019) Multiple Introductions of Mycobacterium tuberculosis Lineage 2–Beijing Into Africa Over Centuries

    Frontiers in Ecology and Evolution 7:112.

    https://doi.org/10.3389/fevo.2019.00112

    • Google Scholar
47. 1. Schopfer K
    2. Rieder HL
    3. Steinlin-Schopfer JF
    4. van Soolingen D
    5. Bodmer T
    6. Chantana Y
    7. Studer P
    8. Laurent D
    9. Zwahlen M
    10. Richner B

    (2015) Molecular epidemiology of tuberculosis in Cambodian children

    Epidemiology and Infection 143:910–921.

    https://doi.org/10.1017/S0950268814001769

    • PubMed
    • Google Scholar
48. 1. Shuaib YA
    2. Khalil EAG
    3. Wieler LH
    4. Schaible UE
    5. Bakheit MA
    6. Mohamed-Noor SE
    7. Abdalla MA
    8. Kerubo G
    9. Andres S
    10. Hillemann D
    11. Richter E
    12. Kranzer K
    13. Niemann S
    14. Merker M

    (2020) Mycobacterium tuberculosis Complex Lineage 3 as Causative Agent of Pulmonary Tuberculosis, Eastern Sudan¹

    Emerging Infectious Diseases 26:427–436.

    https://doi.org/10.3201/eid2603.191145

    • PubMed
    • Google Scholar
49. 1. Sobkowiak B
    2. Banda L
    3. Mzembe T
    4. Crampin AC
    5. Glynn JR
    6. Clark TG

    (2020) Bayesian reconstruction of Mycobacterium tuberculosis transmission networks in a high incidence area over two decades in Malawi reveals associated risk factors and genomic variants

    Microbial Genomics 6:361.

    https://doi.org/10.1099/mgen.0.000361

    • PubMed
    • Google Scholar
50. Thesis

    1. Somphavong S

    (2018)

    Molecular epidemiology of Mycobacterium tuberculosis and antibiotic resistance in Lao PDR

    Université Montpellier).

    • Google Scholar
51. 1. Stimson J
    2. Gardy J
    3. Mathema B
    4. Crudu V
    5. Cohen T
    6. Colijn C

    (2019) Beyond the SNP Threshold: Identifying Outbreak Clusters Using Inferred Transmissions

    Molecular Biology and Evolution 36:587–603.

    https://doi.org/10.1093/molbev/msy242

    • PubMed
    • Google Scholar
52. 1. Torres Ortiz A
    2. Coronel J
    3. Vidal JR
    4. Bonilla C
    5. Moore DAJ
    6. Gilman RH
    7. Balloux F
    8. Kon OM
    9. Didelot X
    10. Grandjean L

    (2021) Genomic signatures of pre-resistance in Mycobacterium tuberculosis

    Nature Communications 12:7312.

    https://doi.org/10.1038/s41467-021-27616-7

    • PubMed
    • Google Scholar
53. 1. Tulu B
    2. Ameni G

    (2018) Spoligotyping based genetic diversity of Mycobacterium tuberculosis in Ethiopia: a systematic review

    BMC Infectious Diseases 18:1–10.

    https://doi.org/10.1186/s12879-018-3046-4

    • PubMed
    • Google Scholar
54. 1. Vaughan TG
    2. Drummond AJ

    (2013) A stochastic simulator of birth-death master equations with application to phylodynamics

    Molecular Biology and Evolution 30:1480–1493.

    https://doi.org/10.1093/molbev/mst057

    • PubMed
    • Google Scholar
55. 1. Vaziri F
    2. Kohl TA
    3. Ghajavand H
    4. Kargarpour Kamakoli M
    5. Merker M
    6. Hadifar S
    7. Khanipour S
    8. Fateh A
    9. Masoumi M
    10. Siadat SD
    11. Niemann S

    (2019) Genetic Diversity of Multi- and Extensively Drug-Resistant Mycobacterium tuberculosis Isolates in the Capital of Iran, Revealed by Whole-Genome Sequencing

    Journal of Clinical Microbiology 57:e01477-18.

    https://doi.org/10.1128/JCM.01477-18

    • PubMed
    • Google Scholar
56. 1. Volz EM
    2. Siveroni I

    (2018) Bayesian phylodynamic inference with complex models

    PLOS Computational Biology 14:e1006546.

    https://doi.org/10.1371/journal.pcbi.1006546

    • PubMed
    • Google Scholar
57. 1. Walker TM
    2. Merker M
    3. Knoblauch AM
    4. Helbling P
    5. Schoch OD
    6. van der Werf MJ
    7. Kranzer K
    8. Fiebig L
    9. Kröger S
    10. Haas W
    11. Hoffmann H
    12. Indra A
    13. Egli A
    14. Cirillo DM
    15. Robert J
    16. Rogers TR
    17. Groenheit R
    18. Mengshoel AT
    19. Mathys V
    20. Haanperä M
    21. Soolingen DV
    22. Niemann S
    23. Böttger EC
    24. Keller PM
    25. MDR-TB Cluster Consortium

    (2018) A cluster of multidrug-resistant Mycobacterium tuberculosis among patients arriving in Europe from the Horn of Africa: A molecular epidemiological study

    The Lancet. Infectious Diseases 18:431–440.

    https://doi.org/10.1016/S1473-3099(18)30004-5

    • PubMed
    • Google Scholar
58. 1. Walter KS
    2. Dos Santos PCP
    3. Gonçalves TO
    4. da Silva BO
    5. da Silva Santos A
    6. de Cássia Leite A
    7. da Silva AM
    8. Figueira Moreira FM
    9. de Oliveira RD
    10. Lemos EF
    11. Cunha E
    12. Liu YE
    13. Ko AI
    14. Colijn C
    15. Cohen T
    16. Mathema B
    17. Croda J
    18. Andrews JR

    (2022) The role of prisons in disseminating tuberculosis in Brazil: A genomic epidemiology study

    Lancet Regional Health. Americas 9:100186.

    https://doi.org/10.1016/j.lana.2022.100186

    • PubMed
    • Google Scholar
59. 1. Wickham H

    (2011) ggplot2

    Wiley Interdisciplinary Reviews 3:180–185.

    https://doi.org/10.1002/wics.147

    • Google Scholar
60. 1. Wiens KE
    2. Woyczynski LP
    3. Ledesma JR
    4. Ross JM
    5. Zenteno-Cuevas R
    6. Goodridge A
    7. Ullah I
    8. Mathema B
    9. Djoba Siawaya JF
    10. Biehl MH
    11. Ray SE
    12. Bhattacharjee NV
    13. Henry NJ
    14. Reiner RC Jr
    15. Kyu HH
    16. Murray CJL
    17. Hay SI

    (2018) Global variation in bacterial strains that cause tuberculosis disease: a systematic review and meta-analysis

    BMC Medicine 16:1–13.

    https://doi.org/10.1186/s12916-018-1180-x

    • PubMed
    • Google Scholar
61. 1. Yang C
    2. Lu L
    3. Warren JL
    4. Wu J
    5. Jiang Q
    6. Zuo T
    7. Gan M
    8. Liu M
    9. Liu Q
    10. DeRiemer K
    11. Hong J
    12. Shen X
    13. Colijn C
    14. Guo X
    15. Gao Q
    16. Cohen T

    (2018) Internal migration and transmission dynamics of tuberculosis in Shanghai, China: an epidemiological, spatial, genomic analysis

    The Lancet. Infectious Diseases 18:788–795.

    https://doi.org/10.1016/S1473-3099(18)30218-4

    • PubMed
    • Google Scholar
62. 1. Yang C
    2. Sobkowiak B
    3. Naidu V
    4. Codreanu A
    5. Ciobanu N
    6. Gunasekera KS
    7. Chitwood MH
    8. Alexandru S
    9. Bivol S
    10. Russi M
    11. Havumaki J
    12. Cudahy P
    13. Fosburgh H
    14. Allender CJ
    15. Centner H
    16. Engelthaler DM
    17. Menzies NA
    18. Warren JL
    19. Crudu V
    20. Colijn C
    21. Cohen T

    (2022) Phylogeography and transmission of M. tuberculosis in Moldova: A prospective genomic analysis

    PLOS Medicine 19:e1003933.

    https://doi.org/10.1371/journal.pmed.1003933

    • PubMed
    • Google Scholar

Decision letter after peer review:

Thank you for submitting your article "Understanding drivers of phylogenetic clustering and terminal branch lengths distribution in epidemics of Mycobacterium tuberculosis" for consideration by eLife. Your article has been reviewed by 2 peer reviewers, and the evaluation has been overseen by a Reviewing Editor and a Senior Editor. The following individuals involved in the review of your submission have agreed to reveal their identity: Jason Andrews (Reviewer #1); Vegard Eldholm (Reviewer #2).

We agree that this paper is interesting, timely and raises important points. As is customary in eLife, the reviewers have discussed their critiques with one another. What follows below is the Reviewing Editor's edited compilation of the essential and ancillary points provided by reviewers in their critiques and in their interaction post-review. Please submit a revised version that addresses these concerns directly. Although we expect that you will address these comments in your response letter, we also need to see the corresponding revision clearly marked in the text of the manuscript. Some of the reviewers' comments may seem to be simple queries or challenges that do not prompt revisions to the text. Please keep in mind, however, that readers may have the same perspective as the reviewers. Therefore, it is essential that you attempt to amend or expand the text to clarify the narrative accordingly.

Essential revisions:

1) Both reviewers had some comments about the simulations for R0 or 1.1 and 0.9 – additional clarity on what impacts R0, and more substantively, can TBL and clustering rate help in estimating R0, if the infectious duration is held fixed? R1 notes "These issues should be addressed and/or the conclusions about lack of ability to infer R0 differences should be tempered a bit".

2) While both reviewers (and I) really like the idea of having an illustrative example showing how this could matter in a scenario with two co-circulating lineages in the same community, the differences seem too pronounced.

3) I have an additional comment (as reviewing editor):

You note that you select simulations with 100-2500 tips – how do you initialize the R0 < 1 simulations to end up with at least 100 tips? The tree size will presumably differ with different R0, as will the number of simulations you run that don't meet the size condition.

Given that you also sample in the last 10 years of the simulation. A very different fraction of the process will potentially end up in the last 10 years under the different parameters. Accordingly, the relationship between R0, TBL and clustering will be complex, and the fact that the results don't seem to be informative about R0 may be partly a result of this "simulation bias" and sampling period.

4) R1 notes that more comments on the sources of data that could help to disentangle the impact of differences in latency vs transmission would be helpful. Presumably that's a population-representative sampling that would enable knowing the fraction of cases that are of the two co-circulating types, over time? Or good estimates of the infectious duration? Or either, or both.

5) Please address minor points on clarity and communication raised by the reviewers. In addition, could you comment on the mutation rate during latency, and what might the result be if mutations were not occurring during latency at the rate that they do in active disease?

Reviewer #1 (Recommendations for the authors):

Overall, the manuscript was fairly clear and straightforward to interpret, with a few comments below. The main concern I had was with the overarching conclusion that TBLs are non-informative about R0. The results show that transmission rates are associated with TBL/clustering, but the argument is that R0 is not because it depends on the ratio between the transmission rate and removal rate, each of which can alter those metrics. But it would be reasonable to make assumptions about the infectious duration, for example within a lineage and/or location, enabling comparison between variants within a lineage, within lineage over time or space, etc. Holding infectious duration constant, comparisons of R0 would then be possible through TBL or clustering. Moreover, in some cases, the infectious duration can be directly estimated (i.e. by performing a prevalence study and comparing with notifications), enabling the relation of R0 to TBL/clustering directly through transmission rates. These issues should be addressed and/or the conclusions about lack of ability to infer R0 differences should be tempered a bit. This wouldn't diminish the importance of the study, as explaining the relationship between all of these parameters and these metrics is the valuable task that this article undertakes.

Reviewer #2 (Recommendations for the authors):

First, this is a cool and timely paper, and I really like the approach.

There are however some things I believe could be clarified and perhaps explored a little bit further.

Page 7 "An example": Drawing up an example of how different characteristics of co-circulating TB-types can result in contra-intuitive findings is great. In the example, a less transmissible type nevertheless has an R0 > 1 whereas a more transmissible variant has R0 < 1. Looking at table 2, it seems evident that the more transmissible Type 2 is contracting (R0<1) due to higher cure/death rates and higher sampling compared to the less transmissible Type 1, whereas it exhibits shorter TBLs and higher clustering due to shorter latency and a slower molecular clock compared to Type 1.

I have a few questions/comments which I hope will clear this up a little bit on my part:

- I think the example would be easier to follow if you explained which of the parameters influenced the R0, TBL and clustering rate, as I have tried to do above. If my summary is not correct, I guess that illustrates that this is a bit complex for many readers.

- This is a bit embarrassing, but what specifically should sampling rate be interpreted as here, and throughout the paper? If I get this right, sampling represents both observation and removal (the patient is sampled, and hence also cured and does not transmit further)?

- I think it would also be cool with a little description of the two types in more biological/medical lingo. If I understand correctly, Type 2 is more transmissible, and also characterized by rapid onset of disease (short latency), followed by rapid death, self-cure or health-seeking. This actually makes sense, and the only reason R0 < 1 is that health-seeking, cure and death rates are even more elevated than transmissibility?

Page 3: If I'm not wrong, the model used is a form of birth-death model? Perhaps this could be spelled out, and if not, explain how it differs from a B-D model.

Suppl figures: I believe the numbering of the suppl figures in relation to the text is wrong in quite a few instances

[Editors' note: further revisions were suggested prior to acceptance, as described below.]

Thank you for resubmitting your work entitled "Understanding drivers of phylogenetic clustering and terminal branch lengths distribution in epidemics of Mycobacterium tuberculosis" for further consideration by eLife. Your revised article has been evaluated by a Senior Editor and a Reviewing Editor.

The manuscript has been improved but there are some remaining issues that need to be addressed, as outlined below:

In response to the question about bias induced by having to choose among your simulations, I fear that you have misunderstood the question and therefore haven't addressed it. (The new exploration of the impact of sampling period is nice to have, though). You write: "Regarding the settings: all simulations are initialized as a single individual, simulations that result in less than 100 tips are discarded. For R0 < 1, a larger number of simulations needs to be discarded, as trees are tendentially smaller. "

This is a source of potentially large bias because the tree you obtain is likely not to meaningfully *have* an R0 of 0.9. To grow from 1 to 100, the mean number of offspring that is realised is definitely above 1 for an extended period, even if the parameter you set in the process was such that the *expected* number was less than 1. (I would be curious to see a birth-death simulation with an R0 of 0.9 that grows from 1 to 100 tips and last over 10 years). The probability of this growth is not zero, but it's small, so the trees you obtain are not at all representative of R0 = 0.9 trees. This means that when the clustering differs or fails to differ, you can't interpret that as being able to identify (or not) trees from R=0.9 vs R=1.1.

The tree sizes will be different too, and from coalescent theory, we know that this alone will change the patterns of genetic diversity (and hence the clustering). I would hazard that if you took simulations from identical parameters (birth, death, sampling, etc all the same) but in one set you chose the rejection criterion to have the trees need to grow bigger or last longer than in the other set, you would also find differences in clustering as a sole function of the cutoff criterion. Some portion of what you find is due probably due primarily to size and you could explore this, as well as not so dramatically biasing the trees you analyse by requiring extremely unlikely events before a tree gets into your sample.

Think of it the other way. If you took trees from a simulation with R0> 1 but you *only* looked at trees that went extinct before reaching 10 taxa, and then you did estimation or some other analysis on those trees, they would look very similar to trees with an R0 < 1 because those also die out. By removing the ones that grew you removed the information that R0 was set to be > 1. Conversely, here, by rejecting the vast majority of trees simulated under R0 < 1 that died out, you are removing the information that R0 was ever < 1.

I have asked for a revision because I think this point is not minor, as quite a bit of the paper is focused on clustering and R0 under simulations, and this same issue could impact many of them. It would be interesting to explore a genuinely declining population (rather than one that is declining in expectation because R0 < 1 but growing in its actual realisation because you start with 1 taxon and branch, rejecting those simulations that do not grow). However, initialising that population is challenging because the results would depend on the initial genetic diversity.

As a side note – under point 4 in the response, you note that "At least, in theory, all the parameters of the model can be estimated with phylodynamic analyses." But this paper finds a fundamental unidentifiability that suggests that these parameters are not identifiable: https://academic.oup.com/mbe/article/38/9/4010/6278301. They use likelihoods, not clustering, but since the likelihoods and portion clustered are fundamentally based on the branching times in the phylogenies it seems that their results would probably carry over.

Essential revisions:

1) Both reviewers had some comments about the simulations for R0 or 1.1 and 0.9 – additional clarity on what impacts R0, and more substantively, can TBL and clustering rate help in estimating R0, if the infectious duration is held fixed? R1 notes "These issues should be addressed and/or the conclusions about lack of ability to infer R0 differences should be tempered a bit".

Regarding the main point (“can TBL and clustering rate help in estimating R0, if the infectious duration is held fixed?”):

In empirical studies, TBL and clustering rates can be used to estimate R0 only if the clock rate, the sampling proportion, the sampling period, the length of the infectious period, and the length of the latency period are ALL known or held fixed. If the infectious period is fixed, but the other factors are not, R0 would not necessarily correlate with clustering rates and TBL. To show this more clearly, I modified the example (Page 6), which now consists of two MTB “types” with identical sampling and death rates (and therefore identical infectious period), but different latency, transmission and clock rates. One type has R0 = 0.9, and the other has R0 = 1.1 (because of different transmission rates). Again, the shrinking type resulted in larger clustering rates and shorter TBL (because of differences in clock rates and latency). This is described in the response to point (2)

I also expanded the description of what influences R0 (Page 4, Ln 146-153):

R0 is the ratio between the rate at which infectious individuals are created (numerator) and the rate at which infectious individuals are removed (denominator). If the numerator is larger than the denominator the epidemic will grow (R0 > 1). Conversely, if the denominator is larger than the numerator the epidemic will shrink (R0 < 1). Specifically, with the epidemiological model used in this study the numerator is the transmission rate (λ), while the denominator is determined by the sum of the sampling rate and the death rate (σ+ε). R0 = λ / (σ +ε).

In non-mathematical terms R0 is determined by the relation between the transmission rate per unit of time and the average length of the infectious period (which is inversely proportional to the removal rate σ +ε). So that a large transmission rate (per unit of time) can result in a shrinking epidemic if the infectious period is very short. Conversely, a low transmission rate per unit of time can result in a growing epidemic if the infectious period is sufficiently long.

I think some confusion on these issues was generated by the order in which I presented the results. To make these points clearer I modified the structure of the Results section (Pages. 3-5). I also amended the abstract (Ln 21) and expanded the discussion (Page 7, Ln 234-242).

2) While both reviewers (and I) really like the idea of having an illustrative example showing how this could matter in a scenario with two co-circulating lineages in the same community, the differences seem too pronounced.

The main criticism of R1 was that in the example the length of the infectious period was drastically different between the two types, and that there is no direct evidence for differences of this magnitude between strains co-circulating in the same setting. I agree with R1 that the evidence for this is indirect. Therefore, I repeated the analysis assuming identical sampling and death rates for the two types (and therefore the same infectious period), also in this case the shrinking population resulted in higher clustering rates and shorter TBL, which were caused by shorter latency and lower clock rate (Page 6).

3) I have an additional comment (as reviewing editor):

You note that you select simulations with 100-2500 tips – how do you initialize the R0 < 1 simulations to end up with at least 100 tips? The tree size will presumably differ with different R0, as will the number of simulations you run that don't meet the size condition.

Given that you also sample in the last 10 years of the simulation. A very different fraction of the process will potentially end up in the last 10 years under the different parameters. Accordingly, the relationship between R0, TBL and clustering will be complex, and the fact that the results don't seem to be informative about R0 may be partly a result of this "simulation bias" and sampling period.

Thank you for this comment. In the first version of the manuscript, I did not consider a factor that could influence the results of the analysis: the sampling period.

To answer this comment, I performed an additional analysis. First, I tested whether different sampling periods can lead to differences in TBL and clustering rates (all else being equal). I tested three different values (5, 10, and 20 years of sampling), and found that indeed, shorter sampling periods resulted in lower clustering rates and longer terminal branches. This is probably because with shorter sampling periods the proportion of samples at the edges of the sampling window is larger. These samples are less likely to be clustered, as only one side of the transmission chain is present in the sample set (Page 3, Ln 108; Appendix 3).

The sampling period could potentially influence the analysis of scenarios with different R0. It is possible that scenarios with R0 < 1 result in generally shorter trees (because the lineage goes extinct within the first 10 years), potentially leading to a difference in the sampling period. Therefore, I repeated the analyses of transmission rates, infectious period, R0, and the example. In these new analyses I added an additional condition for the simulation to be accepted. If the tree height was lower than 10 years, the simulation was rejected. In this way I ensured that all simulations had similar sampling periods (~ 10 years). The results did not change, meaning that the sampling period did not bias the results of these analysis.

Regarding the settings: all simulations are initialized as a single individual, simulations that result in less than 100 tips are discarded. For R0 < 1, a larger number of simulations needs to be discarded, as trees are tendentially smaller.

4) R1 notes that more comments on the sources of data that could help to disentangle the impact of differences in latency vs transmission would be helpful. Presumably that's a population-representative sampling that would enable knowing the fraction of cases that are of the two co-circulating types, over time? Or good estimates of the infectious duration? Or either, or both.

At least in theory, all the parameters of the model can be estimated with phylodynamic analyses. However, this is quite challenging for MTB data sets. In addition to clustering rates and TBL, researchers can use complementary data to understand the epidemiological dynamics. For example, the change through time in the proportion of cases belonging to a certain type can be used as a proxy for the relative reproductive number (i.e., if a type is increasing in proportion within a population, it should have a larger R0 compared to other types). Similarly, if the absolute number of cases is increasing, R0 should be greater than one. However, there can be confounding factors, such as the migration of strains from other regions. Finally, the length of the infectious and latent period can be estimated (with some assumptions) from prevalence surveys and notification data. Usually this is done at the population level, the same analysis on individual lineages, or clones, could provide useful insights on the heterogeneity of MTB populations.

I included this paragraph in the discussion (Page 7, Ln 243-255).

5) Please address minor points on clarity and communication raised by the reviewers. In addition, could you comment on the mutation rate during latency, and what might the result be if mutations were not occurring during latency at the rate that they do in active disease?

I addressed all points raised by the reviewers (see below). Regarding latency and clock rates, first a clarification: here, with latency, I simply refer to the time between being infected and becoming infectious, not to a particular bacterial metabolic state (as often done in the TB literature). In any case, the effect of latency and clock rate are additive. For example: in a type with longer latency and lower clock rate during latency (which would lead to lower clock rate overall), the longer latency would result in lower clustering rates and longer TBL, while the lower clock rate has the opposite effect. The outcome will be determined by the magnitude of the changes: if the variation in clock rates is large enough to compensate for the longer latency, clustering rates will be larger and TBL shorter, otherwise clustering rates will be lower, and TBL longer.

I included this point in the discussion (Page 9, Ln 328-333).

Reviewer #1 (Recommendations for the authors):

Overall, the manuscript was fairly clear and straightforward to interpret, with a few comments below. The main concern I had was with the overarching conclusion that TBLs are non-informative about R0. The results show that transmission rates are associated with TBL/clustering, but the argument is that R0 is not because it depends on the ratio between the transmission rate and removal rate, each of which can alter those metrics. But it would be reasonable to make assumptions about the infectious duration, for example within a lineage and/or location, enabling comparison between variants within a lineage, within lineage over time or space, etc. Holding infectious duration constant, comparisons of R0 would then be possible through TBL or clustering. Moreover, in some cases, the infectious duration can be directly estimated (i.e. by performing a prevalence study and comparing with notifications), enabling the relation of R0 to TBL/clustering directly through transmission rates. These issues should be addressed and/or the conclusions about lack of ability to infer R0 differences should be tempered a bit. This wouldn't diminish the importance of the study, as explaining the relationship between all of these parameters and these metrics is the valuable task that this article undertakes.

I addressed this in the answer to point (1) I report here again the main point.

In empirical studies, TBL and clustering rates could be used to estimate R0 only if the clock rate, the sampling proportion, the sampling period, the length of the infectious period, and the length of the latency period are ALL known or held fixed. If the infectious period is fixed, but the other factors are not, R0 would not necessarily correlate with clustering rates and TBL. To show this more clearly, I modified the example, which now consists of two types with identical infectious periods, but different latency, transmission, and clock rates. One type has R0 = 0.9, and the other has R0 = 1.1 (because of different transmission rates). The shrinking type resulted in higher clustering rates and shorter TBL (because of differences in clock rates and latency). This is described in the response to point (2)

I think some confusion on these issues was generated by the order in which I presented the results. To make these points clearer I modified the structure of the Results section (Pages. 3-5). I also amended the abstract (Ln 21) and expanded the discussion (Page 7, Ln 234-242).

Reviewer #2 (Recommendations for the authors):

First, this is a cool and timely paper, and I really like the approach.

There are however some things I believe could be clarified and perhaps explored a little bit further.

Page 7 "An example": Drawing up an example of how different characteristics of co-circulating TB-types can result in contra-intuitive findings is great. In the example, a less transmissible type nevertheless has an R0 > 1 whereas a more transmissible variant has R0 < 1. Looking at table 2, it seems evident that the more transmissible Type 2 is contracting (R0<1) due to higher cure/death rates and higher sampling compared to the less transmissible Type 1, whereas it exhibits shorter TBLs and higher clustering due to shorter latency and a slower molecular clock compared to Type 1.

I have a few questions/comments which I hope will clear this up a little bit on my part:

- I think the example would be easier to follow if you explained which of the parameters influenced the R0, TBL and clustering rate, as I have tried to do above. If my summary is not correct, I guess that illustrates that this is a bit complex for many readers.

The summary of the original example is correct, there is only one missing aspect: the shorter TBL and larger clustering rates of Type 2 are indeed caused by shorter latency and lower clock rate. But, the larger transmission rate contributed to these trends as well.

In response to points 2 and 9 I changed the settings of the example, which is now simpler. In the updated version, the differences in R0 are due to different transmission rates, because the sampling and death rates are identical between the two types (therefore they have also identical infectious period). The shorter TBL and larger clustering rates of type 2 are now due to its shorter latency and lower clock rate (while the transmission rate is contrasting this trend).

I reworked the explanation of the example, because in the new version the two types are more similar to each other, it should be clearer what factors are influencing R0, TBL and clustering rates (Page 6).

- This is a bit embarrassing, but what specifically should sampling rate be interpreted as here, and throughout the paper? If I get this right, sampling represents both observation and removal (the patient is sampled, and hence also cured and does not transmit further)?

Yes, correct, the sampling rate ε is the rate at which infectious individuals are diagnosed and treated. It is assumed that these individuals stop being infectious immediately, hence they are removed from compartment I.

I added this information to the text (Page 3, Ln 92-94).

- I think it would also be cool with a little description of the two types in more biological/medical lingo. If I understand correctly, Type 2 is more transmissible, and also characterized by rapid onset of disease (short latency), followed by rapid death, self-cure or health-seeking. This actually makes sense, and the only reason R0 < 1 is that health-seeking, cure and death rates are even more elevated than transmissibility?

Yes, this was correct for the original analysis. However, it was pointed by R1 that there is no direct evidence for such divergent sampling and death rates within the same community. This is a fair criticism, I personally think that different durations of the infectious period (together with different transmission rates) could in part explain different TBL and clustering rates between e.g., L1 and L2, but the evidence supporting this is indeed indirect.

In any case, I clarified the differences between the two types in the text (Page 6).

Page 3: If I'm not wrong, the model used is a form of birth-death model? Perhaps this could be spelled out, and if not, explain how it differs from a B-D model.

Yes, it is a BD model, I added this information (Page 3, Ln 88).

Suppl figures: I believe the numbering of the suppl figures in relation to the text is wrong in quite a few instances

Fixed it. Thanks, and sorry for this.

[Editors' note: further revisions were suggested prior to acceptance, as described below.]

The manuscript has been improved but there are some remaining issues that need to be addressed, as outlined below:

In response to the question about bias induced by having to choose among your simulations, I fear that you have misunderstood the question and therefore haven't addressed it. (The new exploration of the impact of sampling period is nice to have, though). You write: "Regarding the settings: all simulations are initialized as a single individual, simulations that result in less than 100 tips are discarded. For R0 < 1, a larger number of simulations needs to be discarded, as trees are tendentially smaller. "

This is a source of potentially large bias because the tree you obtain is likely not to meaningfully *have* an R0 of 0.9. To grow from 1 to 100, the mean number of offspring that is realised is definitely above 1 for an extended period, even if the parameter you set in the process was such that the *expected* number was less than 1. (I would be curious to see a birth-death simulation with an R0 of 0.9 that grows from 1 to 100 tips and last over 10 years). The probability of this growth is not zero, but it's small, so the trees you obtain are not at all representative of R0 = 0.9 trees. This means that when the clustering differs or fails to differ, you can't interpret that as being able to identify (or not) trees from R=0.9 vs R=1.1.

Thank you for clarifying your comment, which indeed I had misunderstood in the first revision. This is an important point, before addressing it (see below), I want to clarify one detail regarding the simulations. To be accepted a simulation does not need to grow to 100 tips existing at the same time. Rather, 100 individuals must have existed over the whole simulated period. Therefore, the simulated lineages are not necessarily growing strongly, in fact some may go extinct before the end of the simulated period. These settings represent a much-reduced source of potential bias compared to what is written in the comment above.

I clarified the description of the simulation setup to avoid misunderstandings (page 10):

“A post-simulation condition was implemented by specifying the minimum and maximum number of tips in the phylogenetic tree necessary to accept the simulation. This threshold regards the total number of tips sampled during the simulated period, not the number of lineages existing at the most recent time point.”

The tree sizes will be different too, and from coalescent theory, we know that this alone will change the patterns of genetic diversity (and hence the clustering). I would hazard that if you took simulations from identical parameters (birth, death, sampling, etc all the same) but in one set you chose the rejection criterion to have the trees need to grow bigger or last longer than in the other set, you would also find differences in clustering as a sole function of the cutoff criterion.

I partially investigated this already in the first version of this manuscript (see Results at pages 3-4, and Appendix 4), although only for R0 = 1. In that analysis I found that different cutoff values on the minimum tree size did not change clustering rates or TBL. I now expanded this analysis and performed it also for different values of R0 confirming that the results are robust to different thresholds on the minimum tree size (see answer to next point).

Some portion of what you find is due probably due primarily to size and you could explore this, as well as not so dramatically biasing the trees you analyse by requiring extremely unlikely events before a tree gets into your sample.

Think of it the other way. If you took trees from a simulation with R0> 1 but you *only* looked at trees that went extinct before reaching 10 taxa, and then you did estimation or some other analysis on those trees, they would look very similar to trees with an R0 < 1 because those also die out. By removing the ones that grew you removed the information that R0 was set to be > 1. Conversely, here, by rejecting the vast majority of trees simulated under R0 < 1 that died out, you are removing the information that R0 was ever < 1.

I expanded the analyses investigating transmission and sampling rates (Appendix 5). I focused on these because they are the ones exploring the broadest range of R0 values (from 0.8 to 1.2). With the original settings I accepted only simulations that resulted in at least 100 tips over the whole simulated period. I repeated these analyses using 50 and 25 as alternative minimum tips thresholds, with these settings the potential bias should be reduced. In addition, I included a scenario with a larger threshold (200 tips), to observe how the results changed when increasing the potential source of bias.

For both transmission and sampling rates I confirmed the results of the original analysis for all thresholds. Increasing transmission rates (and therefore R0) resulted in larger clustering rates and shorter terminal branches. This trend was clear for all thresholds, and it was monotone over the entire range of R0 values (0.8-1.2). Conversely, the sampling rate did not strongly affect clustering rates and TBL. There might be a weak trend towards lower clustering rates for lower sampling rates, this was similar for all thresholds, and it was mostly visible for R0 values between 0.8 and 1 (Appendix 5 – figures 1-8, Appendix 5 – tables 1-8).

To explore the influence of different thresholds more directly, I compared the results obtained with the same settings and different thresholds. To limit the number of supplementary tables and figures I focused on the scenarios with R0 = 0.8, 1, and 1.2. When I compared the results obtained for R0 = 0.8, the scenario with potentially the largest bias, I found that the clustering rates and TBL were not identical when using different thresholds on the minimum tree size. There was a weak trend towards lower clustering rates for larger thresholds (on the minimum tree size), although this was not monotone for all SNP thresholds (Appendix 5 – figures 9-10, Appendix 5 – table 9-10). Additionally, for one of the analyses with R0 = 0.8 (λ = 1 and ε = 1.25) I found that larger thresholds also resulted in longer terminal branches, although the difference between the mean of the TBL distributions was smaller than 0.01 SNPs (0.3729, 0.3768, 0.3773, and 0.3815 for a threshold of 25, 50, 100, and 200 respectively). These trends were not present for larger values of R0, indicating that they are indeed a bias due to the minimum tree size, which affect scenarios with shrinking populations (Appendix 5 – figures 11-13, Appendix 5 – tables 11-13).

These results suggest that the threshold on the minimum tree size can bias the outcome of some analyses. However, with the settings used in this study this bias was negligible, and the results were robust to different thresholds. For the transmission rate, the trend towards larger clustering rates for higher transmission rates was clear also for scenarios with R0 equal or greater than one, these scenarios are scarcely affected by the bias, further indicating that this is a genuine trend due to the epidemiological process (Appendix 5 – figures 1 and 3-5). Conversely, for the sampling rate the observed trend was evident mostly for scenarios with R0 smaller than one, suggesting that it was caused by the threshold on the tree size (Appendix 5 – figures 2 and 6-8). In any case the differences were rather small, and I conclude that if there is any effect of the sampling rate on clustering and TBL, this is extremely weak.

I report all these new analyses in Appendix 5 (page 21). I also modified the manuscript (pages 4-5) to include the latest results.

I have asked for a revision because I think this point is not minor, as quite a bit of the paper is focused on clustering and R0 under simulations, and this same issue could impact many of them.

I agree that these aspects needed to be explored before publication. I would also like to point out that there are several factors that are not related to transmission, or R0, but influence the results of clustering and TBL (the most relevant are latency and the molecular clock rate). Therefore, this study would reach the same conclusion (i.e., that clustering rates and TBL cannot be used to identify differences in transmission or evolutionary success, unless all other factors are known or identical between sub-populations), also without including the analyses of scenarios with R0 < 1. Of course, it is much better to include these analyses, as they allow to explore the effect of the transmission and sampling rates, and therefore R0, on clustering rates and TBL, and I’m glad I could show that they are robust to different thresholds on the minimum tree size.

It would be interesting to explore a genuinely declining population (rather than one that is declining in expectation because R0 < 1 but growing in its actual realisation because you start with 1 taxon and branch, rejecting those simulations that do not grow). However, initialising that population is challenging because the results would depend on the initial genetic diversity.

The editor already identified a fundamental challenge. In a stable population (R0 = 1), 50% of the lineages existing at any time point are expected to die without transmitting further, with R0 = 0.8 this proportion rises to 60%. Therefore, if we initialize a population and consider all individuals, the clustering rate will be in large part determined by the initial genetic diversity, which has nothing to do with the simulated epidemiological processes. It would be possible to repeat this analysis several times with different “starting populations”, characterized by different underlying genealogies. However, the set of possible starting conditions is extremely large, and I believe that these analyses could not alter the main findings of this study.

As a side note – under point 4 in the response, you note that "At least, in theory, all the parameters of the model can be estimated with phylodynamic analyses." But this paper finds a fundamental unidentifiability that suggests that these parameters are not identifiable: https://academic.oup.com/mbe/article/38/9/4010/6278301. They use likelihoods, not clustering, but since the likelihoods and portion clustered are fundamentally based on the branching times in the phylogenies it seems that their results would probably carry over.

Thank you for pointing this out. I removed that sentence from the manuscript, and I amended the conclusions (pages 9-10), to take into account these findings:

“Phylodynamic methods that estimate the parameters of an epidemiological model from genomic data are becoming available (Kühnert et al. 2016, Didelot et al. 2017, Volz and Siveroni 2018), however, these analyses have limits that cannot be overcome exclusively with phylogenetic data (Louca et al. 2021), and in any case they are challenging with current MTB data sets (Kühnert et al. 2018, Walter et al. 2022). Methodological advances, the integration of different types of data, and more complete and longer sampling series of MTB epidemics will allow to study epidemiological dynamics more accurately in the future. In the meantime, the results of clustering and TBL analyses should not be over-interpreted.”

Author details

1. Fabrizio Menardo

   Department of Plant and Microbial Biology, University of Zurich, Zurich, Switzerland

   Contribution

   Conceptualization, Formal analysis, Funding acquisition, Investigation, Methodology, Software, Visualization, Writing - original draft, Writing - review and editing

   For correspondence

   fabrizio.menardo@uzh.ch

   Competing interests

   No competing interests declared

   "This ORCID iD identifies the author of this article:" 0000-0002-7885-4482

• 1,003

  Page views

• 147

  Downloads

• 6

  Citations

Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.