Surely there are many: these are all polynomials in one variable, so every two of them are algebraically dependent because of the transcendence degree argument :-)
However, I am sure that this is not what you wanted to hear, so here you are a nice argument showing how to guess your formula and obtain other formulas somewhat similar to it. Note that there is a remarkable symmetry property $P_k(-1-N)=(-1)^{k+1} P_k(N)$ for $k>0$. (Basically, for $k>0$ the polynomial $P_k(x)$ is the only polynomial of degree $k+1$ solving the functional equation $f(x)-f(x-1)=x^k$ together with the condition $f(0)=0$, and then you can show that $Q_k(x)=(-1)^{k+1} P_k(-1-x)$ satisfied exactly the same conditions, which proves the symmetry property without any annoying computations.) If we re-define $P_0(N)=N+\frac12$ (and assume $P_{-1}=1$), this symmetry will hold in general.
Now, the polynomial $P_1^2$, as a polynomial of degree $4$, should be a rational combination of $P_0$, $P_1$, $P_2$ and $P_3$ (and such a combination is clearly unique - you yourself observed that they form a basis), and because of the type of symmetry it possesses, it is actually a combination of $P_1$ and $P_3$ (because other polynomials change sign under the symmetry $N\mapsto -1-N$, and this would contradict the linear independence), and looking at it carefully we observe that the $P_1$-coefficient is equal to zero, and the $P_3$-coefficient is equal to~$1$, which is your formula.
For the same reason, the product $P_mP_n$ is expressed as a linear combination of $P_l$ where $l\le m+n+1$, $l\equiv m+n+1\pmod{2}$, - half of the terms disappear for free! (And, because of vanishing at~$0$, the redefined $P_0=N+\frac12$ and $P_{-1}=1$ will not show up in such a combination if $m+n>0$.)
Some examples: $6P_1P_2=5P_4+P_2$, $3P_2^2=2P_5+P_3$, $12P_2P_3=7P_6+5P_4$, $2P_3^2=P_7+P_5$, $60P_3P_4=27P_8+35P_6-2P_4$ (this last one is a bit disappointing!) etc.
