Why are gluons believed to be massless?

Simply put, the Higgs isn't charged under the strong force. It doesn't have standard electrical charge either. The $W^{\pm},Z$ bosons aquire a mass through the Higgs mechanism because the Higgs itself is charged under the weak force. Leptons aquire masses through the Higgs mechanism because they too interact with the Higgs.

No Higgs interaction means no effective mass.

You ask why the Higgs does not interact with gluons. It has to do with the quantum numbers (charges) of the fundamental particles in the standard model. It turns out that you aren't allowed to freely choose quantum numbers for the different particles. If you made a bad choice, you'd violate gauge invariance and have an inconsistent theory. This puts relatively strict constraints on the allowed quantum numbers. check out Gauge Anomaly for more details.

Basically, the known quantum numbers of the other standard model particles constrains the allowed quantum numbers of the Higgs, specifically prohibiting a gluon-higgs interaction. if you wanted to add that interaction, you would necessarily imply the existence of other particles in order to balances all the charges of the theory. I don't know if that would be possible, but its a matter of simple algebra to figure it out.

One has to keep clear in mind that the standard model, which has the SU(3)xSU(2)xU(1) structure is called standard because it describes extremely well an enormous number of particle data and is successful in its predictions for new experiments. As Luc J. Bourhis explores in his answer this does not stop theorists from exploring beyond the standard model theories.

The standard model is a quantum field theoretical model, which means it has precise predictions for evaluating Feynman diagrams which will be used for fitting or predicting experimental data. In its mathematical structure there is no simple vertex, i.e. lowest order, ( higher order loops can give an interaction) between the electroweak bosons (W, Z, γ) of the table, and the gluon .

Thus by construction of the mathematics of the model the gluon does not "see" the Higgs field. As the very existence of the concept of a gluon depends on the mathematics of the standard model, our "belief" in the standard model means massless gluons. Other theories exist beyond the standard model , which may give a coupling, ( example ) but they are beyond the standard model. Be sure that experiments will be testing for any discrepancy with the standard model that could be due to a massive gluon.

Question in comment:

Could you please elaborate how you jumped from the weak bosons to the Higgs field

It has to do with the fact that Feynman diagrams are formulated with specific rules on the fields of all the elementary particles in the table, including the Higgs field. So in a simple diagram of e-e- scattering,

the electron creation/annihilation operators operate sequentially on the electron field, when a created electron interacts with the photon field and photon creation operators create a virtual photon, which interacts with the electron field and creates an outgoing electron.

Thus in field theory the vertex coupling constant should exist with the field so that the interaction can happen and, in this case, a photon can be generated.

Because of the zero mass there is no coupling constant for the vertex gluon-higgs field, to generate a virtual higgs and continue from there. Only higher quark loops can act , i.e. glue to q q_bar, virtual higgs loop Higgs meson from gluon fusion.

Edit after edit of OP:

This is ok for the Higgs boson but not for the Higgs field

The Higgs has a weak charge and therefore interacts with the W and Z bosons thus giving them mass. The Higgs does not have the electric or color charge and therefore does not interact with the photon or gluons thus leaving them massless.

The mass the elementary particles of the SM table acquire comes from the interaction with the Higgs FIELD, not the Higgs Boson. The Higgs boson is just another massive particle in the elementary particles table:

These are the only masses the Higgs Field generates. including the Higgs boson itself. All other masses come from special relativity relations, the invariant mass of the added fourvectors of composite particles.

The standard model has a Lagrangian describing the SU(3)xSU(2)xU(1) experimentaly observed symmetries in particle data and any comparison depends on feynman diagram calculations within this model, where there are strict rules for the exchange vertices.

Each elementary particle in the table defines a field in four dimensional space, and the particle is considered and excitation on this field. The electron is an excitation of the electron field, the Higgs boson is the exitation of the Higgs field.

The elementary particles in the table do not acquire mass by exchanging anything in the form of feynman diagrams. The mass is acquired once, at the time of symmetry breaking of the electroweak interactions, where the three couplings approach each other :

when the couplings of the weak and electromagnetic change due to the Higgs field. One has to study mathematically this to be convinced, but the fact is , that the standard model as now known describes practically all particle physics data and is very predictive of new ones, as the LHC experience shows.

The gluons are massless by construction the way the photon is, and the way the Z and W were massless before symmetry breaking..

Once again, it is the Higgs Field that gives the masses to the elementary particles , not the boson.

I disagree with @Mr.Weather. Eventually all the theoretical arguments don't hold much water. We could very well have a very small gluon mass which would badly break the mathematical beauty of the of the Standard Model but with little enough consequences that we haven't noticed yet. The ultimate judge is experiment. Increasingly precise experiments. Let's look at how we can give gluons a mass then.

