That's in fact a very good question. I've never thought about this before. I'm not sure whether the following is mathematically rigorous.
I'd start with the simple case a=1. The standard way to evaluate the integral is by setting
I=\int_{-\infty}^{\infty} \exp(-x^2)>0.
Then
I^2=\int_{\mathbb{R}^2} \mathrm{d} x \mathrm{d} y exp(-x^2-y^2).
Then we use polar coordinates in the xy plane to get
I^2=2 \pi \int_0^{\infty} \mathrm{d} r r \exp(-r^2)=-\pi \exp(-r^2)|_{r=0}^{\infty}=\pi.
Since I>0 we uniquely get
I=\sqrt{\pi}
with the usual positive square root for a positve real number.
For real a>0 then you get by substitution y=\sqrt{a} x
I_a=\int_{-\infty}^{\infty} \mathrm{d} x \exp(-a x^2) = \frac{1}{\sqrt{a}} \int_{-\infty}^{\infty} \exp(-y^2)=\frac{\pi}{\sqrt{a}},
where again, I used the usual positive square root of a positive real number. You could as well have substituted y=-\sqrt{a} x, but then the infinite boundaries change signs, and thus the integral also flips signs, so that you get the same positive result as it must be.
Now, for a \in \mathbb{C}, for convergence you obviously should have \mathrm{Re} \; a > 0. Again we can take both roots in the substitution above, but it's most convenient to write
a=|a| \exp(\mathrm{i} \varphi),
where \varphi \in (-\pi/2,\pi/2) (this is one of many possible choices of the argument for a complex number with positive real part) and then use the square root as
\sqrt{a}=\sqrt{|a|} \exp(\mathrm{i} \varphi/2), \quad \sqrt{|a|}>0. \qquad (*)
Now consider the Integral I_{a} as an integral in the complex x plane along the real axis. Then again we substitute x=\sqrt{a} t with the meaning of \sqrt{a} given by (*).
The real integration path in the complex x plane then maps to an integration path in the complex t plane, which is a straight line through the origina, running from the lower left quadrant into the upper right (for \varphi>0) or from the upper left to the lower right quadrant (for \varphi<0). In both cases, you can define a closed path by adding the real axis (run from right to left in both cases) and two vertical parts at infinity. If you integrate \exp(-t^2) along that closed path you get 0 due to Cauchy's integral theorem, and this means that instead to integrate along the straight line in the t plane you can as well integrate along the real axis (in the normal positive sense). Thus we find
I_a=\frac{\sqrt{\pi}}{|a|} \exp(-\mathrm{i} \varphi/2),
where the angle \varphi \in (-\pi/2,\pi/2).
