diff --git a/doc/changes/latest.inc b/doc/changes/latest.inc
index 85ac44237a5..ca4557a04a4 100644
--- a/doc/changes/latest.inc
+++ b/doc/changes/latest.inc
@@ -198,6 +198,8 @@ Bug
 
 - Fix bug with :func:`mne.setup_volume_source_space` when ``volume_label`` was supplied where voxels slightly (in a worst case, about 37% times ``pos`` in distance) outside the voxel-grid-based bounds of regions were errantly included, by `Eric Larson`_
 
+- Fix bug with :ref:`mne coreg` where reverse scalings were applied to ``src[0]['src_mri_t']`` for volume source spaces, so morphing and plotting did not work correctly by `Eric Larson`_
+
 - Fix bug with :func:`mne.io.read_raw_ctf` when reference magnetometers have the compensation grade marked by `Eric Larson`_
 
 - Fix bug with `mne.SourceSpaces.export_volume` with ``use_lut=False`` where no values were written by `Eric Larson`_
diff --git a/mne/coreg.py b/mne/coreg.py
index 872c9afeed6..002f566ea62 100644
--- a/mne/coreg.py
+++ b/mne/coreg.py
@@ -1091,10 +1091,10 @@ def scale_source_space(subject_to, src_name, subject_from=None, scale=None,
         ss['subject_his_id'] = subject_to
         ss['rr'] *= scale
         # additional tags for volume source spaces
-        if 'vox_mri_t' in ss:
+        for key in ('vox_mri_t', 'src_mri_t'):
             # maintain transform to original MRI volume ss['mri_volume_name']
-            ss['vox_mri_t']['trans'][:3, :3] /= scale
-            ss['src_mri_t']['trans'][:3, :3] /= scale
+            if key in ss:
+                ss[key]['trans'][:3] *= scale[:, np.newaxis]
         # distances and patch info
         if uniform:
             if ss['dist'] is not None:
diff --git a/mne/tests/test_coreg.py b/mne/tests/test_coreg.py
index 53647ecb0c2..c1c99769d44 100644
--- a/mne/tests/test_coreg.py
+++ b/mne/tests/test_coreg.py
@@ -113,10 +113,13 @@ def test_scale_mri(tmpdir, few_surfaces, scale):
     assert os.path.isfile(os.path.join(tempdir, 'flachkopf', 'surf',
                                        'lh.sphere.reg'))
     vsrc_s = mne.read_source_spaces(spath % 'vol-50')
-    pt = np.array([0.12, 0.41, -0.22])
-    assert_array_almost_equal(
-        apply_trans(vsrc_s[0]['src_mri_t'], pt * np.array(scale)),
-        apply_trans(vsrc[0]['src_mri_t'], pt))
+    for vox in ([0, 0, 0], [1, 0, 0], [0, 1, 0], [0, 0, 1], [1, 2, 3]):
+        idx = np.ravel_multi_index(vox, vsrc[0]['shape'], order='F')
+        err_msg = f'idx={idx} @ {vox}, scale={scale}'
+        assert_allclose(apply_trans(vsrc[0]['src_mri_t'], vox),
+                        vsrc[0]['rr'][idx], err_msg=err_msg)
+        assert_allclose(apply_trans(vsrc_s[0]['src_mri_t'], vox),
+                        vsrc_s[0]['rr'][idx], err_msg=err_msg)
     scale_labels('flachkopf', subjects_dir=tempdir)
 
     # add distances to source space after hacking the properties to make
diff --git a/mne/viz/_brain/_brain.py b/mne/viz/_brain/_brain.py
index 49afb196ac6..8eb8567f86a 100644
--- a/mne/viz/_brain/_brain.py
+++ b/mne/viz/_brain/_brain.py
@@ -705,9 +705,10 @@ def _add_volume_data(self, hemi, src, volume_options):
             vertices = self._data[hemi]['vertices']
             assert self._data[hemi]['array'].shape[0] == len(vertices)
             # MNE constructs the source space on a uniform grid in MRI space,
-            # but let's make sure
+            # but mne coreg can change it to be non-uniform, so we need to
+            # use all three elements here
             assert np.allclose(
-                src_mri_t[:3, :3], np.diag([src_mri_t[0, 0]] * 3))
+                src_mri_t[:3, :3], np.diag(np.diag(src_mri_t)[:3]))
             spacing = np.diag(src_mri_t)[:3]
             origin = src_mri_t[:3, 3] - spacing / 2.
             scalars = np.zeros(np.prod(dimensions))