The probability of primordial black hole (PBH) formation is known to be
boosted during the Quantum Chromodynamics (QCD) crossover due to a slight
reduction of the equation of state. This induces a high peak and other features
in the PBH mass distribution. But the impact of this variation during the PBH
formation has been so far neglected. In this work we simulate for the first
time the formation of PBHs by taking into account the varying equation of state
at the QCD epoch, compute the over-density threshold using different curvature
profiles and find that the resulting PBH mass distributions are significantly
impacted. The expected merger rate distributions of early and late PBH binaries
is comparable to the ones inferred from the GWTC-3 catalog for dark matter
fractions in PBHs within $0.1<{f}_{\mathrm{P}\mathrm{B}\mathrm{H}}<1$ . The distribution of
gravitational-wave events estimated from the volume sensitivity could explain
mergers around $30-50{M}_{\odot}$ , with asymmetric masses like GW190814, or in the
pair-instability mass gap like GW190521. However, none of the considered cases
leads to a multi-modal distribution with a secondary peak around $8-15{M}_{\odot}$ , as suggested by the GWTC-3 catalog, possibly pointing to a mixed
population of astrophysical and primordial black holes.