On the relation between the asymptotics of a Dirichlet series' coefficients and the series' analytic continuability There is a wonderful series of articles by Flajolet et. al. about Mellin Transforms and the asymptotic analysis of generating functions. In particular, on page 45 of the article Mellin Transforms and Asymptotics: Harmonic Sums, they state the following result:

Proposition 6 (Growth of Special Dirichlet Series): Let $\lambda_{k}$
  and $\mu_{k}$ admit asymptotic expansions in descending powers of $k$
  as
  $$\lambda_{k}\sim\sum_{r=0}^{\infty}\frac{a_{r}}{k^{\alpha_{r}}}$$
$$\mu_{k}\sim k^{w}\left(1+\sum_{r=1}^{\infty}\frac{b_{r}}{k^{\beta_{r}}}\right)$$
  Then the Dirichlet series $\sum_{k}\lambda_{k}\mu_{k}^{-s}$ can be
  continued to a meromorphic function $\Lambda\left(s\right)$ in the
  whole of the complex plane.

Now, let $V$ be an arbitrary set of infinitely many positive integers, and let: $$\zeta_{V}\left(s\right)\overset{\textrm{def}}{=}\sum_{v\in V}\frac{1}{v^{s}}$$
Enumerating the elements of $V$ in increasing order as $v_{1},v_{2},\ldots$, recall that one way of defining the natural density of $V$ (denoted $d\left(V\right))$ is: $$d\left(V\right)=\lim_{n\rightarrow\infty}\frac{n}{v_{n}}$$ In the case where $d\left(V\right)$ exists and is positive, this gives the asymptotic $d\left(V\right)v_{n}\sim n$. As such, for: $$\frac{\zeta_{V}\left(s\right)}{\left(d\left(V\right)\right)^{s}}=\sum_{n=1}^{\infty}\frac{1}{\left(d\left(V\right)v_{n}\right)^{s}}$$ I have that $$\lambda_{k}=1$$ and $$\mu_{k}=d\left(V\right)v_{k}\sim k$$ which is obtained by taking:$$\mu_{k}\sim k^{w}\left(1+\sum_{r=1}^{\infty}\frac{b_{r}}{k^{\beta_{r}}}\right)$$ setting $w=1$, and letting all the $b_{r}$s be $0$. Hence, unless I am mistaken, Proposition 6 implies that $\zeta_{V}\left(s\right)$ extends to a meromorphic function on $\mathbb{C}$ whenever $V$ has positive natural density.
However, consider the following. Let $\mathbb{P}$ denote the set of prime numbers. Then, $\zeta_{\mathbb{P}}\left(s\right)$ is the so-called Prime Zeta Function, which is known to have a natural boundary on the imaginary axis. On the other hand, since $d\left(\mathbb{P}\right)=0$, it follows that $\mathbb{N}/\mathbb{P}$ has a well-defined natural density of $1$, and thus, by Proposition 6, that $\zeta_{\mathbb{N}/\mathbb{P}}\left(s\right)$ is meromorphic on $\mathbb{C}$. However, since I can write:$$\zeta_{\mathbb{P}}\left(s\right)=\zeta\left(s\right)-\zeta_{\mathbb{N}\backslash\mathbb{P}}\left(s\right)$$ it follows that $\zeta_{\mathbb{P}}\left(s\right)$ is the difference of two meromorphic functions, which forces $\zeta_{\mathbb{P}}\left(s\right)$ to be meromorphic on $\mathbb{C}$, which is obviously not correct. 
So, where's the error, and how (if at all) can it be rectified? In particular, when, if ever, does the existence of $d\left(V\right)$ imply that $\zeta_{V}\left(s\right)$ is meromorphic on $\mathbb{C}$?
 A: Proposition 6 does not imply that if $V$ has positive natural density then $\zeta_V(s)$ extends to a meromorphic function. This is because the density assumption is much weaker than the assumptions in Proposition 6. Indeed, if the elements of $V$ are $v_1<v_2<\dotsb$, then the density assumption says that $v_k$ is asymptotically a constant times $k$, while Proposition 6 requires a full asymptotic expansion of $v_k$ in descending powers of $k$. The latter means that, for any nonnegative integer $R$, there are exponents $w_1>w_2>\dotsb>w_R$ and coefficients $c_1,c_2,\dotsc,c_R\in\mathbb{R}$ such that
$$v_k=c_1 k^{w_1}+c_2 k^{w_2}+\dotsb+c_{R-1}k^{w_{R-1}}+(c_R+o(1))k^{w_R}.$$
The density assumption implies that we have such a relation for $R=1$, namely $w_1=1$ and $c_1=1/d(V)$ work, but it does not imply the relation for $R=2$. And in fact, for the sequence of non-primes $\mathbb{N}\setminus\mathbb{P}$, the above approximation only holds for $R=1$.
A: 
$$F(s)=\sum_{k}\lambda_{k}\mu_{k}^{-s}$$ extends meromorphically with finitely many poles in each strip because $$\Gamma(s)F(s)= \int_0^\infty f(t)t^{s-1}dt, \qquad f(t)=\sum_{k}\lambda_{k} e^{-t \mu_{k}}$$ where $f$ has an expansion in powers of $t$ near $0$ which makes clear the exact conditions you need for it to hold. 



*

*The expansion $\mu_k= k^{w}\left(1+\sum_{r\le R}b_{r}k^{-\beta_r} + o(k^{-\beta_R})\right)$  leads to (expanding the $\exp$) $$e^{-\mu_k t}= e^{-k^wt} \exp(-  \sum_{r\le R} b_r k^{(w-\beta_r) t}) (1+o(k^{-\beta_R})) \\= e^{-k^wt} \sum_{j\le J} c_j t^{d_j} k^{e_j}+ o(t^{d_J}k^{-e_j} e^{-k^wt})$$ 

*together with the expansion $\lambda_k = \sum_{r \le R} a_r k^{-\alpha_r} + o(k^{-\alpha_R})$ it leads to $$f(t) = \sum_{l\le L} C_l t^{D_l} \sum_k k^{E_l} e^{-k^wt}+ o(t^{d_L} \sum_k k^{e_L} e^{-k^wt})$$
For every $r$, for $L$ large enough we have $t^{d_L} \sum_k k^{e_L} e^{-k^wt} = o(t^r)$ so the error term won't be a problem. 
And from our knowledge of $\Gamma(s) \zeta(Bs+C)$, of Mellin inversion and the residue and Tauberian theorem we know that  $\sum_k k^{E_l} e^{-k^wt}$ has an expansion in powers of $t$.
Thus so does $f$ $$f(t) = \sum_{m \le M} u_m t^{v_m} + o(t^{v_M}), \qquad v_m \to \infty$$ from which we have our meromorphic continuation to $\Re(s) > -v_M$ $$\Gamma(s)F(s)= \int_0^\infty (f(t)-\sum_{m \le M} u_m t^{v_m} 1_{t < 1})t^{s-1}dt+ \sum_{m \le M}\frac{u_m }{s+v_m}$$
