Apologies to the authors for my delayed response.

This previous substantive concern has not been addressed: "The main conclusion of the paper is that the epidemic in new hosts is growing faster than that in classical hosts, however the confidence with which this conclusion can be made is not stated. The doubling time in classical hosts since 1997 is estimated to be 0.58 - 10.13 years, and the doubling time in new hosts 0 - 3.51 years. The relevant quantity for the conclusion is the posterior for the ratio of these two parameters. The authors do present a ratio comparing the two hosts with regards to the reproduction number, and find that the confidence intervals do not exclude 1; if the same is true of the doubling time, which seems plausible given the similar parameter dependencies of R_0 and doubling time, the main conclusion is not supported. The same point applies to the other host parameters inferred to be different: assortativity and recovery/removal rate. The parameters themselves need not be redefined, but the posteriors of their ratios should be examined to support claims of differences" (and the precision with which these differences have been measured). This applies both to the difference between the two host types, and also to the difference in the effective reproduction number before and after the third generation tests, where a difference is currently presented despite strongly overlapping confidence intervals in the two separate values. If the ratio of doubling times does indeed have a credibility interval excluding 1, I would suggest reporting it as 10 [x-y] times smaller for new hosts instead of an order of magnitude lower, to also communicate the precision of the estimate.

I am uncertain, but think my concern about differences in sampling has not been addressed. My concern was that higher sampling in the new group than in the classical group might lead to the data observed, and to the resulting conclusion, even if there were no difference in transmission dynamics in the two groups. I think the authors' reply can be summarised as saying that if the new group was under-sampled, spread there would be even faster than currently estimated. But the same is also true in the classical group, and if it were true there to a greater extent - i.e. lower sampling there - then the true difference in transmission dynamics between the groups would be less than that estimated here. In the absence of any information of the true underlying population sizes, the most defensible assumption on sampling might be that equal proportions of both populations end up diagnosed, and not that equal proportions of both populations end up included in the present study, given that the classical group seems to have been deliberately downsampled and by a known factor.

Unless I have misunderstood, the authors' choice of the prior Unif(1, 10) for the parameter nu means it was impossible to come to any conclusion other than a greater transmission rate for new hosts. The prior for nu should have equal weight above and below 1 if we are to learn how the data inform the ratio of transmission rates; here it has zero weight below 1. I think that treating the two hosts equally a priori (i.e. imposing a host-type interchangeability symmetry) implies that the prior for log(nu) should be symmetrical about 0. Thank you for the suggestion. We now run the same ABC analyses with a prior of \nu in Unif(0,10) to include both hypotheses. Re-articulating my previous concern, this was that 100% of the prior probability mass was for nu greater than 1, when it should be 50% if the prior is to be uninformative of differences between the host types. In modifying the prior from Unif(1,10) to Unif(0,10), it is now the case that 90% of the prior probability mass is greater than 1; this is still far from 50%. More specifically, for parameter priors to be agnostic about the hosts, they should be unchanged by relabelling the two host types new <-> classical. For parameters describing a multiplicative difference, the prior p should be such that it is equally likely that the one group's parameter is x times larger as it that the other group's parameter is x larger. In other words p should satisfy int{1/a}^{1/b} p(x) dx = int{b}^{a} p(x) dx for any a and b. Equivalently, probability mass should be symmetric about zero on a logarithmic scale for x, i.e. F(log(x)) dlog(x) = F(log(x))/x dx for any symmetric function F. The simplest example is F=1, i.e. a prior of 1/x, truncated above a point L and below a point 1/L to give a proper prior. Given the previously plotted shape of the posterior for nu, I am confident that the authors' qualitative result will be unchanged with this more appropriate prior. However the paper's conclusion would be considerably more persuasive if it arises from an uninformative prior than from a prior geared 90% towards that conclusion. Also judging by the posteriors in Figure 2, the aforementioned truncation point L for the parameter nu looks like it ought to be greater than 10; similarly, priors allowing for higher values of gamma2 and R01^{t2} seem justified given the data (posteriors squashing to the prior boundary).

Minor suggestions:

For a number of points where I requested clarification, this was provided in the response without indicating whether it was also added to the manuscript. For the authors to check and consider, to taste.

With the shift in nomenclature from basic reproductive number to effective reproductive number, it would be appropriate to drop the subscript 0 in R_0.

146 onwards: I don't understand the hypothesis being tested here. I think it is that previously the epidemic consisted only of people who were diagnosed a long time after infection, but more recently there is a separate epidemic of people who were diagnosed shortly after infection. This only makes sense to me if diagnosis time is merely a proxy for another characteristic that is causal for difference in dynamics. Here, we wanted to validate the fact that the differences in the structure and the labels in the phylogeny are not due to the infection phase (acute vs. chronic), but rather to the epidemiological profile (classical hosts vs new hosts). This is because new/MSM hosts are all detected during acute infections. I think I now understand this: in the alternative model, nu is the difference in infectiousness for the same person in acute and chronic phases, not for two different people who tend to be diagnosed in different stages, correct? Including the ODE for this model in the supplementary material would clarify this.

