One way is to consider that Lorentz transformations apply more fundamentally to momentum/ energy, than to space/time. So, with a boost transformation along the $z$ axis, we will have : 

$\begin {pmatrix} p'_z \\E' \end{pmatrix} = \gamma(v)\begin {pmatrix} 1 & \frac{v}{c} \\\frac{v}{c}  &1\end{pmatrix}\begin {pmatrix} p_z \\E \end{pmatrix}$

It is not difficult to see that, if, $|\large \frac{p_zc}{E}|=1$, then $|\large  \frac{p'_zc}{E'}|=1$

Now, by dimensional analysis, we have : $\frac{\vec Pc}{E} = \frac{\vec V}{c}$, where $\vec V$ has the dimension of a velocity. The most natural possibility is that $\vec V$ is the velocity of the particle. So, we see, that a particule with speed $|\vec V|=c$ in a galilean frame, has also  $|\vec V'|=c$ in another galilean frame. We could also check that the quantity $E^2-\vec p^2c^2$ is conserved by Lorentz transformations, and call this quantity $m^2c^4$, where $m$ is the mass. So, particles who have $|\frac{\vec Pc}{E}| = |\frac{\vec V}{c}|=1$, are massless particles.