The estimate you seek is reminiscent of H. Weyl's tube formula. I will give you some pointers referring for more details to section 9.3.5. of these lectures.
Denote by $r$ the distance to $\newcommand{\pa}{\partial}$ $\pa \Omega$ $\newcommand{\bn}{\boldsymbol{n}}$ and by $\bn$ the innner pointing unit normal of $\pa \Omega$. $\newcommand{\bp}{\boldsymbol{p}}$ There exists $r_0>0$ such that the map $\newcommand{\bR}{\mathbb{R}}$
$$
\pa \Omega \times [0,r_0)\ni(\bp,r)\stackrel{\Phi}{\longmapsto} \bp+r\bn\in \bR^n
$$
is a diffeomorphism onto the region $\DeclareMathOperator{\dist}{dist}$
$$A_{\rho_0}:=\big\{\; x\in \bar{\Omega};\;\;\dist(x,\pa \Omega)<\rho_0\;\big\}.$$
Then
$$\int_{A_{\rho_0}} e^{-r^2/t} dx=\int_{\pa\Omega\times [0,r_0)}\Phi^*\big(e^{-\rho^2/t} dx\Big)=\int_{\pa\Omega\times [0,r_0)}e^{-\rho^2/t}\Phi^*(dx).
$$
The pullback of the Euclidean volume form $dx$ via the map $\Phi$ is described explicitly in the above reference. It has the form $\newcommand{\eQ}{\mathscr{Q}}$ $\DeclareMathOperator{\tr}{tr}$
$$
\Phi^*(dx)= \eQ_\bp(r)dV_{\pa \Omega}(\bp) dr,
$$
where, for each $\bp\in \pa \Omega$ $(r)$, $\eQ_\bp$ is a polynomial of degree $n-1$ in $r$
$$
\eQ_\bp(r)=\sum_{j=0}^{n-1}c_j(\bp) r^j=\det\big(1-r S_\bp\big),
$$
where $S_{\bp}$ is the second fundamental form of the hypersurface $\pa \Omega$ at the point $\bp$ defined in terms of the inner normal $\bn$. More precisely if $(x^i)$ are local coordinates on $\pa \Omega$ near $\bp$, then
$$
S_\bp(\pa_{x^i},\pa_{x^j})=\big(\; \bn(\bp),\pa^2_{x^ix^j}\bn(\bp)\;\big),
$$
where $(-,-)$ denotes the canonical inner product on $\bR^n$.
Thus
$$\int_{A_{\rho_0}} e^{-r^2/t} dx=\int_0^{r_0} e^{-r^2/t}\left(\int_{\pa\Omega} \eQ_{\bp}(r)dV_{\pa \Omega}(\bp)\right) dr. $$
The integral
$$
K(r):=\int_{\pa\Omega} \eQ_{\bp}(r)dV_{\pa \Omega}(\bp)
$$
appears in Weyl's tube formula and, more precisely $\newcommand{\bom}{\boldsymbol{\omega}}$ (see Eq. (9.3.18) in the above reference)
$$
K(r)=\sum_{k=0}^{n-1} (-1)^{n-1-k}\bom_{n-k}r^{n-k}\mu_k(\Omega),
$$
$$
=\bom_0\mu_{n-1}(\Omega)r-\bom_1\mu_{m-2}(\Omega)r^2+\cdots +(-1)^{n-1}\bom_n\mu_0(\Omega)r^n,
$$
where $\bom_m$ denotes the volume of the $m$-dimensional Euclidean unit ball and $\mu_k(\Omega)$ is the curvature measure of degree $k$. (You need to be careful about various sign conventions. In the above reference the second fundamental form is defined using the outer normal.)$\DeclareMathOperator{\vol}{vol}$ For example
$$
\mu_{n-1}(\Omega)= \frac{1}{2} \vol_{n-1}(\pa\Omega),
$$
$$
\mu_{n-2}(\Omega)=\frac{1}{2\pi}\int_{\pa \Omega} \tr S_\bp dV_{\pa\Omega}(\bp),
$$
where $\tr S_\bp$ is the mean curvature of $\pa\Omega$ at $\bp$. Also, $\mu_0(\Omega)$ is the Euler characteristic of $\Omega$.
The asymptotics of
$$
J(t)=\int_0^{r_0} e^{-r^2/t}K(r) dr,
$$
can be determined easily by making the change in variables $s=r^2/t$, $r=\sqrt{st}$ so that
$$
J(t)=\frac{\sqrt{t}}{2}\int_0^{r_0^2/t} e^{-s}K(\sqrt{st}) s^{-1/2} ds
$$
$$
=\frac{1}{2}\sum_{k=0}^{n-1}(-1)^{n-1-k}t^{\frac{n-k}{2}}\mu_k(\Omega)\int_0^{r_0^2/t} e^{-s} s^{\frac{n-k}{2}-1} ds.
$$
Observe that as $t\searrow 0$
$$
\int_0^{r_0^2/t} e^{-s} s^{\frac{n-k}{2}-1} ds=\;\underbrace{\int_0^{\infty} e^{-s} s^{\frac{n-k}{2}-1} ds}_{\Gamma\big( \frac{n-k}{2}\big)}\;+ O\big(t^{N}\big),\;\;\forall N>0.
$$
Hence
$$
J(t)=\frac{1}{2}\sum_{k=0}^{n-1}(-1)^{n-1-k}t^{\frac{n-k}{2}}\Gamma\Big(\;\frac{n-k}{2}\;\Big)\mu_k(\Omega) + O\big(t^{N}\big),\;\;\forall N>0.
$$
Finally
$$
|I(t)-J(t)|=\int_{\dist(x,\pa \Omega)>r_0} e^{-\dist(x,\pa\Omega)^2/t} dx
$$
$$
\leq e^{-r_0^2/t}\vol(\Omega)=O(t^{N}),\;\;\forall N>0.
$$
Hence
$$ I(t)=\frac{1}{2}\sum_{k=0}^{n-1}(-1)^{n-1-k}t^{\frac{n-k}{2}}\Gamma\Big(\;\frac{n-k}{2}\;\Big)\mu_k(\Omega) + O\big(t^{N}\big),\;\;\forall N>0.$$
The leading term of this yields the estimate indicated by Carlo Beenakker.
About the tube formula The tube formula for smooth domains or convex domains or more generally sets with positive reach are both special cases of the very general kinematic formulas. The Brunn-Minkowski formula is also a special case of the kinematic formula.
