  We present a study of the two-photon-exchange ($2\gamma$-exchange)
corrections to the $S$-levels in muonic ($\mu$D) and ordinary (D) deuterium
within the pionless effective field theory (pionless EFT). Our calculation
proceeds up to next-to-next-to-next-to-leading order (N3LO) in the pionless EFT
expansion. The only unknown low-energy constant entering the calculation at
this order corresponds to the coupling of a longitudinal photon to the
nucleon-nucleon system. To minimise its correlation with the deuteron charge
radius, it is extracted using the information about the hydrogen-deuterium
isotope shift. We find the elastic $2\gamma$-exchange contribution in $\mu$D
larger by several standard deviations than obtained in other recent
calculations. This discrepancy ameliorates the mismatch between theory and
experiment on the size of $2\gamma$-exchange effects, and is attributed to the
properties of the deuteron elastic charge form factor parametrisation used to
evaluate the elastic contribution. We identify a correlation between the
deuteron charge and Friar radii, which can help one to judge how well a form
factor parametrisation describes the low-virtuality properties of the deuteron.
We also evaluate the higher-order $2\gamma$-exchange contributions in $\mu$D,
generated by the single-nucleon structure and expected to be the most important
terms beyond N3LO. The uncertainty of the theoretical result is dominated by
the truncation of the pionless EFT series and is quantified using a Bayesian
approach. The resulting extractions of the deuteron charge radius from the
$\mu$D Lamb shift, the $2S-1S$ transition in D, and the $2S-1S$
hydrogen-deuterium isotope shift, with the respective $2\gamma$-exchange
effects evaluated in a unified pionless EFT approach, are in perfect agreement.
