The study of Higgs boson properties is one of the main tasks of contemporary high-energy physics. Higgs boson interaction with gluons is interesting since it is a dominant production mechanism at LHC and since it is sensitive to unknown additional elementary particles. One of the major sources of uncertainty in the theoretical description of $ggH$ coupling originates from mixed QCD-electroweak contributions. NLO QCD corrections to these contributions are currently known only in the limit of large mass for electroweak bosons, and an exact computation is desirable. I will present the analytic evaluation of NLO mixed QCD-electroweak three-loop virtual corrections to $gg \to H$ in the framework of differential equations for master integrals.