Preliminaries and Notation
I can proof your result without using the explicit transition density. The approach could generalize to the case, where the drift $\mu$ is a potential, so set $\mu(x)= -\nabla H(x)$. Then in the case of the OU-Process $H(x)=\frac{1}{2}x^2$. So we consider a process of the form $$ \mathrm{d} X_t = - \nabla H(X_t) \mathrm{d}t + \sqrt{2} \mathrm{d} W_t . $$ Then this process under some growth assumptions on $H$ has a unique invariant measure absolute continuous wrt. the Lebesgue measure $\mathrm{d} x$ given by $\mu(\mathrm{d} x) = \frac{1}{Z} \exp(- H(x))\mathrm{d}x $.
The evolution of the density $p_t$ is given by the Fokker-Planck-Equation $$ \partial_t p_t = \nabla\cdot ( \nabla p_t + p_t \nabla H) . $$ If we start the process not with a Dirac measure, but with a measure $p_0$ absolute continuous wrt. to $\mu$. We find a $\rho_t$ with $p_t = \rho_t \mu$. This relative density satisfies the equation $$ \partial_t \rho_t = \Delta \rho_t - \nabla H \cdot \nabla \rho_t =: L \rho_t . $$ The operator $-L$ is symmetric in $L^2(\mu)$ and under some growth condition on $H$ also elliptic.
Formal solution
Assume that the first eigenvalues of $-L$ are non-degenerated, i.e. $\lambda_0=0 < \lambda_1 < \lambda_2 \leq \lambda_3 \dots$. Then for a initial density $p_0$ absolute continuous wrt. to $\mu$ we have the expansion for $\rho_t$ $$ \rho_t = 1 + \langle \rho_0, v_1 \rangle_{L^2(\mu)} v_1 e^{-\lambda_1 t} + O(e^{-\lambda_2 t}) , $$ where $v_1$ is the eigenfunction of $L$ to the eigenvalue $\lambda_1$. Now we find $$ \begin{split} \mathbb{E}_{\rho_t \mu} f- \mathbb{E}_{\mu} f & = \int f(x) (\rho_t - 1)\mathrm{d}\mu \\ &= \int f(x) \langle \rho_0,v_1\rangle_{L^2(\mu)} e^{-\lambda_1 t} v_1 \mathrm{d}\mu + O(e^{-\lambda_2 t}) \end{split}$$ Multiplying by $e^{\lambda_1 t}$ we find $$ e^{\lambda_1 t} (\mathbb{E}_{\rho_t \mu} f- \mathbb{E}_{\mu} f) = \langle \rho_0 , v_1 \rangle_{L^2(\mu)} \int f(x) v_1 \mathrm{d} \mu + O(e^{-(\lambda_2-\lambda_1)t}). \tag{1}$$
OU-semigroup
Now setting $H(x) = \frac{1}{2} x^2$ we have in 1d $$ L = \partial_{xx} - x \partial_x .$$ The invariant measure is Gaussian $\mu = \frac{1}{\sqrt{2 \pi}} e^{-\frac{1}{2} x^2}$ and the eigenfunctions are the probabilistic Hermite polynomials with eigenvalues $0,-1,-2,\dots$. Hence $\lambda_1=1$ and $v_1(x)=x$ and $(1)$ takes the form $$ e^{t} (\mathbb{E}_{\rho_t \mu} f- \mathbb{E}_{\mu}f) = \langle \rho_0, x\rangle_{L^2(\mu)} \int f(x) x \mathrm{d}\mu + O(e^{-t})$$ Integrating by parts and using the fact $\mu'(x) = -x \mu(x)$ yields $$ e^{t} (\mathbb{E}_{\rho_t \mu} f- \mathbb{E}_{\mu}f) = \langle \rho_0, x\rangle_{L^2(\mu)} \int f'(x)\mathrm{d}{\mu} + O(e^{-t}).$$
Now with an approximation argument we take a series of inital distribution $\rho_0^n$ such that $\rho_0^n \mu \to \delta_{\xi}$ and find that $\langle \rho_0^n , x\rangle_{L^2(\mu)}\to \xi$, which shows your formula.
Problems in the general case
We used the spectral gap of the operator $L$, which holds in a wide class of potential $H$, but is hard to find explicit. There are estimates to bound it from below. Also sharp asymptotics for small noise are known.
Further we need the first eigenfunction, which corresponds to the metastable state of the system. Again for general $H$ explicit formulas are not known, but there exists again small noise expansions.
For your formula in the case of the OU-semigroup we heavily used the fact, that $-v_1(x) \mu(x) = \mu'(x)$, which is also not true in general. So the general formula would be $(1)$, which maybe can again be simplified, depending on the form of $v_1$.