  We perform a calculation of the glueball spectrum for $N_f=4$ degenerate
dynamical fermions with masses corresponding to light pions. We do so by making
use of ensembles produced within the framework of maximally twisted fermions by
the Extended Twisted Mass Collaboration (ETMC). We obtain masses of states that
fall into the irreducible representations of the octahedral group of rotations
in combination with the quantum numbers of charge conjugation $C$ and parity
$P$; the above quantum numbers result in 20 distinct irreducible
representations. We implement the Generalized Eigenvalue Problem (GEVP) using a
basis that consists only of gluonic operators. The purpose of this work is to
investigate the effect of light dynamical quarks on the glueball spectrum and
how this compares to the statistically more accurate spectrum of $SU(3)$ pure
gauge theory. Given that glueball states may have broad widths and thus need to
be disentangled from all the relevant mixings, we use large ensembles of the
order of ${\sim {~\cal O}}(20 {\rm K})$ configurations. Despite the large
ensembles, the statistical uncertainties allow us to extract the masses for
only a few irreducible representations; namely $A_1^{++}$, $A_1^{-+}$, $E^{++}$
as well as $T_2^{++}$. The results for the scalar $A_1^{++}$ representation
show that an additional state appears as the lightest state in the scalar
$A_1^{++}$ channel of the glueball spectrum, while the next two excited states
are consistent with the lightest two states of the pure gauge theory. To
further elucidate the nature of this additional state we perform a calculation
using $N_f=2+1+1$ configurations and this demonstrates that it possesses a
large quark content. Finally, the ground states of the $E^{++}$ and $T_2^{++}$
tensor channels and of the $A_1^{-+}$ pseudoscalar channel show, at most, minor
effects due to the inclusion of dynamical quarks.
