Supercongruence for the sporadic sequences B and F

Let $p\ge 5$ be a prime, and let $\\{u_n\\}_{n\ge 0}$ be one of the sequences $\mathbf{B}$ and $\mathbf{F}$, defined by

$$u_n=\sum_{k=0}^n(-1)^k3^{n-3k}\binom{n}{3k}\binom{3k}{2k}\binom{2k}{k}$$

or

$$u_n=\sum_{k=0}^n(-1)^k8^{n-k}\binom{n}{k}\sum_{j=0}^k\binom{k}{j}^3.$$

Here $\left(\frac{p}{3}\right)$ is the Legendre symbol, $B_{p-2}(x)$ is the Bernoulli polynomial, and $\\{\mathcal{U}_n\\}_{n\ge 1}$ is an integer sequence. The supercongruence for $\mathbf{B}$ and $\mathbf{F}$. For all positive integers $n$, the sequence satisfies

$$u_{np}\equiv u_n+\frac{1}{2}p^2\left(\frac{p}{3}\right)B_{p-2}\left(\frac{1}{3}\right)\mathcal{U}_n.$$

The sequence $\\{\mathcal{U}_n\\}_{n\ge 1}$ is independent of $p$. This is a numerical conjecture about supercongruences for two sporadic Apéry-like sequences, with the correction term expressed using a Bernoulli polynomial; the displayed values of $\mathcal{U}_n$ in the source provide computational evidence but no proof or resolution is supplied.