Is $p_t(y,x)$ the kernel of the time reversal of the diffusion $X$, for $p_t(x,y)$ the kernel of $X$? Short version. If $X$ is a diffusion with generator $L$ and the Lebesgue measure is invariant for $X$, then $L^*$ has no term of order zero and it corresponds to another diffusion $X^*$. Denoting by $p$ and $p^*$ their kernels,

do we have $p^*_t(x,y)=p_t(y,x)$?

I am interested in full answers, but I would assume there is a reference for this.
In details. Let $X$ be a diffusion, solution to
$$\mathrm dX_t = b(X_t)\mathrm dt + \sum_i\sigma_i(X_t)\mathrm dW^i_t.$$
I am in a somewhat regular case: $b$ and all the $\sigma_i$ are smooth, the diffusion is Feller (the semigroup acts on continuous functions vanishing at infinity) and hypoelliptic, i.e. the vector fields satisfy the parabolic Hörmander condition when written in Stratonovich form (so the kernel is well-defined and smooth). I am also interested in cases with stronger hypotheses.
This diffusion has generator $L:f\mapsto \sum_kb^k\partial_kf+\frac12\sum_i\sum_{k\ell}\sigma_i^k\sigma_i^\ell\partial_{k\ell}f$, with (formal) adjoint
$$ L^*:f\mapsto
- \sum_k\partial_k\big(b^kf\big)+\frac12\sum_i\sum_{k\ell}\partial_{k\ell}\big(\sigma_i^k\sigma_i^\ell f\big)
= \frac12\sum_i\sum_{k\ell}\sigma_i^k\sigma_i^\ell\partial_{k\ell}f + (\text{lower order}).$$
If $L^*1=0$, i.e. there are no terms of order zero (the Lebesgue measure is preserved), this $L^*$ is the generator of another diffusion
$$\mathrm dX^*_t = b^*(X_t)\mathrm dt + \sum_i\sigma_i(X_t)\mathrm dW^i_t$$
with the same noise components. (I am abusing notation a bit, in the sense that $b^*$ is not the adjoint of $b$, since the former depends on the $\sigma_i$.)
Since $L$ is hypoelliptic, $X$ started at $x$ admits a smooth kernel $(t,y)\mapsto p_t(x,y)$, i.e.
$$\mathbb P_x(X_t\in A)=\int_Ap_t(x,y)\mathrm dy.$$
We can check (easier in Stratonovich form) that hypoellipticity of $L$ equivalent to that of $L^*$, so $X^*$ also has a smooth kernel $p^*$.

Question.
Is it true that for any $f,g$ smooth with compact support (say), we have
$$ \int f(x)p^*_t(x,y)g(y)\mathrm dx\mathrm dy = \int f(x)p_t(y,x)g(y)\mathrm dx\mathrm dy? $$

I call this a question about time reversals because it should be true that $t\in[0,T]\mapsto X_{T-t}$ has the same distribution as $t\in[0,T]\mapsto X^*_t$ when the distribution of the initial point is the Lebesgue measure (whatever that means), so $p_t(y,x)$ is the time reversal of the kernel $p_t(x,y)$ of $X$, and $p^*$ is the kernel of the time reversal $X^*$ of $X$.
Haussmann and Pardoux (open access) prove results about time reversals that can be used to answer my question in the positive if we have instead a finite equilibrium measure $\mu$ (i.e. $L^*\mu=0$ when we consider the adjoint with respect to $\mu$).
Formal proof. At the formal level, I see the following reasoning as evidence that the result should be true.
Let $f,g$ be nice functions, and consider the derivative of the integral
$$I:s\in(0,t)\mapsto \int p^*_s(x,y)p_{t-s}(z,y)g(z)\mathrm dy.$$
The Fokker-Planck equation asserts that $\partial_tp = L^*_yp$ and $\partial_tp^* = L^{**}_yp^* = L_yp^*$ ($y$ is the second variable), so
$$I'(s) = \int\Big((L_yp^*)_s(x,y)p_{t-s}(z,y) - p^*_s(x,y)(L^*_yp)_{t-s}(z,y)\Big)\mathrm dy = 0 $$
since $L^*$ is the adjoint of $L$ (this is more or less the same as saying that $s\mapsto p_{t-s}(z,X_s)$ is a local martingale, for any initial condition $X^*_0=x$ and point $z$, although this latter fact is actually true and easy to prove).
Hence we must have $I(0)=I(t)$ with
$$ I(0) = \int \delta_x(y)p_t(z,y)\mathrm dy
        = p_t(z,x) $$
