We perform a Bayesian analysis of the maximum mass of neutron stars with a quark core, incorporating the observational data from tidal deformability of the GW170817 binary neutron star merger as detected by LIGO/Virgo and the mass and radius of PSR J0030+0451 as detected by \nicer. The analysis is performed under the assumption that the hadron-quark phase transition is of first order, where the low-density hadronic matter described in a unified manner by the soft QMF or the stiff DD2 equation of state (EOS) transforms into a high-density phase of quark matter modeled by the generic “Constant-sound-speed” (CSS) parameterization. The mass distribution measured for the pulsar, MSP J0740+6620, is used as the lower limit on . We find the most probable values of the hybrid star maximum mass are () for QMF (DD2), with an absolute upper bound around , to the 90% posterior credible level. Such results appear robust with respect to the uncertainties in the hadronic EOS. We also discuss astrophysical implications of this result, especially on the post-merger product of GW170817, short gamma-ray bursts, and other likely binary neutron star mergers.