On a congruence involving harmonic series and Bernoulli numbers

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)$.