I am grateful to the authors for taking my previous comments seriously. I think that the manuscript is much improved and remains interesting. Nevertheless, looking at the new Figure S7, I still have concerns about the legitimacy of estimating Ne from S_bar. If the authors agree, I think that the concerns should be emphasized more strongly. I also think that the authors could do more to show that the results could in principle be real. They need to show e.g. that a past bottleneck could lead to very large differences in the Ne values that apply to different statistics. Neither of these suggestions would require large changes.

1. Can Sbar be used to estimate Ne when the ri vary?

In the very welcome new section 460-474, the authors come close to stating that their key statistic, Sbar has no clear meaning.
They argue that, with highly variable ri, “what exactly \theta measures in this case is unclear”. I think \theta does have a clear meaning: the mutation rate multiplied by the length of the terminal branches in the genealogies. But the meaning of S_bar really is unclear, because it affects predictions for all frequency classes.

The authors then argue that the realized ~3-fold variation in the ri estimates (0.5-1.5; Figure S7) "does not suggest to us that the ri's pose a major problem of comparability in this analysis." This statement seemed too confident to me. Because the ri do not apply to different epochs, but to branches across the entirely genealogy, I think that the 3-fold differences imply quite major departures from a standard neutral genealogy, and it is not clear how Sbar will behave in those cases. For example, I think Figure S7 implies lots of high frequency polymorphisms in the flies, which is suggestive of balancing selection or structure.

Currently, it is suggested that the reader can skip sections 1-3 of the discussion, but, unless the authors think my comments above are mistaken, I think that sections 4-5 need to clearly acknowledge the possibility that S_bar does not provide a meaningful estimate of Ne.

Smaller comments/suggestions:

1. The simulations.
   The simulations should prove to the reader that non-equilibrium demography could lead to very large differences in the range Ne values as estimated from piS and Sbar. Currently, this is not very explicit from Figure 4. Also, I could not find a description of the bottleneck depth at generation 15,000.

2. The introduction
   The introduction should probably acknowledge the classical results showing that populations can be characterized by different values of Ne when the Wright-Fisher assumptions are violated (e.g “inbreeding effective size” vs “variance effective size” etc.). These differences are not very surprising in themselves. What is surprising is the very large (~10-fold) differences estimated in this work.

3. “To our knowledge, the Gamma + lethal model has been tried in two studies before this one.” Didn’t Nielsen and Yang 2003 MBE also use this model?

4. In my previous comment 2.3, I asked whether the authors might attempt to estimate s_bar (i.e the mean selection coefficient unscaled by Ne). I think that previous authors have done this, and it relates to the interesting discussion in lines 631-635.

5. What are the CIs in Figure S7?

6. Is Kaiser and Charlesworth 2009 TIG relevant to the discussion of linked selection and mutation load? (lines 621-625).

7. l445: “can be interpreted as the relative effective mutation rate”. Could r_i be more clearly described as the relative length of the portions of the genealogies that lead to i tips?

8. There are two different legends for Figure 4 (pp. 35 and 39).

Grammar etc.
l38: "which implies to also infer" should be "which implies also inferring".
l51: "distribution of fitness effect[s]".
l72: "can as well be" should be "can also be".
l101: "by jointly analysis" should be "by jointly analyzing".
l242 and Figure 4. Does 2.10^4 mean 20,000? If so, notation could be clearer.
l243: 2.2 10^{-7} (missing multiplication sign).
l269: "we rather fitted" should be "we instead fitted"
l303 and l310: "probability to be observed" should be "probability of being observed"
l331: "Setting plth [set]".
l338: "colons" should be "columns".
l440: “parameters of nuisance” should be “nuisance parameters”.
l464: “in [the] absence of”.
l536: “implies to accommodate” should be “implies accommodating”.
l556: “did not change much the picture” should be “did not change the picture much”.

The effective size of a population is the census size of a Wright-Fisher population that would give the same value of some statistic, related to drift. When the Wright-Fisher assumptions are violated, different statistics can imply different values of Ne. Differences between different methods of estimating Ne (e.g. from mean time to pairwise coalescence, vs. total genealogy length) are interesting, because they tell us about violations of the Wright-Fisher assumptions, and because the different Ne values might affect different things (e.g. total neutral diversity vs. current efficacy of purifying selection).

This preprint compares the extent of differences in Ne, estimated in two ways. First, it uses piS, the mean pairwise differences at synonymous sites, where E(piS)=4Nemu. Second, it estimates S=4Nes using the non-synonymous and synonymous SFS. Estimates are obtained for a wide variety of animal species, but focussing on primates vs. Drosophila. Results show that S/min(S) is 1-2 orders of magnitude more variable than piS/min(piS). This is a surprising result, especially if the DFE is identical between species.

