  Let $m$, $r$ and $n$ be positive integers. We denote by ${\bf k}\vdash n$ any
tuple of odd positive integers ${\bf k}=(k_1,\dots,k_t)$ such that
$k_1+\dots+k_t=n$ and $k_j\ge 3$ for all $j$. In this paper we prove that for
every sufficiently large prime $p$ $$ \sum_{\substack{l_1+l_2+\cdots+l_n=mp^r
p\nmid l_1 l_2 \cdots l_n }} \frac1{l_1l_2\cdots l_n} \equiv p^{r-1} \sum_{{\bf
k}\vdash n} C_{m,{\bf k}} B_{p-{\bf k}} \pmod{p^r} $$ where $B_{p-{\bf
k}}=B_{p-k_1}B_{p-k_2}\cdots B_{p-k_t}$ are products of Bernoulli numbers and
the coefficients $C_{m,{\bf k}}$ are polynomials of $m$ independent of $p$ and
$r$. This generalizes previous results by many different authors and confirms a
conjecture by the authors and their collaborators.