and
$$ I(t) = \int p^*_t(x,y)\delta_z(y)\mathrm dy
        = p^*_t(x,z). $$
Maybe there are standard arguments that can prove the derivative and limits taken above are valid, but I have not been able to find them either myself or in the literature.
Edit: The dual process approach.
If I am not mistaken, one approach suggested by Mateusz Kwaśnicki in his answer is to go in the direction $X\mapsto\hat X\mapsto \hat L$ rather than $X\mapsto L^*\mapsto X^*$. It would be sufficient to carry through the following steps.

*

*Prove that there exists a dual process $\hat X$ (more or less the resolvent of the processes are in duality), in the sense of section 13.2 of [CW].

*Prove that the semigroup of $\hat X$ is given by $P_tf(x)=\int p_t(y,x)f(y)\mathrm dy$.

*Prove that $\hat X$ satisfies the dual SDE, so in fact $\hat X=X^*$.

These all come with some level of difficulty, but there are a few results in [CW] that may help. For instance, for 1. we can use Theorem 15.5, and all that remains to show is that the Lebesgue measure is invariant (in the strong sense $\int\mathbb P_x(X_t\in A)\mathrm dx=\mathrm{leb}(A)$) and that there exists a “cofine topology” in the sense of Definition 15.4. Point 2. is very close to the duality of the semigroups, as in Proposition 13.6, and point 3., once we have point 2., reduces to having enough functions in the domain of $\hat L$.
Although the route seems sound, altogether the results in [CW] look very measure-theoretic to me (unsurprisingly) and I am having trouble relating them to the smooth quantities I am considering.
[CW] Chung, Walsh, Markov Processes, Brownian Motion, and Time Symmetry.
 A: This is called duality of Markov processes, and it is well-known that the dual process is essentially the time-reversal. A standard reference is

*

*Chung, Walsh, Markov Processes, Brownian Motion, and Time Symmetry, DOI:10.1007/0-387-28696-9
Although it is a truly brilliant book, it is, unfortunately, not an easy read. Additionally, it avoids the use of generators, so it does not immediately answer your question. Thus, I'll try to briefly describe the picture.

Assume that we know that $X_t$ is a "nice" Markov process, the Lebesgue measure is invariant for $X_t$ (and proving this property is probably the hardest part of the story), and that $X_t$ has continuous transition densities $p_t(x, y)$.
Side note: more generally, it is enough to assume that there is an excessive measure (with full support, and perhaps some additional conditions) and that the $\lambda$-potential kernel $u_\lambda(x, A) = \int_0^\infty e^{-\lambda t} p_t(x, E) dt$ has a "nice" density function $u_\lambda(x, y)$ for every $\lambda > 0$.

Under the above assumptions, the dual transition kernel:
$$\hat p_t(x,y) = p_t(y,x)$$
corresponds to another "nice" Markov process, the dual process, denoted by $\hat X_t$. Furthermore, the time-reversed process $X_{T-t}$ (or, strictly speaking, its càdlàg version) is related to $\hat X_t$ in a way discussed in Section 13.4 of the Chung–Walsch; but in fact Chapters 10, 11 and 13 in this book deal mostly with time-reversal related concepts. Here let me only mention that, by definition,
$$ \iint f(x) p_t(x, y) g(y) dx dy = \iint g(y) \hat p_t(y, x) f(x) dx dy ,$$
and hence
$$ \int f(x) \mathbb E^x g(X_t) dx = \int g(y) \hat{\mathbb E}^y f(\hat X_t) dy .$$
Combined with the Markov property, this relatively easily leads to the following theorem: for every measurable functional $\Phi$ on the path space, we have
$$ \int \mathbb E^x \Phi((X_{T-t} : t \in [0, T])) dx = \int \hat{\mathbb E}^x \Phi((\hat X_t : t \in [0, T])) dx . $$
For a detailed and rigorous statement, see Theorem 13.18 in the book. (Note that, in the statement of this result, $r_t$ stands for time-reversal, defined just before Remark 13.16, and $\zeta$ is the life-time of the process.)

