The congruences that you stated are not correct; correct versions are given in Ofir Gorodetsky's answer. The related congruences in Meštrović's paper are (37)–(41); see also his Remark 21 on page 13.

A related paper is doi.org/10.37236/1592.

Here's a generating function version of this argument. Let $p_n$ be the number of monic polynomials of degree $n$, and let $s_n$ be the number of square-free monic polynomials of degree $n$. Then the generating function for monic polynomials that are squares is $\sum_n p_n x^{2n}$, so the unique factorization of an arbitrary monic polynomial as a square-free polynomial times a square gives $\sum_{n}p_n x^n = \left(\sum_{n}s_n x^n\right)\left(\sum_n p_n x^{2n}\right)$. Since $p_n = q^n$, this gives $\sum_n s_n x^n = (1-qx^2)/(1-qx)=1+qx + \sum_{n=2}^\infty (q^n - q^{n-1})x^n$.

One thing you might try is multiply (3) by $z^s/s!$ and sum on $s$, and apply Taylor's theorem. This will give you a simple closed form for the generating function in $s$ of (1).