We characterized MRI isocenter variation at various gantry positions in two 0.35 T MRgRT systems using two independent methods. First, image center-based quantification was employed on 3D volumetric and 2D cine images of a 24 cm diameter spherical phantom at various gantry positions in the MRI QA mode. The center of the phantom images was identified to quantify the variation of the imaging center at each gantry position. Second, image registration-based quantification was used in radiotherapy mode. 3D volumetric MRIs of a cylindrical phantom were acquired and corresponding image registration from MRI to planning CT was performed. The shifts of the couch were identified to quantify the variation of the imaging center. For verification of noticeable MRI isocenter variation, star-shot pattern measurements with five beams were delivered on the radio-chromic film inserted into the phantom after the couch was shifted. The center of the star-shot pattern was identified to quantify the variation of the imaging center. The proposed methods for measuring MRI isocenter variation were demonstrated with MR-LINAC and MR-60Co systems. Both of the MRgRT systems had field inhomogeneities <5 ppm over a 24 cm diameter spherical volume (DSV) and spatial integrity distortion: <1 mm within 100 mm radius and <2 mm within 175 mm radius. The MRI isocenter of the MR-LINAC system showed noticeable 3D variation (max magnitude: 1.8 mm) compared to that of MR-60Co system (max magnitude: 0.9 mm) relative to the reference gantry positions. In addition, 2D variations (max magnitude) of the MRI isocenter from sagittal cine images were 0.9 mm for the MR-LINAC system and 0.5 mm for the MR-60Co system. Two proposed methods quantified the MRI isocenter variation for various gantry positions in two 0.35 T MRgRT systems. The results of significant isocenter variation in the MR-LINAC system requires further investigation to determine the cause.