  In 2003, Zhao discovered a curious congruence involving harmonic series and
Bernoulli numbers: for any odd prime $p$, $$\sum_{\substack{i,j,k\ge
1\\\gcd(ijk,p)=1\\i+j+k=p}}\frac{1}{ijk}\equiv -2B_{p-3} \pmod{p},$$ where
$B_n$ is the $n$-th Bernoulli number. This congruence was generalized by Wang
and Cai in 2014, and Cai, Shen and Jia in 2017 by replacing the odd prime $p$
in the summation and modulus with an odd prime power, and a product of two odd
prime powers, respectively. In particular, Cai, Shen and Jia proposed a
conjectural congruence: for any positive integer $n$ with an odd prime factor
$p$ such that $p^r \parallel n$ where $r\ge 1$, $$\sum_{\substack{i,j,k\ge
1\\\gcd(ijk,n)=1\\i+j+k=n}}\frac{1}{ijk}\equiv -2B_{p-3}\cdot \frac{n}{p}\cdot
\prod_{\substack{\text{prime $q\mid n$}\\q\ne
p}}\left(1-\frac{2}{q}\right)\left(1-\frac{1}{q^3}\right) \pmod{p^r}.$$ In this
paper, we establish the following generalization of their conjecture: for any
positive integer $n$ with an odd prime factor $p$ such that $p^r \parallel n$
where $r\ge 1$, $$\begin{aligned} \sum_{\substack{i,j,k\ge
1\\\gcd(ijk,n)=1\\a_1 i+a_2 j+a_3 k=An}}\frac{1}{ijk}&\equiv -2B_{p-3}\cdot
\frac{n}{p}\cdot \frac{Ag^3}{3}\left(\frac{1}{a_1^2 g_1^2}+\frac{1}{a_2^2
g_2^2}+\frac{1}{a_3^2 g_3^2}\right)\\ &\quad\times \prod_{\substack{\text{prime
$q\mid n$}\\q\ne p}}\left(1-\frac{2}{q}\right)\left(1-\frac{1}{q^3}\right)
\pmod{p^r}, \end{aligned}$$ where $a_1$, $a_2$ and $a_3$ are positive integers
coprime to $p$, and $A$ is a positive common multiple of $a_1$, $a_2$ and
$a_3$. Also, $g_1=\gcd(a_2,a_3)$, $g_2=\gcd(a_3,a_1)$, $g_3=\gcd(a_1,a_2)$ and
$g=\gcd(a_1,a_2,a_3)$.
