maybe we start with untangling some confusion.
The probability generating function is your function $G(z)$ and it is actually defined via expectation: $\mathbb E[z^X]$. However, what you are looking for is the moment generating function, which is similarly defined as follows:
$$M_X(t):=\mathbb E[exp(Xt)]=\int_{\mathbb R}exp(xt)f_X(t)dt$$
We can indeed use the tower property that you mentioned.
We first note: $$M_Z(t)=(1-p+p\cdot exp(t))^N$$
if $p$ where NOT random. How do we get rid of the randomness? That is doable via tower property:
$$M_Z(t)=\mathbb E[exp(Zt)]=\mathbb E[\mathbb E[exp(Zt)|Y]]$$
To compute that look at the quantity: $\mathbb E[exp(Zt)|Y=y]$. This has the form as mentioned above, but with p defined as given. Now this is a function $g(y)$ of y, a number for each fixed y. To compute the outer expectation with respect to the given exponential distribution, just integrate with respect to the exponential density:
$$M_Z(t)=\int_{\mathbb R^+}g(y)\lambda_0 \exp(-\lambda_0 y)dy.$$
If you need more explicit help, feel free to ask. Hope this helps without spoiling the fun ;)
EDIT: I just add my computation for the moment generating function for $N=1$. We keep in mind, how the exponential distribution is defined and we will utilize that to keep us from computing any integrals explicitly.
$$M_Z(t)=\mathbb E_Y[\left(1+(1-e^{-\lambda_1 Y})(e^t-1)\right)]$$
$$=\int_0^\infty\left(1+(1-e^{-\lambda_1 y})(e^t-1)\right)\lambda_0e^{-\lambda_0 y}dy$$
$$=\int_0^\infty\lambda_0e^{-\lambda_0y}dy+(e^t-1)\int_0^\infty(1-e^{-\lambda_1y})\lambda_0e^{-\lambda_0 y}dy$$
$$=1+(e^t-1)\left(1-\frac{\lambda_0}{\lambda_0+\lambda_1}\cdot\int_0^\infty(\lambda_0+\lambda_1)e^{-(\lambda_0+\lambda_1) y}dy\right)$$
$$=1+(e^t-1)\left(1-\frac{\lambda_0}{\lambda_0+\lambda_1}\right)
=1+(e^t-1)\cdot\frac{\lambda_1}{\lambda_0+\lambda_1}$$
$$=\frac{\lambda_0}{\lambda_0+\lambda_1}+e^t\frac{\lambda_1}{\lambda_0+\lambda_1}=q+e^t(1-q)$$
However using the binomial theorem for other values of $N$, will get us to $$M_Z(t)=(q+(1-q)e^t)^N$$
which just shows, that the $Z$ is binomial again, i.e. $Z\sim Bin(N,1-q)$. In other words, randomizing the success parameter of a binomial distribution with an exponential parameter does not change the type of distribution.
This explains also, why the inverse Fourier-transform yields something constant. The random variable $Z$ is just discrete and thus has no (absolutely continuous) density (with respect to Lebesgue-measure).
Consider the complement event. For instance, we can use the union bound as follows:
$$
\mathrm{P}(|X-\mu|\geq a) = \mathrm{P}(\{X-\mu\geq a\} \cup \{X-\mu\leq -a\} )\leq \mathrm{P}(X-\mu\geq a)+ \mathrm{P}(X-\mu\leq -a).
$$
This is a generic approach and you may use other tricks for particular problems.
Best Answer
There are several common notions of tightness of bounds, below is perhaps the simplest one.
Denote the Chernoff bound as $B(x) \equiv \frac{ \lambda }{ \lambda - r} e^{- rx}$ for the exponential function, which tail probability (complement CDF) is $P(X > x) = 1 - F_X(x) = e^{-\lambda x}$.
The error $B(x) - e^{-\lambda x}$ is readily shown to be always positive over $x \in \mathbb{R}^{+}$ for all $0 < r < \lambda$.
Consider the total error over $\mathbb{R}^{+}$, which is an routine integral: \begin{align} \Delta(r) \equiv \int_{ x = 0}^{\infty} \left( \frac{ \lambda }{ \lambda - r} e^{- r x} - e^{-\lambda x} \right)\,\mathrm{d}x = \frac{ \lambda}{ r\, (\lambda - r)} - \frac1{\lambda} \end{align} Note that the relative error $(B(x) - e^{-\lambda x})/e^{-\lambda x}$ doesn't converge as is. You may consider various transformation to make it work if the notion of relative error is really important to you.
Moving on. One can already guess that the minimum of $\Delta(r)$ occurs at $r_0 = \lambda / 2$ in the middle because the function is mirror symmetric (exchange $r \longleftrightarrow \lambda-r$). Anyway, checking that with the first and second derivative of $\Delta(r)$ is routine.
As $r$ varies, the curve $B(x)$ may get closer or further to $P(X > x)$ point-wise at any given $x$, depending on the relation between $r$ and $\lambda$. As a whole, one may consider the $r = r_0 = \lambda / 2$, which minimizes $\Delta(r)$ as defined above, as producing the "tightest" bound.
Sometimes, depending on the applications, one might care only about a specific interval of $a < x < b$. One then optimizes $r$ by minimizing the corresponding error over that interval only. It will be a different $r_0$ for different situation.
BTW, here the error over $\mathbb{R}^{+}$ is minimized to $\min \{\Delta(r)\} = \Delta(\frac{\lambda}2) = 3/\lambda$.
If your so-called Markov bound refers to the common $P( X > x) < \frac{ E[X] }x$, then it is just $B_M(x) = \frac1{\lambda x}$.
Note that the Markov inequality is very loose and the corresponding error doesn't converge, as in it's basically $\int \frac1{x} dx$ that blows up.
The Markov bound as such is a simple function, and one can certainly compare the functional forms of the Chernoff bound $B(x)$ with $B_M(x)$, either in general, or for given $\lambda = 1$ and discuss the numerical details at specific points of $x$.
I'm afraid not much (or too much) can be said without further context.