I suppose that all matrices in question are $p \times p$. It is clear that no closed-form result for the integral can be obtained if $p=1$, so it is doubtful that such an expression exists for general $p$, especially since there is no such expression for the generalized hypergeometric function ${}_0F_1$ of matrix argument.
Next, I note that in multivariate statistical analysis, where this integral was derived, the notation $\color{blue}{\mathbf{\Theta\Sigma^{-1}}}$ is shorthand for a symmetric matrix $\mathbf{M}$ that has the same eigenvalues as the matrix $\color{blue}{\mathbf{\Theta\Sigma^{-1}}}$.
I provide an approach to the integral that is derived from the paper, "Integral transform methods in goodness-of-fit testing, II: The Wishart distributions," Ann. Inst. Statist. Math. 72 (2020), 1317-1370; see p. 1328, Proposition 2 in that paper.
Let us rewrite the integral in the form
$$
f(\mathbf{M}_1,\mathbf{M}_2) := \int_{\mathbf{V}>\mathbf{O}} \mathbf{V} \, |\mathbf{V}|^\alpha \exp(-{\rm{tr}} \, \mathbf{M}_1^{-1}\mathbf{V}) \, {}_0F_1(\beta; \mathbf{M}_2 \mathbf{V}) \frac{\mbox{d}\mathbf{V}}{|\mathbf{V}|^{(p+1)/2}},
$$
where $\alpha$, $\beta$, $\mathbf{M}_1$, and $\mathbf{M}_2$ are trivially expressible in terms of your original notation, e.g., $\mathbf{M}_1 = 2 \mathbf{\Sigma}$, etc. With this notation, $\mathbf{M}_1$ is a positive definite (and symmetric) matrix, and $\mathbf{M}_2$ is a symmetric matrix.
Make the substitution $\mathbf{V} \to \mathbf{M}_1^{1/2} \mathbf{V} \mathbf{M}_1^{1/2}$. It is well-known that the measure $\mbox{d}\mathbf{V}/|\mathbf{V}|^{(p+1)/2}$ is invariant under this substitution; so, after simplification, we obtain
\begin{align*}
f(\mathbf{M}_1,\mathbf{M}_2) = |\mathbf{M}_1|^\alpha \mathbf{M}_1^{1/2} \int_{\mathbf{V}>\mathbf{O}} \mathbf{V} |\mathbf{V}|^\alpha \exp(-{\rm{tr}} \, \mathbf{V}) \, {}_0F_1\left(\beta; \mathbf{M}_1^{1/2} \mathbf{M}_2 \mathbf{M}_1^{1/2} \mathbf{V} \right) \frac{\mbox{d}\mathbf{V}}{|\mathbf{V}|^{(p+1)/2}} \mathbf{M}_1^{1/2}.
\end{align*}
That is,
$$
f(\mathbf{M}_1,\mathbf{M}_2) = |\mathbf{M}_1|^\alpha \mathbf{M}_1^{1/2} f(\mathbf{I}_p,\mathbf{M}) \mathbf{M}_1^{1/2},
$$
where $\mathbf{I}_p$ denotes the $p \times p$ identity matrix, and $\mathbf{M} = \mathbf{M}_1^{1/2} \mathbf{M}_2 \mathbf{M}_1^{1/2} $.
Now consider the integral,
$$
g(\mathbf{M}) := f(\mathbf{I}_p,\mathbf{M}) = \int_{\mathbf{V}>\mathbf{O}} \mathbf{V} \, |\mathbf{V}|^\alpha \exp(-{\rm{tr}} \, \mathbf{V}) \, {}_0F_1(\beta; \mathbf{M} \mathbf{V}) \frac{\mbox{d}\mathbf{V}}{|\mathbf{V}|^{(p+1)/2}},
$$
where $\mathbf{M}$ is any $p \times p$ symmetric matrix. Observe that $g(\mathbf{M})$ is a symmetric $p \times p$ matrix, each of whose entries is an integral. That the integral defining $g(\mathbf{M})$ converges absolutely for all $\alpha > (p-1)/2$ and all $\mathbf{M}$ can be proved using the Poisson integral for the Bessel functions of matrix argument; see Herz (Ann. Math., 61 (1955), 474--523).
Denote by $O(p)$ the group of all $p \times p$ orthogonal matrices. Note that if $\mathbf{H} \in O(p)$ then
$$
g(\mathbf{H} \mathbf{M} \mathbf{H}^{-1}) = \int_{\mathbf{V}>\mathbf{O}} \mathbf{V} \, |\mathbf{V}|^\alpha \, \exp(-{\rm{tr}} \, \mathbf{V}) \, {}_0F_1\left(\beta; \mathbf{H} \mathbf{M} \mathbf{H}^{-1} \mathbf{V} \right) \frac{\mbox{d}\mathbf{V}}{|\mathbf{V}|^{(p+1)/2}},
$$
Making the substitution $\mathbf{V} \to \mathbf{H} \mathbf{V} \mathbf{H}^{-1}$, using the invariance of the measure $\mbox{d}\mathbf{V}/|\mathbf{V}|^{(p+1)/2}$, and simplifying, we obtain
$$
g(\mathbf{H} \mathbf{M} \mathbf{H}^{-1}) = \mathbf{H} \int_{\mathbf{V}>\mathbf{O}} \mathbf{V} |\mathbf{V}|^\alpha \exp(-{\rm{tr}} \, \mathbf{V}) \, {}_0F_1(\beta; \mathbf{M} \mathbf{V}) \frac{\mbox{d}\mathbf{V}}{|\mathbf{V}|^{(p+1)/2}} \mathbf{H}^{-1}.
$$
In short, we have shown that, for all $p \times p$ symmetric matrices $\mathbf{M}$ and all $\mathbf{H} \in O(p)$,
$$
g(\mathbf{H} \,\mathbf{M} \,\mathbf{H}^{-1}) = \mathbf{H} \,g(\mathbf{M}) \,\mathbf{H}^{-1}.
$$
This property has appeared earlier in the 1970's and 1980's papers of K. I. Gross and R. A. Kunze on generalized Bessel functions of matrix argument. Gross and Kunze called it an orthogonal "covariance" property to distinguish it from the property of invariance under $O(p)$.
Since every symmetric matrix can be diagonalized by a transformation of the form $\mathbf{M} \to \mathbf{H} \mathbf{M} \mathbf{H}^{-1}$, the covariance property also shows that, in calculating $g(\mathbf{M})$, it suffices to assume that $\mathbf{M}$ is diagonal.
The above covariance property is the best that one can do for general $p$. For $p=2$, I suggest that you try to evaluate $g(\mathbf{M})$ by using an explicit formula for the ${}_0F_1$ function when $p=2$ (see Muirhead's book, Aspects of Multivariate Statistical Theory, Wiley, New York, 1982) and then calculating each entry of $g(\mathbf{M})$ term-by-term.
For the case in which $\mathbf{M} = \mathbf{I}_p$ (or a multiple of the identity), there is some hope that the matrix $g(\mathbf{M})$ can be calculated explicitly for general $p$. For $\mathbf{M} = \mathbf{I}_p$, the covariance property reduces to
$$
g(\mathbf{I}_p) = \mathbf{H} \, g(\mathbf{I}_p) \, \mathbf{H}^{-1}
$$
for all $\mathbf{H} \in O(p)$. By Schur's lemma, it follows that $g(\mathbf{I}_p) = c \, \mathbf{I}_p$ for some constant $c$. Therefore, all off-diagonal entries of $g(\mathbf{I}_p)$ are equal to zero. As for the diagonal entries, by taking traces, we obtain
\begin{align*}
c p = {\rm{tr}} \,(c \,\mathbf{I}_p)
&= {\rm{tr}} \, g(\mathbf{I}_p) \\
&= {\rm{tr}} \, \int_{\mathbf{V}>\mathbf{O}} \mathbf{V} \, |\mathbf{V}|^\alpha \, \exp(-{\rm{tr}} \, \mathbf{V}) \, {}_0F_1(\beta; \mathbf{V}) \frac{\mbox{d}\mathbf{V}}{|\mathbf{V}|^{(p+1)/2}} \\
&= \int_{\mathbf{V}>\mathbf{O}} ({\rm{tr}} \, \mathbf{V}) \, |\mathbf{V}|^\alpha \, \exp(-{\rm{tr}} \, \mathbf{V}) \, {}_0F_1(\beta; \mathbf{V}) \frac{\mbox{d}\mathbf{V}}{|\mathbf{V}|^{(p+1)/2}} \\
&= - \frac{\partial}{\partial t} \int_{\mathbf{V}>\mathbf{O}}|\mathbf{V}|^\alpha \, \exp(- t \, {\rm{tr}} \, \mathbf{V}) \, {}_0F_1(\beta; \mathbf{V}) \frac{\mbox{d}\mathbf{V}}{|\mathbf{V}|^{(p+1)/2}}\Bigg|_{t=1}.
\end{align*}
Expanding the function ${}_0F_1(\beta; \mathbf{V})$ in a series of zonal polynomials $Z_\kappa$, integrating term-by-term using a formula for the Laplace transform of $Z_\kappa$, and then differentiating term-by-term, one obtains a final result in terms of an infinite series of zonal polynomials.
Best Answer
Proceeding along Barabara's idea:
Let $f(t) = \log (A + tT) $. We have \begin{align} &\left\| \frac{f(h)-f(0)}{h} - \int_0^\infty (xI + A)^{-1}T(xI + A)^{-1} \mathrm{d} x\right\|\\ =\ & \left\| \int_0^\infty \frac{(xI + A)^{-1} - (xI+A+hT)^{-1} }{h} - (xI + A)^{-1}T(xI + A)^{-1} \mathrm{d} x\right\|. \end{align} Denote $B = xI + A$. We have \begin{align} B^{-1} - (B+hT)^{-1} &= B^{-1}[I - B(B+hT)^{-1}] \\ &= B^{-1}[I - (I + hTB^{-1})^{-1}]\\ &= B^{-1}[I - (I + hTB^{-1})^{-1}(I + hTB^{-1} - hTB^{-1})]\\ &= h B^{-1}(I + hTB^{-1})^{-1} TB^{-1}. \end{align} So, we have \begin{align} &\frac{B^{-1} - (B+hT)^{-1}}{h} - B^{-1}TB^{-1}\\ =\ & B^{-1}(I + hTB^{-1})^{-1} TB^{-1} - B^{-1}TB^{-1}\\ =\ & B^{-1}[(I + hTB^{-1})^{-1} - I]TB^{-1}\\ =\ & B^{-1}[(I + hTB^{-1})^{-1}(I + hTB^{-1} - hTB^{-1}) - I]TB^{-1}\\ =\ & - h( B + hT)^{-1}T B^{-1}TB^{-1}. \end{align} So, we have \begin{align} &\left\| \frac{f(h)-f(0)}{h} - \int_0^\infty (xI + A)^{-1}T(xI + A)^{-1} \mathrm{d} x\right\|\\ =\ & |h|\left\|\int_0^\infty (xI + A + hT)^{-1}T (xI + A)^{-1}T(xI + A)^{-1} \mathrm{d} x\right\| \to 0. \end{align} Remark: One need to prove that $\int_0^\infty (xI + A)^{-1}T (xI + A)^{-1}T(xI + A)^{-1} \mathrm{d} x$ is finite. Omitted here.