A search for dark matter is performed looking for events with large missing transverse momentum and a Higgs boson decaying either to a pair of bottom quarks or to a pair of photons. The data from proton-proton collisions at a center-of-mass energy of 13 TeV, collected in 2015 with the CMS detector at the LHC, correspond to an integrated luminosity of 2.3 fb−1. Results are interpreted in the context of a Z′-two-Higgs-doublet model, where the gauge symmetry of the standard model is extended by a U(1)Z ′ group, with a new massive Z′ gauge boson, and the Higgs sector is extended with four additional Higgs bosons. In this model, a high-mass resonance Z′ decays into a pseudoscalar boson A and a light SM-like scalar Higgs boson, and the A decays to a pair of dark matter particles. No significant excesses are observed over the background prediction. Combining results from the two decay channels yields exclusion limits in the signal cross section in the m Z ′ - m A phase space. For example, the observed data exclude the Z′ mass range from 600 to 1860 GeV, for Z′ coupling strength g Z ′ = 0.8, the coupling of A with dark matter particles g χ = 1, the ratio of the vacuum expectation values tan β = 1, and m A = 300 GeV. The results of this analysis are valid for any dark matter particle mass below 100 GeV.