As a matter of fact, it isn't hard to construct a multiplicative sequence $a_n$ such that $f(z)$ is an entire function without zeroes. Unfortunately, it is completely useless for the questions that you brought up as "motivation".
Here is the construction.
Claim 1: Let $\lambda_j\in [0,1]$ ($j=0,\dots,M$). Assume that $|a_j|\le 1$ and $\sum_ja_j\lambda_j^p=0$ for all $0\le p\le P$. Then, if $P>2eT$, we have $\left|\sum_j a_je^{\lambda_j z}\right|\le (M+1)(eT/P)^P$
Proof: Taylor decomposition and a straightforward tail estimate.
Claim 2: Let $P$ be large enough. Let $\Delta>0$, $M>P^3$ and $(M+1)\Delta<1$. Let $I_j$ ($0\le j\le M$) be $M+1$ adjacent intervals of length about $\Delta$ each arranged in the increasing order such that $I_0$ contains $0=\lambda_0$. Suppose that we choose one $\lambda_j$ in each interval $I_j$ with $j\ge 1$. Then for every $|a_0|\le 1$, there exist $a_j\in\mathbb C$ such that $|a_j|\le 1$ and $\sum_{j\ge 0} a_j\lambda_j^p=0$ for $0\le p\le P$.
Proof: By duality, we can restate it as the claim that $\sum_{j\ge 1}|Q(\lambda_j)|\ge |Q(0)|$ for every polynomial $Q$ of degree $P$. Now, let $I$ be the union of $I_j$. It is an interval of length about $M\Delta$, so, by Markov's inequality, $|Q'|\le CP^2(M\Delta)^{-1} A\le CP^{-1}\Delta^{-1}A$ where $A=\max_I |Q|\ge |Q(0)|$. But then on the 5 intervals $I_j$ closest to the point where the maximum is attained, we have $|Q|\ge A-5\Delta CP^{-1}\Delta^{-1}A\ge A/2$. The rest is obvious.
Claim 3: Suppose that $a_0$ is fixed and $a_j$ ($j\ge 1$) satisfy $\sum_{j\ge 0} a_j\lambda_j^p=0$ for $0\le p\le P$. Then we can change $a_j$ with $j\ge 1$ so that all but $P+1$ of them are exactly $1$ in absolute value and the identities still hold.
Proof: As long as we have more than $P+1$ small $a_j$, we have a line of solutions of our set of $P+1$ linear equations. Moving along this line we can make one of small $a_j$ large. Repeating this as long as possible, we get the claim.
Now it is time to recall that the logarithm of the function $f(z)$ is given by
$$
L(z)=\sum_{n\in\Lambda}a_ne^{-z\log n}
$$
where $\Lambda$ is the set of primes and prime powers and $a_n=m^{-1}a_p^m$ if $n=p^m$. We are free to choose $a_p$ for prime $p$ in any way we want but the rest $a_n$ will be uniquely determined then. The key point is that we have much more primes than prime powers for unit length.
So, split big positive numbers into intervals from $u$ to $e^\Delta u$ where $\Delta$ is a slowly decaying function of $u$ (we'll specify it later). Formally we define the sequence $u_k$ by $u_0=$something large, $u_{k+1}=e^{\Delta(u_k)}u_k$ but to put all those backward apostrophes around formulae is too big headache, so I'll drop all indices. Choose also some slowly growing functions $M=M(u)$ and $P=P(u)$ to be specified later as well.
We need a few things:
1) Each interval should contain many primes. Since the classical prime number theorem has the error term $ue^{-c\sqrt{\log u}}$, this calls for $\Delta=\exp\{-\log^{\frac 13} u\}$
Then we still have at least $u^{4/5}$ primes in each interval (all we need is to beat $u^{1/2}$ with some margin).
2) We should have $M\Delta\ll 1$, $M>P^3$, and $u\left(\frac{eT}P\right)^P\le (2u)^{-T-3}$ for any fixed $T>0$ and all sufficiently large $u$. This can be easily achieved by choosing $P=\log^2 u$ and $M=\log^6 u$.
3) At last, we'll need $M(P+\sqrt u)\ll u^{4/5}$, which is true for our choice.
Now it is time to run the main inductive construction. Suppose that $a_n$ are already chosen for all $n$ in the intervals up to $(u,e^{\Delta}u)$ and we still have almost all primes in the intervals following the current interval free (we'll check this condition in the end). We want to assign $a_p$ for all $p$ in our interval for which the choice hasn't been made yet or was made badly. We start with looking at all $a_p$ that are not assigned yet or assigned in a lame way, i.e., less than one in absolute value. Claim 3 (actually a small modification of it) allows us to upgrade all of them but $P+1$ to good ones (having absolute value $1$) at the expense of adding an entire function that in the disk of radius $T$ is bounded by $(2u)^{T} u\left(\frac{eT}P\right)^P\le u^{-3}$ to $L(z)$. Now we are left with at most $\sqrt u$ powers of primes and $P+1$ lame primes to take care of. We need the prime powers participate in small sums as they are and we need the small coefficients to be complemented by something participating in small sums too. For each of them, we choose $M$ still free primes in the next $M$ intervals (one in each) and apply Claim 2 to make a (lame) assignment so that the corresponding sum is again bounded by $u^{-3}$ in the disk of radius $T$. We have at most $u$ such sums, so the total addition will be at most $u^{-2}$. This will finish the interval off. Now it remains to notice that we used only about $\sqrt u+P$ free primes in each next interval and went only $M$ intervals ahead. This means that in each interval only $M(\sqrt u+P)$ free primes will ever be used for compensating the previous intervals, so we'll never run out of free primes. Also, the sum of the blocks we constructed will converge to an entire function. At last, when $\Re z>1$, we can change the order of summation and exponentiate finally getting the Dirichlet series representation that we need.
The end.
A.F. Leont'ev continued to work on general Dirichlet series well into 1980s (until his death in 1987). Actually, he published three monographs on the subject from 1976 to 1983! He made a short summary of his earlier results for the 1974 ICM in Vancouver (a free preview of the lecture is available here).
A.F. Leont'ev obtained in some sense final results on the representation of analytic functions by general Dirichlet series of the form
$$f(s)=\sum\limits_{n=1}^{\infty}a_n e^{-\lambda_n s},\quad s\in D\subset \mathbb C.$$
He studied Dirichlet series in bounded and unbounded convex domains (including half-planes). The problem is that his results may not be directly applicable to the `ordinary' Dirichlet series with $\lambda_n=\ln n$. A typical Leont'ev's theorem for half-planes is as follows (link to the original article in Russian).
Theorem. For every $\rho>1$, there is a sequence $\lambda_n>0$, $n\in\mathbb N$, satisfying the condition
$$\lim\limits_{n\to\infty}\frac{n}{\lambda_n^\rho}=\tau,\quad 0<\tau<\infty,$$
such that any function $f$, which is analytic in the right half-plane $\Re z > 0$ , can be represented in the form
$$f(z)=\sum\limits_{n=1}^{\infty}a_n e^{-\lambda_n z}+\Phi(z),\qquad \Re z > 0,$$
where $\Phi$ is entire.
This obviously doesn't cover the case $\lambda_n=\ln n$, $n\in\mathbb N$.
Anyway, if you're interested and if you have a colleague who speaks Russian I can send you a couple of original articles by Leont'ev (PDF files). (Edit: sent.) By the way, the papers you've mentioned both deal with the case of convergence in a bounded domain $D$.
Best Answer
It should be possible to make a Dirichlet series whose domain of meromorphicity is as screwy as you want. Notice that $\zeta(s - 1 - \alpha) = \sum n^{1+\alpha}/n^s$, so $\zeta(s - 1-\alpha)$ is a Dirichlet series, with pole at $\alpha$. Let $\gamma$ be a curve dividing the complex plane into two pieces, one of which contains all $z$ with $\Re(z)$ sufficiently large; and choose $\alpha_i$ a sequence of points of $\gamma$ which is dense in $\gamma$. Consider $$\sum c_i \zeta(s -1- \alpha_i)$$ where $c_i$ is some sequence which goes to zero very fast.
As long as the $c_i$ go to zero fast enough, this should converge absolutely on away from $\gamma$; and will be represented by a Dirichlet series on the half plane to the right of the rightmost point of $\gamma$. The dense set of poles will ensure that you can't continue past $\gamma$. Details are left to the reader. :)