The Higgs pair production in gluon fusion is a sensitive probe of beyond-Standard Model (BSM) phenomena and its detection is a major goal for the LHC and higher energy hadron collider experiments. In this work we reanalyze the possible modifications of the Higgs pair production cross section within low energy supersymmetry models. We show that the supersymmetric contributions to the Higgs pair production cross section are strongly correlated with the ones of the single Higgs production in the gluon fusion channel. Motivated by the analysis of ATLAS and CMS Higgs production data, we show that the scalar superpartners contributions may lead to significant modification of the di-Higgs production rate and invariant mass distribution with respect to the SM predictions. We also analyze the combined effects on the di-Higgs production rate of a modification of the Higgs trilinear and top-quark Yukawa couplings in the presence of light stops. In particular, we show that due to the destructive interference of the triangle and box amplitude contributions to the di-Higgs production cross section, even a small modification of the top-quark Yukawa coupling can lead to a significant increase of the di-Higgs production rate.