Yes it's possible.
Suppose that $X\sim \mathrm{Gamma}(\alpha, \beta)$ where $\alpha$ is the shape parameter and $\beta$ the scale parameter.
Define the incomplete gamma function as
$$
\Gamma(a, z)=\int_{z}^{\infty}t^{a-1}\mathrm{e}^{-t}\;\mathrm{d}t
$$
and the generalized incomplete gamma function as
$$
\Gamma(a, z_0, z_1)=\int_{z_0}^{z_1}t^{a-1}\mathrm{e}^{-t}\;\mathrm{d}t = \Gamma(a,z_0)-\Gamma(a,z_1).
$$
Further, define the generalized regularized/normalized incomplete gamma function as
$$
Q(a,z_0,z_1)=\frac{\Gamma(a, z_0, z_1)}{\Gamma(a)}
$$
Finally, the inverse gamma regularized/normalized function is the solution for $z_1$ in $s = Q(a, z_0, z_1)$
The $q$-quantile of a Gamma distribution is given by $\beta\;\mathrm{InverseGammaRegularized(\alpha, 0, q)}$. Assuming that the $q$-quantile is $k$, solving for $\beta$ we have
$$
\beta = \frac{k}{\mathrm{InverseGammaRegularized(\alpha, 0, q)}}
$$
For your example
$$
\beta = \frac{130}{\mathrm{InverseGammaRegularized(2, 0, 0.9)}}=33.4214
$$
I'm not sure if this function is implemented in Julia.