  We compute the gamma-ray output of axion-mediated dark matter and derive the
corresponding constraints set by recent data. In such scenarios the dark matter
candidate is a Dirac fermion that pair-annihilates into axions and/or scalars.
Provided that the axion decays (at least partly) into photons, these models
naturally give rise to a box-shaped gamma-ray spectrum that may present two
distinct phenomenological behaviours: a narrow box, resembling a line at half
the dark matter mass, or a wide box, spanning an extensive energy range up to
the dark matter mass. Remarkably, we find that in both cases a sizable
gamma-ray flux is predicted for a thermal relic without fine-tuning the model
parameters nor invoking boost factors. This large output is in line with recent
Fermi-LAT observations towards the Galactic centre region and is on the verge
of being excluded. We then make use of the Fermi-LAT and H.E.S.S. data to
derive robust, model-independent upper limits on the dark matter annihilation
cross section for the narrow and wide box scenarios. H.E.S.S. constraints, in
particular, turn out to match the ones from Fermi-LAT at hundreds of GeV and
extend to multi-TeV masses. Future Cherenkov telescopes will likely probe
gamma-ray boxes from thermal dark matter relics in the whole multi-TeV range, a
region hardly accessible to direct detection, collider searches and other
indirect detection strategies.