0. Caveat: fundamental mass vs effective mass

Gluon, as quarks, do acquire an effective mass through the dynamic of QCD: naively they dress up with clouds of gluons. Taking into account such a gluon mass, which is of the order of 1 GeV, has proven important to explain experimental results for the decay $J/Psi\to gg\gamma$ and similar decays of $\Upsilon$. This is different from a fundamental mass, which appear at the level of the Lagrangian of the theory (and actually this effective mass depends on the particular process under scrutiny). You are clearly asking about the latter but I wanted to clear this up.

1. Bare mass

Before I address your question about the Higgs mechanism, I need to mention that we could add a mass term to QCD "by hand". Theorists despise this solution because we break local gauge symmetry and renormalisability but again Nature would be the ultimate judge on that. But in all fairness, renormalisability is a big deal. Theorists do a lot of perturbative computations. Let's say the coupling constant of a theory is $\alpha$. One can compute terms of order $\alpha$, then $\alpha^2$, $\alpha^3$, etc with increasingly much more work involved but that's another story: this is perturbative computation. In a non-renormalisable theory, typically, each order will introduce new constants which have to be fitted to experimental data. And yes, that means that eventually, there are an infinite number of unknown constants in the theory. This is not so in renormalisable theory: there is only a finite, and small, number of constants which we know about upfront and that does not change as we move up the perturbative ladder order by order. Thus a renormalisable theory has typically vastly more predictive power. So yes, the Standard Model has all the couplings with the Higgs which are unknown, and the mixing matrices between quarks and neutrino flavours which are also unknown but at least we have them all, once and for all.

2. Coloured Higgs

The Higgs mechanism can be tweaked to give a mass to gluon. There are actually quite a few possibilities. They are far too technical for me to try to explain. Suffice to say that these solutions have more than one Higgs field which carry a colour number. Some of them predicts both massive and massless gluons. Such models are a priori more appealing because they preserve renormalisability.

One of the problem with such models is the following: the so-called running of QCD coupling constant, denoted $\alpha_S$. In simple terms, $\alpha_S$ depends on the energy scale $\mu$ of the process: it grows as $\mu$ decreases. This explains why the interaction is so strong at low energy in bound states like mesons and baryons, and why it is quite weak at high energy, at LHC energies e.g., a phenomenon called asymptotic freedom (i.e. the quarks become increasingly free as the energy scale increases). How much it decreases depends on the number of particles that can feel the strong interaction. It is very well verified experimentally that the variation of $\alpha_S$ is compatible with 5 quarks above the mass of the quark $b$ (the top quark decays too fast to play any role) and 4 quarks between that of the quark $c$ and $b$. But now we are adding Higgs bosons to that count. That rules out light Higgs. Worse, the coupling of Higgs with itself grows with the energy scale, and since quarks can interact through those Higgs now, it could easily break the very asymptotic freedom of QCD, which is fantastically well tested.

3. Composite Higgs

There are also models where the Higgs is a composite. A condensate of fermions which feels another gauge interaction than those in the Standard Model. So the gauge symmetry would be the Standard Model $SU(3)_c\times SU(2)_L\times U(1)_Y$ but then $\times SU(2)_{c'}$ and there would be two new fermions which would be single for $SU(2)_L$ but at least of them couple to gluon like normal quarks, and couple to the new gauge interaction of course. One example can be found in [2], where it is concluded that it predicts too high a gluon mass actually.

4. Experiment

In any case, measurements have the final say. There are a couple of effects that have been looked for. First, with a mass $m_g$ for the gluon, the potential for the strong interaction becomes essentially 0 for distances $r>1/m_g$ whereas below that it grows linearly with $r$. There is therefore a potential barrier at an energy $E_*$ proportional to $m_g$. With collisions energetic enough, we should therefore be able to get the quarks in protons go over that barrier. Or the quarks could tunnel through the barrier, making the proton unstable. The final argument is cosmological: in the early big bang, quarks are roaming free. When energy goes below $E_*$, some quarks will get confined but some will stay outside of the potential barriers, if they were farther than $1/m_g$ from another one. So we should have free quark relics. I am quoting an article from Yndurain from memory here [1], which is the one quoted in the Particle Data Group reviews. It is quite old and it predates the theoretical ideas I developed in 2) and 3) above. A quick look at the preprint database arXiv turns [2] which sets itself the goal to update on [1] actually, and which does review the models I wrote about in 2).

[1] F.J. Ynduráin. Limits on the mass of the gluon. Physics Letters B, 345(4):524 – 526, 1995.

[2] https://arxiv.org/abs/1005.0850v1