We present the mixed QCD-electroweak corrections to Higgs boson pair production in the quark-antiquark channel. The virtual amplitudes are computed fully analytically using the method of differential equations. We determine the integration constants by matching our expressions to the large mass expansion limit of the canonical integrals. We implement the results in the POWHEG-BOX framework for phenomenological studies. The corrections are found to have a significant impact on the shapes of differential cross sections, reaching up to +10 % for the invariant mass distribution of the Higgs boson pair near the production threshold. This channel has not been considered before in calculations of the next-to-leading order electroweak corrections to Higgs boson pair production.