The authors relate their finding to Lewontin's paradox: the finding that genetic diversity between species varies much less than would be expected under an equilibrium neutral model, assuming that effective population sizes are proportional to current census sizes. The authors implicate non-equilibrium demography, because piS is strongly influenced by past periods of low population size, while piN/pi_S equilibrates more rapidly, and so should be less influenced by past demography. This point is illustrated with simulations.

The preprint presents interesting work on an important topic, but I think some major points need to be clarified and developed before the results can be properly understood.

1- What is being estimated?

This preprint compares two different types of Ne, but I think the authors could be much clearer about what they measure.
In the equation E(pi_S)=4Nemu, 2Ne describes the average time to coalescence of two randomly chosen sequences.
The method of Galtier (2016) estimates Ne in two different ways, via the compound parameters theta=4Nemu and S=4Ne*s. Only the second measure is reported, and I could not tell when/whether the two Ne values are expected to differ from each other, and whether the second measure (which is the only one reported) is a valid measure of Ne.

From eqs. 6 and 10 of Galtier 2016 with r1=1, it looks like 4Ne estimated from theta/mu measures the mean length of the terminal branches on the genealogy. With non-equilibrium demography, this could differ from the mean time to pairwise coalescence. The lengths of the internal branches are fit via the nuisance parameters (ri).

I found the other Ne (estimated from S/4s) much more difficult to interpret. S is estimated from the complete SFS of non-synonymous mutations. But it does not use existing results for the SFS with selection and non-equilibrium demography or linkage effects (e.g. Evans et al 2007 TPB; Good et al 2014 PLoS Genet). Instead, it uses the standard equilibrium formula for the population allele frequency (eq. 4 of Galtier 2016), while also allowing for non-equilibrium effects on the shape of the underlying genealogy, by allowing for arbitrary variation in the r_i.

For this reason, I struggled to understand what these Ne estimates meant, and whether all of the possible variation in the non-synonymous SFS could be modelled via variation in Ne, even assuming that the gamma+method DFE is accurate.

Because the variation in these Ne values is the major result of the paper, I think this needs to be clarified.

2- Model adequacy

Estimates of S are highly model-dependent, and the authors present some suggestive evidence that the standard gamma model underfits the SFS data. They use a gamma+lethal model and show that the extra parameter has a substantial influence on the estimates.

This is an interesting result, but I think that the authors might be too quick to assume that the gamma+lethal model solves the problem of model adequacy, and to conclude that there is little difference in the DFE across their data sets.

Previous authors have examined the adequacy of the gamma model, with different sorts of data, including Nielsen & Yang (2003, MBE) and Loewe and Charlesworth (2006, Biol Lett). Nielsen and Yang also used a gamma+lethal model (and a normal + lethal model), while Loewe and Charlesworth argued for a lognormal model. Other authors have combined a gamma distribution with a class of purely neutral mutations (Loewe et al 2006, Genetics; Betancourt et al. 2012 Evolution). It would be useful to know if results are robust to using other distributions like this. While Table 2 contains many good robustness analyses, I think more formal work on model adequacy could be presented to make the headline results really convincing.

Second, the authors argue that their estimated likelihood surface "suggests that the DFE perhaps does not differ so dramatically between primates and fruit flies", but doesn't this flat likelihood surface suggest instead that the parameters are non-identifiable with data of this kind?

Third, as the authors say, the results for Ne depend heavily on the very strong assumption that s_bar has the same value for one set of genes in Drosophila and a different set of genes in primates. Chen et al. 2017 and Loewe et al. 2006 present methods for estimating the strength of selection directly. Could these be used to test this strong assumption?

3- Simulations

Figure 4 presents simulation runs, aiming to show that piN/piS is less affected by demographic events than pi_S. This is an important part of the paper, but I was not sure how well the simulations related to the results reported.

First, the S estimates used the full SFS for non-synonymous and synonymous polymorphisms, while the simulations report piN/piS. These quantities might equilibrate differently (as is evident from Tajima's D).

Second, if I have understood correctly, the "recovery phase" shows a gain in genetic diversity starting from a large but genetically uniform population. Is this realistic? Why did the authors not explicitly simulate the recovery from a bottleneck?

Third, to place their results in context, the authors cite Brandvain and Wright 2016 regarding the different equilibration times of piN and piS. But I think the explanation of the simulation results might be that ratios of pairwise diversity equilibrate more rapidly than raw diversity measurements. This is seen in neutral Fst, for example (Pannell and Charlesworth 2002).

Finally, if the authors suspect that non-equilibrium effects explain some of their results, why do they not test for these effects directly in their data (e.g. by reporting the r_i from Galtier 2016, or Tajima's D for synonymous sites etc.)?

Minor suggestions:

1. The gamma shape parameter is defined in different ways, so it might be helpful to include an equation.

2. The simulation methods state that the negative gamma distribution has a mean of -2.5. Does this apply to 4Ns, and if so, what happens when N changes?

3. data sets with few polymorphisms, and all multi-allelic sites were excluded. Are the authors confident that this does not introduce biases?

4. The heatplots would be clearer with scale bars.