Long Comment. Here are three more representations of the integral
$$ I_N(p) = I_N(p; a_1, \ldots, a_N) := \int_{0}^{\infty} e^{-px^2} \prod_{n=1}^{N} \operatorname{erf}(a_nx) \, \mathrm{d}x, $$
where $p$ and $a_n$'s are a priori assumed to satisfy $p > 0$ and $a_1, \ldots, a_N > 0$ and then extended to a larger domain by the principle of analytic continuation, if available. (For example, the integral converges absolutely and locally uniformly on the region $\operatorname{Re}(p) > 0$ and $\left| \arg(a_n) \right| < \frac{\pi}{4}$ for all $n$, hence $I_N$ defines an analytic function in that region.)
Representation 1. The starting point is the following simple formula:
$$ \operatorname{erf}(x) = \frac{2x}{\sqrt{\pi}} \int_{0}^{1} e^{-x^2 t^2} \, \mathrm{d}t. $$
Plugging this into the defining integral of $I_N$,
$$\begin{align*}
I_N(p)
&= \int_{0}^{\infty} e^{-px^2} \prod_{n=1}^{N} \left[ \frac{2a_n x}{\sqrt{\pi}} \int_{0}^{1} e^{-a_n^2 x^2 t_n^2} \, \mathrm{d}t_n \right] \, \mathrm{d}x \\
&= \left( \frac{2}{\sqrt{\pi}} \right)^N \int_{[0,1]^N} a_1 \cdots a_N \int_{0}^{\infty} x^N e^{-(p+\sum_{n=1}^{N} a_n^2 t_n^2) x^2} \, \mathrm{d}x \mathrm{d}t_1 \cdots \mathrm{d}t_N \\
&= \left( \frac{2}{\sqrt{\pi}} \right)^N \Gamma\left(\frac{N+1}{2}\right) \int_{[0,1]^N} \frac{a_1 \cdots a_N}{2( p+\sum_{n=1}^{N} a_n^2 t_n^2 )^{\frac{N+1}{2}}} \mathrm{d}t_1 \cdots \mathrm{d}t_N \tag{1}
\end{align*}$$
The formula $\text{(1)}$ allows extending $I_N$ to a larger domain where $p+\sum_{n=1}^{N} a_n^2 t_n^2 \neq 0$ for all $t_1, \ldots, t_N \in [0, 1]$. Also, when $N = 1$ or $N = 2$, this integral can be computed analytically, yielding the same formula as in OP. (Note that this techinque is actually a disguise of the Feynman's trick, so it is essentially equivalent to @Eli Bartlett's approach.)
Representation 2. For any $x \in \mathbb{C}$, we have
$$ \operatorname{erf}(x) = e^{-x^2} \sum_{k=0}^{\infty} \frac{x^{2k+1}}{(k+\frac{1}{2})!}. $$
Here, the factorial of the fractional argument is computed via the identification $s! = \Gamma(s+1)$. From this, it follows that
$$\begin{align*}
I_N(p)
&= \int_{0}^{\infty} e^{-px^2} \prod_{n=1}^{N} \left[ e^{-(a_n x)^2} \sum_{k_n \geq 0} \frac{(a_n x)^{2k_n+1}}{(k_n+\frac{1}{2})!} \right] \, \mathrm{d}x \\
&= \sum_{k_1, \ldots, k_N \geq 0} \left[ \prod_{n=1}^{N} \frac{a_n^{2k_n+1}}{(k_n+\frac{1}{2})!} \right] \int_{0}^{\infty} x^{2(k_1 + \cdots + k_N) + N} e^{-(p+a_1^2 + \cdots + a_N^2)x^2} \, \mathrm{d}x \\
&= \sum_{k_1, \ldots, k_N \geq 0} \left[ \prod_{n=1}^{N} \frac{a_n^{2k_n+1}}{(k_n+\frac{1}{2})!} \right] \frac{\Gamma(k_1 + \cdots + k_N + \frac{N+1}{2})}{2(p + a_1^2 + \cdots + a_N^2)^{k_1 + \cdots + k_N + \frac{N+1}{2}}}. \tag{2}
\end{align*}$$
However, this series is not suitable for analyzing the limiting behavior or analytic properties of $I_N$.
Representation 3. This time, we invoke the following representation,
$$ \operatorname{erf}(x) = 1 - \frac{1}{\pi} \int_{-\frac{\pi}{2}}^{\frac{\pi}{2}} e^{-x^2\sec^2(\theta)} \, \mathrm{d}\theta,$$
valid for $\left| \arg x \right| < \tfrac{\pi}{4}$. Applying this, we get
$$\begin{align*}
I_N(p)
&= \int_{0}^{\infty} e^{-px^2} \prod_{n=1}^{N} \biggl( 1 - \frac{1}{\pi} \int_{-\frac{\pi}{2}}^{\frac{\pi}{2}} e^{-(a_nx)^2\sec^2(\theta_n)} \, \mathrm{d}\theta_n \biggr)\, \mathrm{d}x \\
&= \sum_{I \subseteq [N]} \left(-\frac{1}{\pi}\right)^{|I|} \int_{0}^{\infty} e^{-px^2} \prod_{n \in I} \biggl(\int_{-\frac{\pi}{2}}^{\frac{\pi}{2}} e^{-(a_nx)^2\sec^2(\theta_n)} \, \mathrm{d}\theta_n \biggr)\, \mathrm{d}x \\
&= \frac{\sqrt{\pi}}{2} \sum_{I \subseteq [N]} \left(-\frac{1}{\pi}\right)^{|I|} J(p ; (a_n)_{n\in I} ), \tag{3}
\end{align*}$$
where
$$ J(p ; \varnothing) = \frac{1}{\sqrt{p}}
\qquad\text{and}\qquad J(p ; (a_n)_{n\in I} ) = \int_{(-\frac{\pi}{2}, \frac{\pi}{2})^I} \frac{\prod_{n \in I} \mathrm{d}\theta_n}{\sqrt{p + \sum_{n\in I} a_n^2\sec^2(\theta_n) }} $$
for $I \neq \varnothing$. When $I$ consists of a single element, $J$ admits a closed-form:
$$ J(p ; a) = \frac{2}{\pi\sqrt{p}} \arctan\left( \frac{\sqrt{p}}{a} \right). $$
Similar to $\text{(1)}$, the formula $\text{(3)}$ also allows extending $I_N$ to a larger domain. Also, although $\text{(3)}$ is more complicated than $\text{(1)}$, I believe $\text{(3)}$ is better suited for analyzing the branch cuts of $I_N$ (and hence the monodromy of $I_N$). For example, when $a_1, \ldots, a_N \neq 0$, only the term $J(p ; \varnothing)$ induces the branch cut from $0$ to $\infty$.