It remains to say why $\hat p_t(x, y)$ is generated by $L^*$. This is in fact quite simple, as long as we have sufficiently many functions in the domains of $L$ and $L^*$. If $f$ is in the domain of $L^*$, $g$ is in the domain of $L$, and both $f$ and $g$ are compactly supported, then
$$ \begin{aligned} \int f(x) L g(x) dx & = \int f(x) \biggl(\lim_{t \to 0^+} \frac{1}{t} \biggl(\int p_t(x, y) g(y) dy - g(x)\biggr)\biggr) dx \\ & = \lim_{t \to 0^+} \frac{1}{t} \int f(x) \biggl(\int p_t(x, y) g(y) dy - g(x) \biggr) dx \\ & = \lim_{t \to 0^+} \frac{1}{t} \biggl(\iint f(x) p_t(x, y) g(y) dx dy - \int f(x) g(x) dx \biggr) \\ & = \lim_{t \to 0^+} \frac{1}{t} \biggl(\iint g(y) \hat p_t(y, x) f(x) dx dy - \int f(y) g(y) dy \biggr) \\ & = \lim_{t \to 0^+} \frac{1}{t} \int g(y) \biggl(\int \hat p_t(y, x) f(x) dx - f(y) \biggr) dy \\ & = \int g(y) \biggl(\lim_{t \to 0^+} \frac{1}{t} \biggl(\int \hat p_t(y, x) f(x) dx - f(y)\biggr)\biggr) dy = \int g(y) \hat L f(y) dy ,\end{aligned} $$
and so $\hat L$ is the formal adjoint of $L$.
A: This is an immediate corollary of Theorem 3.50 in

D. W. Stroock, An introduction to the analysis of paths on a Riemannian manifold (2000).
Mathematical Surveys and Monographs. 74. Providence, RI: American Mathematical Society (AMS). xvii, 269 p.

See also Theorem 4.17 for a manifold version of it.
A: Edit: There is a flaw in the argument, one needs a bit more than stated. I am not sure how to fix the argument, but apparently there is a simple reference anyway. — Pierre PC
I will prove the result assuming a bit more regularity: $X$ and $X^*$ are both Feller (their semigroups map $\mathcal C_0$ to itself) and complete (defined for all times).
Associated to $X$ is a semigroup $t\mapsto P_t$, sending $f\in\mathcal C_0$ to $x\mapsto\mathbb E_x[f(X_t)]$, and a generator $L:D(L)\subset\mathcal C_0\to\mathcal C_0$ defined as
$$ Lf:=\mathcal \lim_{t\to0}\frac1t(P_tf-f) $$
whenever the limit exists (in the $\mathcal C_0$ topology). The same objects $P^*$ and $L^*$ are well-defined for $X^*$. The question can be rephrased as such: do we have
$$ \langle P_tf,g\rangle = \langle f,P^*_tg\rangle $$
for all $f,g\in\mathcal C^\infty_c$, where $\langle f,g\rangle:=\int f(x)g(x)\mathrm dx$?
I believe one can actually show that $X$ and $X^*$ are dual, in the sense of section 13.2 of [CW]. That said, it will not be necessary to get to this level of generality. I will nonetheless consider resolvent operators.
As is well-known (see for instance Chapter 1 in [D]), the norm $f\mapsto|f|_\infty+|Lf|_\infty$ turns $D(L)$ into a Banach space, and $L:D(L)\to\mathcal C_0$ is continuous in this topology. Denoting by $\mathrm{id}:D(L)\to\mathcal C_0$ the inclusion, we define whenever possible $R_z:\mathcal C_0\to D(L)$ as the (double sided) continuous inverse of $z\mathrm{id}-L$. This is at least well-defined for $\Re z>0$, in which case we have classically
$$R_zf = \int_0^\infty\mathsf e^{-zt}P_tf\mathrm dt.$$
The same considerations hold for $X^*$.

Proposition.
For every $z$ with positive real part, and any $f,g\in\mathcal C^\infty_c$, we have
$$ \langle R_zf,g\rangle = \langle f,R^*_zg\rangle. $$

