In relativistic nuclear collisions the production of hadrons with light (u,d,s) quarks is quantitatively described in the framework of the Statistical Hadronization Model (SHM). Since charm quarks are dominantly produced in initial hard collisions but interact strongly in the hot fireball, charmed hadrons can be incorporated into the SHM by treating charm quarks as 'impurities' with thermal distributions, with the total charm content of the fireball fixed by the measured open charm cross section. We demonstrate that this way the measured multiplicities of single charm hadrons in Pb-Pb collisions at LHC energies can be well described with the same thermal parameters as for (u,d,s) hadrons. Furthermore, transverse momentum distributions are computed in a hydrodynamic approach also incorporating resonance decays. The approach is extended to lighter collision systems down to O-O and includes doubly- and triply-charmed hadrons. We show predictions for production probabilities of such states exhibiting a characteristic and rather spectacular enhancement hierarchy.