Notice that in the picture you linked to, the poles moved from $\pm ia$ toward the real axis along counter-clockwise arcs that did not cross the real axis itself (where the integration contour lies). It is this path of approach of the poles toward the real axis that determines how the integration contour needs to be deformed. Informally, the poles make little "dents" in the integration contour as they approach it. If the poles were moved from $\pm ia$ toward the real axis along clockwise arcs, then these dents would be in the other direction and you would get a different sign in front of the residual part of the integral and hence a different answer overall. This shows that the final answer for the integral has branch cuts in its analytic dependence on $a$.
The answer to your main question depends on the way the location of the pole at $x_z=\arcsin(f(z^2)^{-1/2})$ as a function of $z$. You need to explicitly pick a path in the complex $z$-plane that takes you from $a$ to $-ia$ such that the pole $x_z$ does not cross your integration contour. As the pole approaches the real axis, it will make a dent in your integration contour, which will determine the way you need to deform it.
Let's record what is possible, so far. The case $n=2$, from the comments: $$\int_0^1\frac{dy}{x(x+y)}=\frac1{x^2}\int_0^1\frac{dy}{1+\frac{y}x}=\int_0^1\sum_{k\geq0}\frac{(-1)^ky^k}{x^{k+2}}dy=\sum_{k\geq0}\frac{(-1)^k}{k+1}\frac1{x^{k+2}}.$$ Hence, $$\int_1^{\infty}\sum_{k\geq0}\frac{(-1)^k}{k+1}\frac{dx}{x^{k+2}}=\sum_{k\geq0}\frac{(-1)^k}{(k+1)^2}=\frac{\zeta(2)}2.$$
The case $n=3$: $$\begin{aligned}\int_0^1\frac{dz}{x(x+y)(x+y+z)}&=\frac1{x(x+y)^2}\int_0^1\frac{dz}{1+\frac{z}{x+y}}\\&=\frac1x\sum_{k\geq0}\int_0^1\frac{(-1)^kz^kdz}{(x+y)^{k+2}}\\&=\frac1x\sum_k\frac{(-1)^k}{(k+1)(x+y)^{k+2}}.\end{aligned}$$ Hence,
$$\frac1x\int_0^1\sum_k\frac{(-1)^kdy}{(k+1)(x+y)^{k+2}}=\frac1x\sum_k\frac{(-1)^k}{(k+1)^2}\left[\frac1{x^{k+1}}-\frac1{(x+1)^{k+1}}\right].$$ Now, integrate with respect to $x$ (standard):
$$\int_1^{\infty}\frac{dx}{x^{k+2}}=\frac1{k+1} \qquad \text{and} \qquad \int_1^{\infty}\frac{dx}{x(x+1)^{k+1}}=\sum_{j=k+1}^{\infty}\frac1{j\cdot 2^j}.$$
Therefore, we compute the two series:
$$\sum_{k=1}^{\infty}\frac{(-1)^{k-1}}{k^3}=\frac34\zeta(3) \qquad \text{and} \qquad \sum_{k=1}^{\infty}\frac{(-1)^{k-1}}{k^2}\sum_{j=k}^{\infty}\frac1{j\cdot 2^j}=\frac{13}{24}\zeta(3).$$
We arrived at the valued predicted by Peter Mueller,
$$\int_1^{\infty}\int_0^1\int_0^1\frac{dz\,dy\,dx}{x(x+y)(x+y+z)}=\frac5{24}\zeta(3).$$
Caveat. One may anticipate higher values of $n$ to scale up the challenge.
UPDATE. Regarding Fedor's question, one contention is a follows: the sum in question is a weight 3 polylog, so it is a rational combination of $\zeta(3), \zeta(2)\log 2$ and $\log^3(2)$. Since a large numerical agreement verifies equality with only $\frac{13}{24}\zeta(3)$, it must be the exact evaluation.
UPDATE. I like to address the request from GH from MO directed to Agno.
Then, a direct answer to Fedor's question.
This time, we start integrating with respect to $x$:
$$\int_1^{\infty}\frac{dx}{x(x+y)(x+y+z)}
=\frac{y\,\log(1+y)+z\,\log(1+y)-y\,\log(1+y+z)}{yz(y+z)}.$$
Next, integrate in the variable $y$:
$$\int_0^1\frac{y\,\log(1+y)+z\,\log(1+y)-y\,\log(1+y+z)}{yz(y+z)}\,dy
=\frac{\text{Li}_2(z+2)-\text{Li}_2(z+1)-\text{Li}_2(2)}z;$$
where $\text{Li}_2(z)$ is the dilogarithm function
$$\text{Li}_2(z)=\int_1^z\frac{\log t}{1-t}\,dt.$$
Finally, we integrate in the last variable $z$:
$$\begin{align} \int_0^1\frac{\text{Li}_2(z+2)-\text{Li}_2(z+1)-\text{Li}_2(2)}z\,dz
&=\int_0^1\left(\int_1^{z+2}\frac{\log t}{1-t}\,dt-\int_1^{z+1}\frac{\log t}{1-t}\,dt-\int_1^2\frac{\log t}{1-t}\,dt
\right)\frac{dz}z \\
&=\int_0^1\left(\int_{z+1}^{z+2}\frac{\log t}{1-t}\,dt-\int_1^2\frac{\log t}{1-t}\,dt
\right)\frac{dz}z \\
&=\int_0^1\left(\int_2^{z+2}\frac{\log t}{1-t}\,dt-\int_1^{z+1}\frac{\log t}{1-t}\,dt\right)\frac{dz}z \\
&=\int_2^3\frac{\log t}{1-t}\left(\int_{t-2}^1\frac{dz}z\right)dt-\int_1^2\frac{\log t}{1-t}\left(\int_{t-1}^1\frac{dz}z\right)dt \\
&=\int_2^3\frac{\log t\,\log(t-2)}{t-1}\,dt-\int_1^2\frac{\log t\,\log(t-1)}{t-1}\,dt \\
&=\int_0^1\frac{\log (t+2)\,\log t}{t+1}\,dt-\int_0^1\frac{\log(t+1)\,\log t}t\,dt \\
&=\int_0^1\frac{\log (t+2)\,\log t}{t+1}\,dt+\frac34\zeta(3).
\end{align}$$
For the first integral in the last equality, write $\log(t+2)=\log 2+\log(1+\frac{t}2)$ and apply Taylor series:
$$\begin{align}
\int_0^1\frac{\log (t+2)\,\log t}{t+1}\,dt
&=\log 2\int_0^1\frac{\log t}{t+1}\,dt+\sum_{n\geq1}\frac{(-1)^{n-1}}{2^nn}\int_0^1\frac{t^n\log t}{t+1}\,dt \\
&=-\frac12\zeta(2)\,\log2+\sum_{n\geq1}\frac{(-1)^{n-1}}{2^nn}\int_0^1\frac{t^n\log t}{t+1}\,dt \\
&=-\frac12\zeta(2)\,\log2+\sum_{n\geq1}\frac{(-1)^{n-1}}{2^nn}\left[\frac{(-1)^{n-1}}2\zeta(2)+(-1)^{n-1}\sum_{k=1}^n\frac{(-1)^k}{k^2}\right] \\
&=-\frac12\zeta(2)\,\log2+\frac12\zeta(2)\sum_{n\geq1}\frac1{2^nn}+\sum_{n\geq1}\frac1{2^nn}\sum_{k=1}^n\frac{(-1)^k}{k^2} \\
&=-\frac12\,\zeta(2)\,\log2+\frac12\,\zeta(2)\,\log2+\sum_{n\geq1}\frac1{2^nn}\sum_{k=1}^n\frac{(-1)^k}{k^2} \\
&=\sum_{k\geq1}\frac{(-1)^k}{k^2}\sum_{n=k}^{\infty}\frac1{2^nn}.
\end{align}$$
The above derivations indicate we do not need Agno's $\log(e^x\pm1)$ integrals, instead we got Robert Z's $\log$-integral which sends us to his useful link evaluating as $\frac{13}{24}\zeta(3)$.
Best Answer
Here are a few contours of diverse complexity (click on the image for the link to the source publication).