Let us define the Bernoulli numbers $B_{i}$, for all natural numbers $i \geq 0$, by
(1) $\frac {X}{e^{X} - 1} = \sum_{i=0}^{\infty} \frac{B_{i}}{i!}X^{i}$,
where $X$ is a variable. (The first member is a well defined formal power series in $X$, let us say over $\mathbb{C}$, thus this definition makes sense.)
A well known theorem of von Staudt has the following special case :
(2) For every prime number $p$, $p B_{p-1}$ is a $p$-integral rational number congruent with -1 modulo $p$.
I got the following strengthening as a by-product :
(3) For every prime number $p$, $p B_{p-1}$ is a $p$-integral rational number congruent with $p + (p-1) !$ modulo $p^{2}$.
(In view of Wilson's theorem, (3) implies (2).)
The proof that I can extract from my work is very heavy, so I wonder if someone could know or could find a more straightforward proof. Here is the mine. It is true for $p = 2$, thus we can assume that $p$ is odd.
The formal power series $Log (1-X)$ is divisible by $X$ (in the ring $\mathbb{C}[[X]]$ of formal power series in $X$), thus we can replace $X$ by $Log (1-X)$ in (1). We find
(4) $\frac{Log (1-X)}{-X} = \sum_{i=0}^{\infty} \frac{B_{i}}{i!}(Log (1-X))^{i}$.
Let $p$ be a fixed odd prime number. For every rational integer $i$, let us define $f_{i}(X)$ (the $i$-th Mirimanoff polynomial with regard to $p$) as
$f_{i}(X) = \sum_{r=1}^{p-1} r^{i-1} X^{r}$.
Then
(5) $Log (1-X) \equiv - f_{0}(X) - \frac{X^{p}}{p} (mod\ X^{p+1})$
in the ring $\mathbb{C}[[X]]$. Dividing both members and the modulus of this congruence by $X$ gives
(6) $\frac{Log (1-X)}{X} \equiv - \frac{f_{0}(X)}{X} - \frac{X^{p-1}}{p} (mod\ X^{p}).$
Thus (4) gives
(7) $\frac{f_{0}(X)}{X} + \frac{X^{p-1}}{p} \equiv \sum_{i=0}^{\infty} \frac{B_{i}}{i!}(Log (1-X))^{i} (mod\ X^{p})$
in the ring $\mathbb{C}[[X]]$.
Relation (5) implies the weaker relation
$Log (1-X) \equiv - f_{0}(X) (mod\ X^{p})$
thus (7) gives
(8) $\frac{f_{0}(X)}{X} + \frac{X^{p-1}}{p} \equiv \sum_{i=0}^{\infty} \frac{B_{i}}{i!} (-1)^{i} f_{0}(X)^{i} (mod\ X^{p})$
in the ring $\mathbb{C}[[X]]$. Since this is a congruence modulo $X^{p}$ and $f_{0}(X)^{i}$ is divisible by $X^{i}$, we can neglect the values $\geq p$ of $i$, so we obtain a congruence in the polynomial ring $\mathbb{Q}[X]$ :
$\frac{f_{0}(X)}{X} + \frac{X^{p-1}}{p} \equiv \sum_{i=0}^{p-1} \frac{B_{i}}{i!} (-1)^{i} f_{0}(X)^{i} (mod\ X^{p})$
in the polynomial ring $\mathbb{Q}[X]$. Isolating the values 0 and $p-1$ of $i$, we put this under the form
(9) $\frac{f_{0}(X)}{X} + \frac{X^{p-1}}{p} \equiv 1 + \sum_{i=1}^{p-2} \frac{B_{i}}{i!} (-1)^{i} f_{0}(X)^{i} + \frac{B_{p-1}}{(p-1)!} f_{0}(X)^{p-1} (mod\ X^{p})$
in the polynomial ring $\mathbb{Q}[X]$.
Now we need the following fact :
(10) for $1 \leq j \leq p-2$, the coefficient of $X^{p-1}$ in $f_{0}(X)^{j}$ is congruent to $- j !$ modulo $p$.
This can be deduced from the following result of Mirimanoff :
for $1 \leq j \leq p-2$,
$f_{0}(X)^{j} \equiv (-1)^{j-1} j! f_{p-j}(1-X) (mod\ p \mathbb{Z}_{(p)}[X] + X^{p} \mathbb{Z}_{(p)}[X])$,
where $\mathbb{Z}_{(p)}$ denotes the ring of $p$-integral rational numbers.
(See D. Mirimanoff, "L'équation indéterminée $x^{l} + y^{l} + z^{l} = 0$ et le critérium de Kummer", Journal für die reine und angewandte Mathematik, vol. 128, 1905, p. 45-68, spec. p. 60-61. It is on https://gdz.sub.uni-goettingen.de/ )
In view of (10) and the fact that $B_{1}, B_{2}, ..., B_{p-2}$ are $p$-integral, considering the terms in $X^{p-1}$ in (9) gives
$\frac{1}{p} \equiv - \sum_{i=1}^{p-2} (-1)^{i}B_{i} + \frac{B_{p-1}}{(p-1)!} (mod\ p)$
(defining that two non $p$-integral numbers are congruent modulo $p$ if their difference is $p$-integral and divisible by $p$ in $\mathbb{Z}_{(p)}$).
This result can also be written
(11) $\frac{1}{p} \equiv 1 - \sum_{i=0}^{p-2} (-1)^{i}B_{i} + \frac{B_{p-1}}{(p-1)!} (mod\ p)$.
But $\sum_{i=0}^{p-2} (-1)^{i}B_{i} \equiv 0 (mod\ p)$. (Note that $B_{1} = -1/2$ and $B_{i} = 0$ if $i$ is odd and > 1, and use Congruence between Bernoulli numbers )
Thus (11) gives
$\frac{B_{p-1}}{(p-1)!} \equiv \frac{1}{p} - 1 (mod\ p).$
Multiplying by $p (p-1)!$, we get a congruence modulo $p^{2}$ :
$p B_{p-1} \equiv (p-1)! - p (p-1)! (mod\ p^{2})$.
In view of Wilson's theorem, $ (p-1)! \equiv - 1 (mod\ p)$, thus $p (p-1)! \equiv -p (mod\ p^{2})$ and (3) is proven.
Let me indicate without proof that the following can be deduced from (3) :
If $p$ is a prime number $\geq 5$, if $j$ is a natural number $\geq 1$, then
$p B_{j(p-1)} \equiv j (1 + (p-1)!) + p - 1 (mod\ p^{2})$.
In view of Wilson's theorem, it is a strengthening of von Staudt's theorem (but only valid for $p \geq 5$).
As I said, I find that my proof of (3) is very heavy and I would be interested by a simpler proof, and also by a reference to the literature. Thanks in advance.
By Faulhaber's formula $$\sum_{k=1}^n k^{p-1}=\frac{n^p}{p} +\frac{n^{p-1}}{2}+\cdots+B_{n-1}n$$ where the omitted terms have coefficients that are $p$-integers (von Staudt). This means that $$B_{p-1}p\equiv\sum_{k=1}^p k^{p-1}\equiv\sum_{k=1}^{p-1} k^{p-1}\pmod{p^2}$$ So all we need is $$p+(p-1)!\equiv\sum_{k=1}^{p-1}k^{p-1}\pmod{p^2}.\tag{*}$$ This congruence was proved by Lerch. See this paper by Sondow and MacMillan.
Here's a summary of their proof. We can write $k^{p-1}=1+pa_k$ for $p\nmid k$. Then $$\sum_{k=1}^{p-1}k^{p-1}=-1+p+p(a_1+a_2+\cdots+a_{p-1})$$ and $$(p-1)!^{p-1}=(1+pa_1)(1+pa_2)\cdots (1+pa_{p-1}) \equiv1+p(a_1+a_2+\cdots+a_{p-1})\pmod {p^2}.$$ Write $(p-1)!=pb-1$ with $b\in\Bbb Z$. Then $$(p-1)!^{p-1}=(1-pb)^{p-1} \equiv1+bp\equiv(p-1)!+2\pmod{p^2}.$$ We conclude that $$\sum_{k=1}^{p-1}k^{p-1}\equiv-2+p+(p-1)!^{p-1}\equiv p+(p-1)! \pmod{p^2}.$$