You draw a picture in the y-z plane, and you rotate it by 45 degrees and blow up by a factor of $\sqrt{2}$. This gives you a triangle region, where the vertices are determined to be the image of the old vertices under the rotation: (x,x),(a,x),(a,a). The images in the $\beta,\gamma$ plane are $(2x,0),(a+x,a-x),(2a,0)$.
But I understand that you want a formal algorithm, so that you can do it mechanically. First, you make a one-variable transformation, where you trade in z for $\gamma= z-y$. If you really insist on doing the rotation (it's not the right thing), you can then trade in $y$ for $\beta$, using $\beta= 2y + \gamma$.
The answer for the first shift is that the z range is from y to a, so translating z by -y goes from 0 to a-y.
$$\int_0^a dx f(x) \int_x^a dy \int_0^{a-y} d\gamma g(\gamma) $$
To complete the mechanical change of variable, introduce an indicator function of one dimension $\phi(x)$, which is 1 for $x>0$ and 0 for $x<0$. Then you write the integral as follows:
$$ \int_0^a dx f(x) \int_x^a dy \int_0^\infty \phi(a-y-\gamma) g(\gamma) d\gamma$$
Rearrange the order of the y and $\gamma$ integrals and perform the y integral explicitly (using the fact that the indefinite integral of $\phi(x)$ is $x\phi(x)$) to get
$$\int_0^a dx f(x) \int_0^\infty (a-x-\gamma)\phi(a-x-\gamma)) g(\gamma) d\gamma $$
The $\phi$ function now gives you the new domain for $\gamma$
$$ \int_0^a dx f(x) \int_0^{a-x} (a-x-\gamma)g(\gamma)$$
and this is the best form. Using explicit indicator functions ($\phi$) this way is useful in many cases for getting explicit answers in closed form. If you really want to do the transformation of variables the way you said it, you introduce the indicator functions for the region, and then transform the region and find the new domain. This is just the same as drawing the triangle and rotating it.
Please see the region shaded in the diagram.
So with change of order, the integral will be
$\displaystyle \int_1^2 \int_y^{y^2} e^{\frac{x}{y}} \ dx \ dy$
This will require you to integrate $\displaystyle y \ e^y$ in the second integral, which can be done by Integration by parts.
Edit: If you combine your second and third integral, and then with the first, you get the same integral I wrote above. You just considered change of order for each sub-region separately and so ended up with three integrals.
Best Answer
First, change the dummy variables in each so that they play the same roles. In this case, I'd change $y = a$ in $I_1$ and $x = \tilde a$ in $I_2$:
$$I_1 = \int_0^\infty b(y)e^{\epsilon a} \int_y^\infty \epsilon e^{(\gamma-\epsilon) z} \int_z^\infty \tilde b(x) \rho e^{-D(x)-qx- \gamma x}\,dxdzdy\\ I_2 = \int_0^\infty \tilde{b}(x) \rho e^{-D(x)-qx - \gamma x} \int_0^x \epsilon e^{(\gamma - \epsilon)z} \int_{0}^{z} b(y) e^{\epsilon y}\,dydzdx$$
If we bring all the integrands inside, we see that they are the same function in both $I_1$ and $I_2$.
In $I_1$, we have $0 \le y$ for the first integral, $y \le z$ from the second integral, and $z \le x$ from the third integral. So the region of integration is $\{(x,y,z)\mid 0 \le y \le z \le x\}$.
In $I_2$, we have $0 \le x, y, z$ from the three integrals, $z \le x$ from the second integral, and $y \le z$ from the third integral. So once again, the region of integration is $\{(x,y,z)\mid 0 \le y \le z \le x\}$.
So both integrals are of the same function over the same region. If the volume integral over the region converges, then these two iterated integrals must also converge to the same value (as would any of the other four iterated integrals over that region).
In general, the volume integral need not converge even if both iterated integrals converge (and if the volume integral does not converge, there is no guarantee the iterated integrals will converge to the same value). However in this case the integrand is the product of separate functions for each variable. And if all three of these 1D functions converge when integrated from $0$ to $\infty$, the volume integral can be shown to converge as well.