One should, however, note that ##\epsilon_0## and ##\mu_0## are just arbitrary conversion constants, which are introduced into physics just for practical convenience within the international system of units. The reason is that you get convenient numbers for electrical engineering. From a theoretical-physics point of view this complicates the description of electromagnetic phenomena very much, and that's why in high-energy physics one uses more convenient units, the socalled Heaviside-Lorentz units, which are rationalized Gaussian units. The only fundamental constant entering the equations of electromagnetism is then the speed of light in vacuum, ##c##, which however in view of special relativity is also only a conversion constant between measurements of time (in seconds) and length (in meters), and that's in fact how it is defined in the SI. This is also quite inconvenient for theoretical purposes and that's why one also introduces units, where ##c=1##. Since HEP physics also deals with quantum phenomena one also sets ##\hbar=1## since also ##\hbar## is nothing than a conversion constant between the arbitrary man-made SI units and natural units.
Now it's an empirical fact that photons are massless, and since we can observe free photons, indeed one can measure the speed of light as the propagation speed of electromagnetic waves in a vacuum.
The case of gluons is much more subtle. While photons are uncharged, the gluons carry color charge, i.e., the non-Abelian charge of the strong interaction. There are 8 differently charged gluons (quarks come in three color charges, usually dubbed "red, green, blue" which is why the entire theory of strong interactions is called quantum chromodynamics, QCD). Since we can treat quantum field theories as QED and QCD analytically only in a perturbative way, the case of QCD is very difficult: In contradistinction to QED with its uncharged gauge field (the electromagnetic field with the asymptotic free single-particle excitations called photons) the gauge field of QCD is charged, and as the perturbative analysis shows, this causes the theory to be asymptotically free. In QFT we always have to renormalize, and this renormalization procedure leads to the introduction of a momentum scale, and the coupling constants become dependent on this scale. In QCD this "running coupling" becomes small at large momentum scales, i.e., only when quarks and gluons scatter with a large momentum transfer you can use perturbation theory to evaluate cross sections.
However, this is problematic either since there's also "confinement", which cannot be understood within the realm of perturbation theory. Perturbation theory only gives a glimpse: At low momentum scales the coupling becomes large, and perturbation theory breaks down. Phenomenologically we have never seen free quarks and gluons, but the observable "asymptotic free degrees of freedom" involving the strong interaction are hadrons, i.e., bound states of quarks and gluons which all carry net-color charge of 0. The standard hadrons are the mesons (effectively describable by as bound states of a quark and an antiquark) and baryons (effectively describable as bound states of three quarks) and their antiparticles. This is known as confinement: The quarks and gluons do not occur as asymptotic free particles and thus you cannot observe the speed of gluons.
One way to see that QCD also describes this confinement correctly is to use lattice-gauge theory, which evaluates certain quantities of QCD with computers on a discrete space-time lattice. Among other things what lattice QCD can calculate well is the mass spectrum of hadrons, and this is done with quite some accuracy, leading to the correct masses of the observed hadrons (+ the prediction of some more not yet observed hadrons).
