diff --git a/docs/src/whatsnew/latest.rst b/docs/src/whatsnew/latest.rst
index a6915c4170..76920c5ad5 100644
--- a/docs/src/whatsnew/latest.rst
+++ b/docs/src/whatsnew/latest.rst
@@ -39,6 +39,12 @@ This document explains the changes made to Iris for this release
    the :class:`~iris.cube.Cube` :attr:`~iris.cube.Cube.mesh_dim` (see
    :ref:`cube-statistics-collapsing`). (:issue:`5377`, :pull:`6003`)
 
+#. `@pp-mo`_ made a MeshCoord inherit a coordinate system from its location coord,
+   as it does its metadata.  N.B. mesh location coords can not however load a
+   coordinate system from netcdf at present, as this needs the 'extended'
+   grid-mappping syntax -- see : :issue:`3388`.
+   (:issue:`5562`, :pull:`6016`)
+
 
 🐛 Bugs Fixed
 =============
diff --git a/lib/iris/experimental/ugrid/mesh.py b/lib/iris/experimental/ugrid/mesh.py
index 2c07daa481..651906b748 100644
--- a/lib/iris/experimental/ugrid/mesh.py
+++ b/lib/iris/experimental/ugrid/mesh.py
@@ -2784,6 +2784,11 @@ def fix_repr(val):
                     )
                     raise ValueError(msg)
 
+        # Don't use 'coord_system' as a constructor arg, since for
+        # MeshCoords it is deduced from the mesh.
+        # (Otherwise a non-None coord_system breaks the 'copy' operation)
+        use_metadict.pop("coord_system")
+
         # Call parent constructor to handle the common constructor args.
         super().__init__(points, bounds=bounds, **use_metadict)
 
@@ -2809,8 +2814,28 @@ def axis(self):
 
     @property
     def coord_system(self):
-        """The coordinate-system of a MeshCoord is always 'None'."""
-        return None
+        """The coordinate-system of a MeshCoord.
+
+        It comes from the `related` location coordinate in the mesh.
+        """
+        # This matches where the coord metadata is drawn from.
+        # See : https://github.com/SciTools/iris/issues/4860
+        select_kwargs = {
+            f"include_{self.location}s": True,
+            "axis": self.axis,
+        }
+        try:
+            # NOTE: at present, a MeshCoord *always* references the relevant location
+            # coordinate in the mesh, from which its points are taken.
+            # However this might change in future ..
+            # see : https://github.com/SciTools/iris/discussions/4438#bounds-no-points
+            location_coord = self.mesh.coord(**select_kwargs)
+            coord_system = location_coord.coord_system
+        except CoordinateNotFoundError:
+            # No such coord : possible in UGRID, but probably not Iris (at present).
+            coord_system = None
+
+        return coord_system
 
     @coord_system.setter
     def coord_system(self, value):
