First observe that on a compact Riemann manifold $(M, g)$ the operator $1+t \Delta$, $t>0$. $\Delta = d^*d: C^\infty(M)\to C^\infty(M)$ has a unique fundamental solution. Jacques Hadamard has constructed very explicit asymptotic expansions for this fundamental solution which lead to convergent series in the case of real analytic manifolds and metrics.
A modern description of Hadamard's construction can be found in volume 3, Sec. 17.4 of L. H"ormander's four volume on linear partial differential operators. (Hadamard's original memoir is also very useful, but harder to penetrate.) I will give a brief description of the fundamental solution $S_r$ of $(r+\Delta)$, $r>0$.
For $\nu=0,1,2,\dotsc $ and $r>0$ denote by $F_{\nu,r}(x)$ the generalized function (a.k.a. distribution) on $\mathbb{R}^n$ described as $\newcommand{\ii}{\boldsymbol{i}}$ a Fourier transform of a temperate distribution.
$$F_{\nu,r}(x)= \nu! (2\pi)^{-n} \int_{\mathbb{R}^n} e^{ \ii\langle x,\xi\rangle} (|\xi|^2+r)^{-\nu-1} d\xi. $$
The function $F_\nu$ can be expressed explicitly in terms of Bessel functions. Note that
$$F_{\nu,r}(x)= r^{\frac{n-\nu-1}{2}} F_\nu(\sqrt{r} x),\;\;F_\nu(x):=F_{\nu,r=1}(x). $$
If $$\Delta=-\sum_k\partial^2_{x_k}$$
denotes the (geometers') Laplacian in $\mathbb{R}^n$$ then
$$(r+\Delta)F_{0,r}=\delta_0,\;\;(r+\Delta)F_{\nu,r}=\nu F_{\nu-1, r},\;\;\forall \nu>0. $$
One can show that the generalized function $F_\nu$ depends only on the distance $|x|$.
Going back to the Riemann manifold $(M,g)$ we denote by $d: M\times M\to \mathbb{R}$ the geodesic distance function.
The Green function $G(x,y)$ then has an asymptotic expansion
$$ G(x,y)\sim \sum_{\nu=0}^\infty U_\nu(x,y) F_{\nu,r}( \; d(x,y)\;)$$
valid for $d(x,y)$ sufficiently small, where the functions $U_\nu(x,y)$ are explicitly described in the above reference. If $(M,g)$ is real analytic, then the above series converges in an appropriate sense.
This asymptotic expansion ought to be enough to investigate your question.
Update I want to add a "philosophical" comment. The question you asked is a special case of the following more general question.
Suppose that $f:\mathbb{R}\to\mathbb{R}$ is a continuous function. For simplicity, let us assume it is also bounded. We can define the bounded symmetric operator $f(\sqrt{\Delta})$ where $\Delta$ is the Laplacian on an $m$-dimensional manifold $M$. Investigate the behavior of $f(\varepsilon\Delta)$ as $\varepsilon \to 0$.Your case corresponds to $f(x)=(1+x^2)^{-1}$. The heat equation problems correspond to $f(x)=e^{-x^2}$. Suppose that $f$ is a symbol of order $k$, where $k$ could be $-\infty$. For example $(1+x^2)^{-k}$ is a symbol of order $-2k$, while $e^{-x^2}$ is a symbol of order $-\infty$.
In any case, when $f$ is a symbol, then $f(\Delta)$ is a pseudodifferential operator, and as such it has a Schwartz kernel which is a distribution on $M\times M$. $\newcommand{\ve}{\varepsilon}$ Your question is about the behavior as $\ve\to 0$ of the Schwartz kernel of $f(\ve \Delta)$ along normal directions to the diagonal of $M\times M$.
If $f$ is rapidly decaying at $\infty$, say $f(x) < (1+x^2)^{-m}$, $m=\dim M$, then the Schwartz kernel of $f(\Delta)$ is given by a continuous function and one can be quite precise about the behavior of the kernel of $f(\ve \Delta)$. In fact, the faster the decay of $f$ at $\infty$, the more accurate one can be about the behavior of the Schwartz kernel of $f(\ve \Delta)$. The radial symmetry you are talking about is then a simple consequence if $f$ decays faster than $|x|^{-N}$, $N$ sufficiently large. (I believe that $N>2m$ ought to do it but I don't want to be too firm.) If $f$ has exponential decay at $\infty$ one can be remarkable accurate and recover the radial symmetry you are mentioning. Your question involves the symbol $(1+x^2)^{-1}$ that isn't decaying fast enough at $\infty$. Translation: your problem requires a bit of care.