diff --git a/package/MDAnalysis/coordinates/TPR.py b/package/MDAnalysis/coordinates/TPR.py
index 1940da4c92..a7f1d71dbe 100644
--- a/package/MDAnalysis/coordinates/TPR.py
+++ b/package/MDAnalysis/coordinates/TPR.py
@@ -127,8 +127,12 @@ def _read_first_frame(self):
             self.ts._pos = np.asarray(
                 tpr_utils.ndo_rvec(data, th.natoms), dtype=np.float32
             )
+            if self.convert_units:
+                self.convert_pos_from_native(self.ts._pos)
         if th.bV:
             self.ts.velocities = np.asarray(
                 tpr_utils.ndo_rvec(data, th.natoms), dtype=np.float32
             )
+            if self.convert_units:
+                self.convert_velocities_from_native(self.ts.velocities)
             self.ts.has_velocities = True
diff --git a/testsuite/MDAnalysisTests/coordinates/test_tpr.py b/testsuite/MDAnalysisTests/coordinates/test_tpr.py
index 8f27c3fdea..fb0aaa0f94 100644
--- a/testsuite/MDAnalysisTests/coordinates/test_tpr.py
+++ b/testsuite/MDAnalysisTests/coordinates/test_tpr.py
@@ -70,6 +70,7 @@
     TPR_gh_5145,
 )
 import MDAnalysis as mda
+from MDAnalysis.coordinates.TPR import TPRReader
 
 
 import pytest
@@ -80,6 +81,9 @@
 @pytest.mark.parametrize(
     "tpr_file, exp_first_atom, exp_last_atom, exp_shape, exp_vel_first_atom, exp_vel_last_atom",
     [
+        # NOTE: expected values are expressed in nm
+        # units and converted to Angstrom in the body
+        # of the test below, before assertions
         # this case is an alanine dipeptide
         # with neural network potential active
         # and nonzero velocities
@@ -460,11 +464,15 @@ def test_basic_read_tpr(
         u = mda.Universe(tpr_file, tpr_file)
     else:
         u = mda.Universe(tpr_file)
-    assert_allclose(u.atoms.positions[0, ...], exp_first_atom)
-    assert_allclose(u.atoms.positions[-1, ...], exp_last_atom)
+    assert_allclose(u.atoms.positions[0, ...], np.asarray(exp_first_atom) * 10)
+    assert_allclose(u.atoms.positions[-1, ...], np.asarray(exp_last_atom) * 10)
     assert_equal(u.atoms.positions.shape, exp_shape)
-    assert_allclose(u.atoms.velocities[0, ...], exp_vel_first_atom)
-    assert_allclose(u.atoms.velocities[-1, ...], exp_vel_last_atom)
+    assert_allclose(
+        u.atoms.velocities[0, ...], np.asarray(exp_vel_first_atom) * 10
+    )
+    assert_allclose(
+        u.atoms.velocities[-1, ...], np.asarray(exp_vel_last_atom) * 10
+    )
     assert_equal(u.atoms.velocities.shape, exp_shape)
 
 
@@ -483,8 +491,8 @@ def test_different_versions():
     # reading topology and positions/velocities
     # for the same system with different TPR
     # file versions
-    exp_first_atom = [3.25000e-01, 1.00400e00, 1.03800e00]
-    exp_last_atom = [-2.56000e-01, 1.37300e00, 3.59800e00]
+    exp_first_atom = [3.25000, 10.0400, 10.3800]
+    exp_last_atom = [-2.56000, 13.7300, 35.9800]
     exp_shape = (2263, 3)
     exp_vel = np.zeros(3)
     u = mda.Universe(TPR2020, TPR2024_4)
@@ -494,3 +502,29 @@ def test_different_versions():
     assert_allclose(u.atoms.velocities[-1, ...], exp_vel)
     assert_equal(u.atoms.positions.shape, exp_shape)
     assert_equal(u.atoms.velocities.shape, exp_shape)
+
+
+@pytest.mark.parametrize("convert_units", [True, False])
+def test_unit_swapping(convert_units):
+    reader = TPRReader(TPR_xvf_2024_4, convert_units=convert_units)
+    ts = reader.ts
+    actual_positions = ts.positions
+    actual_velocities = ts.velocities
+    # nm units:
+    expected_first_pos = np.asarray([3.19900e00, 1.62970e00, 1.54480e00])
+    expected_last_pos = np.asarray([3.39350e00, 3.49420e00, 3.02400e00])
+    expected_first_vel = np.asarray([-0.20668714, 0.26678202, -0.10564042])
+    expected_last_vel = np.asarray(
+        [-3.38010e-02, -3.22064e-01, -1.9863836e-01]
+    )
+    if convert_units:
+        # convert_units=True is the MDA default for readers, so we
+        # expect the nm to get converted to A
+        factor = 10
+    else:
+        # GROMACS "native" nm units are read in
+        factor = 1
+    assert_allclose(actual_positions[0, ...], expected_first_pos * factor)
+    assert_allclose(actual_positions[-1, ...], expected_last_pos * factor)
+    assert_allclose(actual_velocities[0, ...], expected_first_vel * factor)
+    assert_allclose(actual_velocities[-1, ...], expected_last_vel * factor)