We consider equidistant approximations of stochastic integrals driven by Hölder continuous Gaussian processes of order H>12\\documentclass[12pt]{minimal} \\usepackage{amsmath} \\usepackage{wasysym} \\usepackage{amsfonts} \\usepackage{amssymb} \\usepackage{amsbsy} \\usepackage{mathrsfs} \\usepackage{upgreek} \\setlength{\\oddsidemargin}{-69pt} \\begin{document}$$H>\\frac{1}{2}$$\\end{document} with discontinuous integrands involving bounded variation functions. We give exact rate of convergence in the L1\\documentclass[12pt]{minimal} \\usepackage{amsmath} \\usepackage{wasysym} \\usepackage{amsfonts} \\usepackage{amssymb} \\usepackage{amsbsy} \\usepackage{mathrsfs} \\usepackage{upgreek} \\setlength{\\oddsidemargin}{-69pt} \\begin{document}$$L^1$$\\end{document}-distance and provide examples with different drivers. It turns out that the exact rate of convergence is proportional to n1-2H\\documentclass[12pt]{minimal} \\usepackage{amsmath} \\usepackage{wasysym} \\usepackage{amsfonts} \\usepackage{amssymb} \\usepackage{amsbsy} \\usepackage{mathrsfs} \\usepackage{upgreek} \\setlength{\\oddsidemargin}{-69pt} \\begin{document}$$n^{1-2H}$$\\end{document}, which is twice as good as the best known results in the case of discontinuous integrands and corresponds to the known rate in the case of smooth integrands. The novelty of our approach is that, instead of using multiplicative estimates for the integrals involved, we apply a change of variables formula together with some facts on convex functions allowing us to compute expectations explicitly.