We formulate the high temperature expansion in supersymmetric matrix quantum mechanics with 4, 8 and 16 supercharges. The models can be obtained by dimensionally reducing = 1 U(N) super Yang-Mills theory in D = 4,6,10 to 1 dimension, respectively. While the non-zero frequency modes become weakly coupled at high temperature, the zero modes remain strongly coupled. We find, however, that the integration over the zero modes that remains after integrating out all the non-zero modes perturbatively, reduces to the evaluation of connected Green's functions in the bosonic IKKT model. We perform Monte Carlo simulation to compute these Green's functions, which are then used to obtain the coefficients of the high temperature expansion for various quantities up to the next-leading order. Our results nicely reproduce the asymptotic behaviors of the recent simulation results at finite temperature. In particular, the fermionic matrices, which decouple at the leading order, give rise to substantial effects at the next-leading order, reflecting finite temperature behaviors qualitatively different from the corresponding models without fermions.