Let us see first how this implies the result. For $f,g\in\mathcal C^\infty_c$,
$$x\mapsto |P_tf(x)g(x)|\leq|f|_\infty\cdot g(x)$$
is uniformly integrable in $x$, so we can use Fubini's theorem and deduce that the Laplace transforms
$$ \begin{align*}
     z\in(0,\infty)&\mapsto\int_0^\infty\mathsf e^{-zt}\langle P_tf,g\rangle, &
     z\in(0,\infty)&\mapsto\int_0^\infty\mathsf e^{-zt}\langle f,P^*_tg\rangle
   \end{align*} $$
are actually equal, which means that $\langle P_tf,g\rangle_\mu = \langle f,P^*_tg\rangle_\mu$ for almost all $t\geq0$, hence for all $t$ since they are easily seen to be continuous. Since these quantities restrict to continuous bilinear forms over $\mathcal C_0(U)$ for any relatively compact open set $U$, and $\mathcal C^\infty_c(U)\subset\mathcal C_0(U)$ is dense, they actually agree over all such $\mathcal C_0(U)$, hence over $\mathcal C^0_c$.
Proof of the proposition.
Because $X$ is solution to an SDE with smooth coefficients, $\mathcal C^\infty_c\subset D(L)$, so we can define $\mathcal E:=\{(z\mathrm{id}-L)f,f\in\mathcal C^\infty_c\}\subset\mathcal C^\infty_c$ and $\mathcal E^*$ similarly. By definition of $R_z$, $R_zf$ has compact support for $f\in\mathcal E$.
For $f\in\mathcal E$, $g\in\mathcal E^*$, and $z$ with positive real part, we have
$$ \langle R_zf,g\rangle
 = \langle R_zf,(z\mathrm{id}-L^*)R^*_zg\rangle
 = \langle (z\mathrm{id}-L)R_zf,R^*_zg\rangle
 = \langle f,R^*_zg\rangle. $$
Suppose for now that for every $f\in\mathcal C^\infty_c$, we can find a sequence of functions $f_n\in\mathcal E$ such that $(f_n)_{n\geq0}$ and $(R_zf_n)_{n\geq0}$ converge locally uniformly to $f$ and $R_zf$, and similarly for $\mathcal E^*$. Then, given $f,g$ smooth with compact support, we find such approximations $(f_n)_{n\geq0}$ and $(g_m)_{m\geq0}$, and
$$ 0
 = \langle R_zf_n,g_m\rangle - \langle f_n,R^*_zg_m\rangle
 \xrightarrow[n\to\infty]{} \langle R_zf,g_m\rangle - \langle f,R^*_zg_m\rangle
 \xrightarrow[m\to\infty]{} \langle R_zf,g\rangle - \langle f,R^*_zg\rangle. $$
Edit: This is (obviously) not true in general. The second limit has to be justified somehow, for instance showing that $R_zf$ is in $L^1$ and $g_m\to g$ holds in $\mathcal C_0$. We can choose such a sequence $g_m$ if $\mathcal C^\infty_c$ is a core for $L$, i.e. it is dense in $D(L)$. For the latter problem, I don't know of an easy criterion.
It remains to show that $\mathcal E$ is sufficiently dense in $\mathcal C^\infty_c$ in the sense described above. Fix $f\in\mathcal C^\infty_c$. Defining $\phi=R_zf\in D(L)\subset\mathcal C_0$, we know by definition that $(z\mathrm{id}-L)\phi = f$. Since $L$ is hypoelliptic, it means that $\phi$ is actually smooth. Using a sequence of cut-off functions, we find a sequence of functions $\phi_n\in\mathcal C^\infty_0$ such that for any given compact set $K$, $\phi_n$ and $\phi$ eventually agree over $K$. Defining $f_n:=(z\mathrm{id}-L)\phi_n\in\mathcal E$, it means in particular that $f$ and $f_n$ eventually agree over any fixed compact set. Of course $R_zf_n=\phi_n$ and $R_zf=\phi$, so the same holds for $R_zf_n$ and we are done.

[CW] Chung, Walsh, Markov Processes, Brownian Motion, and Time Symmetry.
[D] Dynkin, E. B., Markov processes.
