Heart rate constantly varies under physiological conditions, termed heart rate variability (HRV), and in clinical studies, low HRV is associated with a greater risk of cardiac arrhythmias. Prior work has shown that HRV influences the temporal patterns of electrical activity, specifically the formation of pro-arrhythmic alternans, a beat-to-beat alternation in the action potential duration (APD), or intracellular calcium (Ca) levels. We previously showed that HRV may be anti-arrhythmic by disrupting APD and Ca alternations in a homogeneous cardiac myocyte. Here, we expand on our previous work, incorporating variation in subcellular Ca handling (also known to influence alternans) into a nonlinear map model of a cardiac myocyte composed of diffusively coupled Ca release units (CRUs). Ca-related parameters and initial conditions of each CRU are varied to mimic subcellular Ca heterogeneity, and a stochastic pacing sequence reproduces HRV. We find that subcellular Ca heterogeneity promotes the formation of spatially discordant subcellular alternans patterns, which decreases whole cell Ca and APD alternation for low and moderate HRV, while high subcellular Ca heterogeneity and HRV both promote electromechanical desynchronization. Finally, we find that for low and moderate HRV, both the specific subcellular Ca-related parameters and the pacing sequences influence measures of electromechanical dynamics, while for high HRV, these measures depend predominantly on the pacing sequence. Our results suggest that pro-arrhythmic subcellular discordant alternans tend to form for low levels of HRV, while high HRV may be anti-arrhythmic due to mitigated influence from subcellular Ca heterogeneity and desynchronization of APD from Ca instabilities.