diff --git a/lib/iris/tests/integration/experimental/test_meshcoord_coordsys.py b/lib/iris/tests/integration/experimental/test_meshcoord_coordsys.py
new file mode 100644
index 0000000000..353d56004f
--- /dev/null
+++ b/lib/iris/tests/integration/experimental/test_meshcoord_coordsys.py
@@ -0,0 +1,129 @@
+# Copyright Iris contributors
+#
+# This file is part of Iris and is released under the BSD license.
+# See LICENSE in the root of the repository for full licensing details.
+"""Integration tests for MeshCoord.coord_system() behaviour."""
+
+import pytest
+
+import iris
+from iris.coord_systems import GeogCS
+from iris.experimental.ugrid.load import PARSE_UGRID_ON_LOAD
+from iris.tests.stock.netcdf import ncgen_from_cdl
+
+TEST_CDL = """
+netcdf mesh_test {
+    dimensions:
+        node = 3 ;
+        face = 1 ;
+        vertex = 3 ;
+    variables:
+        int mesh ;
+            mesh:cf_role = "mesh_topology" ;
+            mesh:topology_dimension = 2 ;
+            mesh:node_coordinates = "node_x node_y" ;
+            mesh:face_node_connectivity = "face_nodes" ;
+        float node_x(node) ;
+            node_x:standard_name = "longitude" ;
+        float node_y(node) ;
+            node_y:standard_name = "latitude" ;
+        int face_nodes(face, vertex) ;
+            face_nodes:cf_role = "face_node_connectivity" ;
+            face_nodes:start_index = 0 ;
+        float node_data(node) ;
+            node_data:coordinates = "node_x node_y" ;
+            node_data:location = "node" ;
+            node_data:mesh = "mesh" ;
+}
+"""
+
+BASIC_CRS_STR = """
+        int crs ;
+          crs:grid_mapping_name = "latitude_longitude" ;
+          crs:semi_major_axis = 6371000.0 ;
+"""
+
+
+def find_i_line(lines, match):
+    (i_line,) = [i for i, line in enumerate(lines) if match in line]
+    return i_line
+
+
+def make_file(nc_path, node_x_crs=False, node_y_crs=False):
+    lines = TEST_CDL.split("\n")
+
+    includes = ["node_x"] if node_x_crs else []
+    includes += ["node_y"] if node_y_crs else []
+    includes = " ".join(includes)
+    if includes:
+        # insert a grid-mapping
+        i_line = find_i_line(lines, "variables:")
+        lines[i_line + 1 : i_line + 1] = BASIC_CRS_STR.split("\n")
+
+        i_line = find_i_line(lines, "float node_data(")
+        ref_lines = [
+            # NOTE: space to match the surrounding indent
+            f'            node_data:grid_mapping = "crs: {includes}" ;',
+            # NOTE: this is already present
+            # f'node_data:coordinates = "{includes}" ;'
+        ]
+        lines[i_line + 1 : i_line + 1] = ref_lines
+
+    cdl_str = "\n".join(lines)
+    ncgen_from_cdl(cdl_str=cdl_str, cdl_path=None, nc_path=nc_path)
+
+
+@pytest.mark.parametrize("cs_axes", ["--", "x-", "-y", "xy"])
+def test_default_mesh_cs(tmp_path, cs_axes):
+    """Test coord-systems of mesh cube and coords, if location coords have a crs."""
+    nc_path = tmp_path / "test_temp.nc"
+    do_x = "x" in cs_axes
+    do_y = "y" in cs_axes
+    make_file(nc_path, node_x_crs=do_x, node_y_crs=do_y)
+    with PARSE_UGRID_ON_LOAD.context():
+        cube = iris.load_cube(nc_path, "node_data")
+    meshco_x, meshco_y = [cube.coord(mesh_coords=True, axis=ax) for ax in ("x", "y")]
+    # NOTE: at present, none of these load with a coordinate system,
+    #  because we don't support the extended grid-mapping syntax.
+    #  see: https://github.com/SciTools/iris/issues/3388
+    assert meshco_x.coord_system is None
+    assert meshco_y.coord_system is None
+
+
+def test_assigned_mesh_cs(tmp_path):
+    # Check that when a coord system is manually assigned to a location coord,
+    # the corresponding meshcoord reports the same cs.
+    nc_path = tmp_path / "test_temp.nc"
+    make_file(nc_path)
+    with PARSE_UGRID_ON_LOAD.context():
+        cube = iris.load_cube(nc_path, "node_data")
+    nodeco_x = cube.mesh.coord(include_nodes=True, axis="x")
+    meshco_x, meshco_y = [cube.coord(axis=ax) for ax in ("x", "y")]
+    assert nodeco_x.coord_system is None
+    assert meshco_x.coord_system is None
+    assert meshco_y.coord_system is None
+    assigned_cs = GeogCS(1.0)
+    nodeco_x.coord_system = assigned_cs
+    assert meshco_x.coord_system is assigned_cs
+    assert meshco_y.coord_system is None
+    # This also affects cube.coord_system(), even though it is an auxcoord,
+    #  since there are no dim-coords, or any other coord with a c-s.
+    # TODO: this may be a mistake -- see https://github.com/SciTools/iris/issues/6051
+    assert cube.coord_system() is assigned_cs
+
+
+def test_meshcoord_coordsys_copy(tmp_path):
+    # Check that copying a meshcoord with a coord system works properly.
+    nc_path = tmp_path / "test_temp.nc"
+    make_file(nc_path)
+    with PARSE_UGRID_ON_LOAD.context():
+        cube = iris.load_cube(nc_path, "node_data")
+    node_coord = cube.mesh.coord(include_nodes=True, axis="x")
+    assigned_cs = GeogCS(1.0)
+    node_coord.coord_system = assigned_cs
+    mesh_coord = cube.coord(axis="x")
+    assert mesh_coord.coord_system is assigned_cs
+    meshco_copy = mesh_coord.copy()
+    assert meshco_copy == mesh_coord
+    # Note: still the same object, because it is derived from the same node_coord
+    assert meshco_copy.coord_system is assigned_cs