80: "posterior to 1997" -> "after 1997"

108: "To better apprehend" -> "To better understand" (or similar)

113-116: consistent time units throughout would aid comparison of these numbers.

188-191: I suggest moving this sentence to the introduction. Specifying precisely how the two host groups were defined is helpful for following the rest of the paper (c.f. our earlier confusion).

261 (just before): I think the numerator should be (clarified to) "number of descendant leaves labelled Y"

Figure 2: I suggest spreading these plots over two rows, as the horizontal squashing makes them hard to read. Also using the same x axis scale for the two reproductive numbers would aid comparison.

Best wishes, Chris Wymant

The authors have improved the manuscript and addressed some of the comments from the previous round of reviews, but there are still a few issues with the manuscript.

Major comments

1. How representative are the TMRCAs of the 10 replicate trees of the posterior distribution estimated for the TMRCA is BEAST2? This should be shown. It is of necessity a very small number of replicates (10 samples are usually not enough to sufficiently represent a distribution), and if they are all very close to each other, then this doesn't address the underlying issue - that the target tree's TMRCA distribution has an HPD interval of ~35 years, which is reduced to a single median estimate in the regression-ABC analysis.

2. In the previous review I suggested that including a background set of global HCV genotype 1a sequences sampled over a larger time period will help with reducing the TMRCA uncertainty (adding more sequences over a longer time period means that internal node height estimates are more accurate). The target tree can then be pruned from this bigger tree. Alternatively, I would have to be convinced that the 10 replicates are representative of the posterior TMRCA distribution of the BEAST2 analysis (point 1).

3. Line 140: I agree that if new hosts are undersampled then the analysis is conservative. However, it seems more likely to me that new hosts are oversampled. If the sampling proportion is much higher for new hosts, and high assortativity holds, then we would expect sequences from new hosts to coalesce faster, even if both populations are growing at the same rate (since the population size of new hosts would be smaller). I think the phylogeny is consistent with both the interpretation in the text (equal sampling rates and different growth rates) and with very different sampling rates and equal growth rates.

Minor comments

1. The statement about BEAST2 being used to root a PhyML tree is still in the legend to Figure 1.

2. tf in the methods should be t3

3. Even with a fixed strict clock the target tree TMRCA is still very uncertain with an HPD of ~35 years (very much the same as when co-estimating both parameters). Could the authors please double-check that? Since the TMRCA and the clock-rate usually have a very strong negative correlation to each other, fixing one should also restrict the other. This could mean that the sequences are not very informative about the branch lengths (in substitutions/site) or topology. Perhaps the authors could check if very few nodes in the MCC tree have a high posterior probability (or alternatively look at bootstrap scores in the ML tree).

4. In the coalescent simulation, are branching times drawn from exactly those transmission events in the simulated epidemiological trajectories (step one in the simulation procedure) or are they simply drawn from the probability of two lineages to coalesce given the population sizes from the simulated trajectories?

5. As explained by the editor, a multi-type birth death model is sufficient to model the scenario. However, the revised manuscript still contains a paragraph stating that birth death models are not applicable. Two types of state-change events can be modelled in BEAST2 with BDMM, the default is indeed mutation based, but in the other states change on a transmission event (exactly at branching times).

Summary

This manuscript investigates the roles of different populations in sustaining ongoing HCV epidemics in France using 213 newly sequenced HCV genomes and an ABC phylodynamics approach to infer epidemiological parameters. In particular the authors examine the difference between so-called classical hosts, who spread the disease through needle-sharing and new hosts that spread HCV through sexual transmission.

I found some disconnect between the populations described in the introduction and the ones actually used in the analysis. This may just be a small misunderstanding on my part, but to ensure that everyone is on the same page I wrote down a detailed description of my understanding below.

The weak point of the ABC approach is that all parameter estimates (and some of the model priors) are conditioned on a single independently estimated tree. If this tree is poorly estimated, parameter estimates will likely be biased. From what is presented in the manuscript I am not convinced that the dataset used here has sufficient temporal resolution for the reconstructed tree (in this case the MCC tree from a BEAST analysis) to be correlated to the transmission tree. I also expand on this point below.

I think the application is of high interest and should definitely be investigated. The ABC method is interesting as an alternative to the usual Bayesian phylodynamics approach, since it allows using (almost) arbitrarily complex epidemiological models without the need to calculate the probability of observing a tree under the model. However, it also comes with its own set of limitations, and in this case I think the dataset may not be suitable. In the comments below I offer a few suggestions for addressing this issue. As things stand I trust the estimates of the ABC method, given the tree the authors used, but I don't trust that tree and so I would not extrapolate the results to the real world.

Major comments

