Let's consider the second integral, which can be written in the following form:
$$
I(p, q, i, j, k, l; a, b, c, d)
:=
\int_0^\infty
dt\,
\exp(-p t^2)
t^q
j_i(a t) j_j(b t) j_k(c t) j_l(d t)
$$
where in your case, $i = j = k = 1$, $l = 0$, and $q = -1$.
These kinds of integrals (as well their generalization to a product of arbitrarily many spherical Bessel functions) are discussed in Fabrikant - Elementary exact evaluation of infinite integrals of the product of several spherical Bessel functions, power and exponential, where the main idea is to use the following identity:
\begin{align}
I(p, q, i, j, k, l; a, b, c, d)
&=
(-1)^{i+j+k+l}
a^i
b^j
c^k
d^l
\frac{\partial^i}{(a \partial a)^i}
\frac{\partial^j}{(b \partial b)^j}
\frac{\partial^k}{(c \partial c)^k}
\frac{\partial^l}{(d \partial d)^l}
\biggl[\\
&\int_0^\infty
dt
\exp(-p t^2)
\frac{
j_0(a t) j_0(b t) j_0(c t) j_0(d t)
}
{
t^{i + j + k + l - q}
}
\biggr].
\end{align}
The key point is to now expand the zeroth order spherical Bessel functions into trigonometric functions, and converting the products of the trigonometric functions into sums:
\begin{align}
\sin(ax)
\sin(bx)
\sin(cx)
\sin(dx)
=&
\frac{1}{8}
\biggl\{
\cos[(a + b + c + d)x]
+
\cos[(a + b - c - d)x]
+
\cos[(a - b + c - d)x]\\
&+
\cos[(a - b - c + d)x]
-
\cos[(-a + b + c + d)x]
-
\cos[(a - b + c + d)x]\\
&-
\cos[(a + b - c + d)x]
-
\cos[(a + b + c - d)x]
\biggr\}
\end{align}
followed by the use of the following integral, which is not considered in the reference above, but can be found in Gradshteyn and Ryzhik, 7th ed., formula 3.953.8:
$$
\mathcal{I}(p, s; n)
:=
\int_0^\infty dt\,
t^n
\exp(-p t^2)
\cos(s t)
=
\frac{1}{2}
p^{\frac{-(n + 1)}{2}} \,
e^{-s^2 / 4 p}
\Gamma \left(\frac{1}{2} + \frac{n}{2}\right) \,
_1F_1\left(-\frac{n}{2}; \frac{1}{2}; \frac{s^2}{4 p}\right).
$$
Note that the formal requirement is that $\operatorname{Re}(n) > -1$, but the above result can be understood as an analytic continuation for general values $p, s, n$.
Additionally, it can happen that one of the "angles" above is zero, in which case we have the integral:
$$
\mathcal{I}(p, 0; n)
:=
\int_0^\infty dt\,
t^n
\exp(-p t^2)
=
\frac{1}{2} p^{-\frac{n}{2}-\frac{1}{2}} \Gamma \left(\frac{n+1}{2}\right)
$$
with the same condition on $n$ as above.
The result in your specific case is then:
\begin{align}
I(p, -1, 1, 1, 1, 0; a, b, c, d)
&=
-
a
b
c
\frac{\partial}{(a \partial a)}
\frac{\partial}{(b \partial b)}
\frac{\partial}{(c \partial c)}
\int_0^\infty
dt
\exp(-p t^2)
\frac{
j_0(a t) j_0(b t) j_0(c t) j_0(d t)
}
{
t^4
}\\
&=
-
a
b
c
\frac{\partial}{(a \partial a)}
\frac{\partial}{(b \partial b)}
\frac{\partial}{(c \partial c)}
\int_0^\infty
dt
\exp(-p t^2)
\frac{
\sin(a t) \sin(b t) \sin(c t) \sin(d t)
}
{
a\, b\, c\, d\, t^8
}\\
&=
-
\frac{\partial}{\partial a}
\frac{\partial}{\partial b}
\frac{\partial}{\partial c}
\bigg[
\frac{1}{a\, b\, c\, d}
\int_0^\infty
dt
\exp(-p t^2)
\frac{
1
}
{
t^8
}
\frac{1}{8}
\bigg\{
\\
&\cos[(a + b + c + d)t]
+
\cos[(a + b - c - d)t]
+
\cos[(a - b + c - d)t]
\\
&+
\cos[(a - b - c + d)t]
-
\cos[(-a + b + c + d)t]
-
\cos[(a - b + c + d)t]
\\
&-
\cos[(a + b - c + d)t]
-
\cos[(a + b + c - d)t]
\bigg\}
\bigg]
\\
&=
-
\frac{1}{8}
\frac{\partial}{\partial a}
\frac{\partial}{\partial b}
\frac{\partial}{\partial c}
\biggl\{\frac{1}{a\, b\, c\, d}
\\
&\mathcal{I}(p, a + b + c + d; -8)
+
\mathcal{I}(p, a + b - c - d; -8)
+
\mathcal{I}(p, a - b + c - d; -8)
\\
&+
\mathcal{I}(p, a - b - c + d; -8)
-
\mathcal{I}(p, -a + b + c + d; -8)
-
\mathcal{I}(p, a - b + c + d; -8)
\\
&-
\mathcal{I}(p, a + b - c + d; -8)
-
\mathcal{I}(p, a + b + c - d; -8)
\biggr\}.
\end{align}
The explicit result is fairly cumbersome to fully write out; below is an example Mathematica code which can be used as a starting point to generate the full solution (when $a \pm b \pm c \pm d \neq 0$) and compare it with the numerical result:
numeric[p_, a_, b_, c_, d_] := NIntegrate[
Exp[-p t^2] SphericalBesselJ[1, a t] SphericalBesselJ[1,
b t] SphericalBesselJ[1, c t] SphericalBesselJ[0, d t]/t,
{t, 0, Infinity}
];
integral[p_, s_, n_] :=
1/2 p^(-(n + 1)/2) Exp[-s^2/(4 p)] Gamma[
1/2 + n/2] Hypergeometric1F1[-n/2, 1/2, s^2/(4 p)];
result = -1/8 Table[
Series[
expression,
{epsilon, 0, 0}
] // Normal // D[#/(a b c d), a, b, c] &,
{
expression,
{
integral[p, a + b + c + d, -8 + epsilon],
integral[p, a + b - c - d, -8 + epsilon],
integral[p, a - b + c - d, -8 + epsilon],
integral[p, a - b - c + d, -8 + epsilon],
-integral[p, -a + b + c + d, -8 + epsilon],
-integral[p, a - b + c + d, -8 + epsilon],
-integral[p, a + b - c + d, -8 + epsilon],
-integral[p, a + b + c - d, -8 + epsilon]
}
}
] // Total;
No idea if the solution which the code above generates can be simplified though.