1. The different populations of interest and study populations should be presented more clearly.

   • As I understand it, there are 4 populations of interest:
     1. HIV-positive PWID
     2. HIV-positive MSM
     3. HIV-negative MSM using PrEP
     4. HIV-negative MSM not using PrEP

   • (There is some degree of overlap between populations 1 and 2 and population 4 is mostly unsampled).

   • By the definition in the introduction classical hosts predominantly stem from population 1 and new hosts from the remaining 3 populations, where transmission is sexual.

   • In the genetic dataset new hosts are represented by all MSM patients (which presumably includes representatives from populations 2,3 and 4). Classical hosts are represented by non-MSM, HIV-negative male patients.

   • Did I understand all of the populations correctly?

   • The relationship between classical hosts in the genetic dataset and HIV-positive PWID (population 1) is not clear to me. Are they also PWID? Do you have any information about the main mode of transmission in this population or a reason to believe that it is through needle-sharing?

2. The main informative events for the ABC approach are the branching times in the tree. Since the ABC approach conditions on a single tree, it is extremely important that the timing of the tree is accurate. Judging from the TMRCA estimate, which has an HPD from 1962 to 1997, I think it's likely that many of the other internal nodes have a similar range of uncertainty, which is ignored when conditioning on the MCC tree. The result is that estimates that rely on the timescale (reproduction numbers, infectious periods, origin times) are likely to be biased. As an example, the authors placed a uniform prior from 1962 to 1997 on t0. However, by conditioning on the MCC tree, the upper bound of this prior implicitly becomes 1981, which explains why no t0 estimates are larger than 1981. In addition, the authors divided the tree into 3 epochs, with different reproduction numbers. However, it is likely that the uncertainty associated with internal nodes stretch across multiple epochs, which call into question the estimates reported here. (If the procedure was repeated with different trees drawn from the set of posterior trees parameter estimates would be different).

   • At the least the authors should test for a temporal signal in the data. Since sequences are only sampled over 4 years and HCV has highly variable within-host evolutionary rates, it is possible that there isn't a good clock signal. The authors should also provide more details of the model used for the BEAST analysis, as using a poorly fitting tree prior can lead to inferring biased branching times.

   • Even if there is a temporal signal in the data it may still be impossible to obtain accurate estimates of the branching times, since the sampling period is small compared to the tree height. In this case, the analysis could be performed with a fixed or highy constrained clock rate (based on previous analyses of HCV genotype 1a). Alternatively, more genotype 1a sequences from Europe (spanning a bigger sampling period) could be downloaded from Genbank. This will allow internal nodes to be more accurately estimated. If the Lyon sequences form a monophyletic clade, this clade could be pruned from the tree and used for further analyses by the ABC approach.

3. I believe the authors are estimating the effective reproduction number (Re), not the basic reproduction number (R0).

4. Figure 1 should be improved. It is impossible to see the structure of the tree close to the present or to read the tip labels.

5. Line 213: \nu is not simply a factor of the [mean] number of partners of classical (I1) and new (I2) hosts. The modes of transmission differ between I1 and I2 (blood vs. sexual) and there is also presumably a difference in the per-act likelihood of transmission between the different modes of transmission.
   In addition, \nu is constrained to always be greater than or equal to 1. This appears to be backed up by the data, but is not necessarily always the case. Why is it a priori assumed that there is a higher transmission rate between new hosts (MSM) than between classical hosts (PWID)?

Minor comments

1. What is meant by the phylogeny was estimated in PhyML and then rooted using BEAST? Does this mean that you performed a BEAST analysis on the sequence data using the maximum-likelihood tree as a starting tree?

2. Line 57: I don't know the shape of the posterior distribution for the TMRCA, but regardless of how peaked it is, I wouldn't place much confidence on the TMRCA being in the early 1980s when the HPD interval stretches from 1962 to 1997.

3. Do the authors use the epidemiological data from the Dat'AIDS cohort anywhere? Is there a relation between the cohort and the sequences sampled in Lyon?

4. The AIC is not a good method for Bayesian model comparison. To choose between a strict and relaxed clock it makes more sense to check if the HPD for the coefficient of variation (for the relaxed clock model) excludes 0.

5. Line 242: Please provide more details about how coalescent trees are drawn from the simulated trajectories and sampling dates (are the sampling dates fixed to the truth or are they also from the model?).

6. Using the word "cluster" for heterogeneous and homogeneous clades is misleading as cluster has a specific interpretation in an epidemiological setting. Clade is a better term to use here.

7. Heterogeneous clusters (should be clade) of type Y are defined as clades where more than 70% of leaves are of type Y (and includes clades where all but one leaf are of type Y). What about very heterogeneous clades, where type Y makes up between 50 and 70% of leaves? Are these clades completely ignored? In addition, did you use a greedy algorithm to find the biggest possible clades or do you have nested homogeneous and heterogeneous clades? (the same branch can appear in multiple clades).