diff --git a/NEWS.rst b/NEWS.rst
index d0ed1d1184..4a9e39a24f 100644
--- a/NEWS.rst
+++ b/NEWS.rst
@@ -9,6 +9,25 @@ Upgrade Guide
 Building openPMD-api now requires a compiler that supports C++14 or newer.
 Supported Python version are now 3.6 to 3.9.
 
+Python
+^^^^^^
+
+Reading the ``data_order`` of a mesh was broken.
+The old setter function (``set_data_order``) and read-only property (``data_order``) are now unified in a single, writable property:
+
+.. code-block:: python3
+
+   import openpmd_api as io
+
+   series = io.Series("data%T.h5", io.Access.read_only)
+   rho = series.iterations[0].meshes["rho"]
+   rho.data_order = 'C'  # or 'F'
+
+   print(rho.data_order == 'C')  # True
+
+Note: we recommend using ``'C'`` order since version 2 of the openPMD-standard will simplify this option to ``'C'``, too.
+For Fortran-ordered indices, please just invert the attributes ``axis_labels``, ``grid_spacing`` and ``grid_global_offset`` accordingly.
+
 
 0.12.0-alpha
 ------------
diff --git a/examples/3a_write_thetaMode_serial.py b/examples/3a_write_thetaMode_serial.py
index dcc6d10c97..d9c32c0f79 100755
--- a/examples/3a_write_thetaMode_serial.py
+++ b/examples/3a_write_thetaMode_serial.py
@@ -34,11 +34,11 @@
     E = series.iterations[0].meshes["E"]
     E.set_geometry(io.Geometry.thetaMode)
     E.set_geometry_parameters(geometry_parameters)
-    E.set_data_order(io.Data_Order.C)
     E.set_grid_spacing([1.0, 1.0])
     E.set_grid_global_offset([0.0, 0.0])
     E.set_grid_unit_SI(1.0)
     E.set_axis_labels(["r", "z"])
+    E.data_order = "C"
     E.unit_dimension = {io.Unit_Dimension.I: 1.0,
                         io.Unit_Dimension.J: 2.0}
 
diff --git a/src/binding/python/Mesh.cpp b/src/binding/python/Mesh.cpp
index bb761e8df2..a611dd1b79 100644
--- a/src/binding/python/Mesh.cpp
+++ b/src/binding/python/Mesh.cpp
@@ -51,8 +51,11 @@ void init_Mesh(py::module &m) {
         .def("set_geometry", &Mesh::setGeometry)
         .def_property_readonly("geometry_parameters", &Mesh::geometryParameters)
         .def("set_geometry_parameters", &Mesh::setGeometryParameters)
-        .def_property_readonly("data_order", &Mesh::dataOrder)
-        .def("set_data_order", &Mesh::setDataOrder)
+        .def_property("data_order",
+              [](Mesh const & mesh){ return static_cast< char >(mesh.dataOrder()); },
+              [](Mesh & mesh, char d){ mesh.setDataOrder(Mesh::DataOrder(d)); },
+              "Data Order of the Mesh (deprecated and set to C in openPMD 2)"
+        )
         .def_property_readonly("axis_labels", &Mesh::axisLabels)
         .def("set_axis_labels", &Mesh::setAxisLabels)
         .def_property_readonly("grid_spacing", &Mesh::gridSpacing<float>)
@@ -80,9 +83,4 @@ void init_Mesh(py::module &m) {
         .value("cylindrical", Mesh::Geometry::cylindrical)
         .value("spherical", Mesh::Geometry::spherical)
     ;
-
-    py::enum_<Mesh::DataOrder>(m, "Data_Order")
-        .value("C", Mesh::DataOrder::C)
-        .value("F", Mesh::DataOrder::F)
-    ;
 }
diff --git a/test/python/unittest/API/APITest.py b/test/python/unittest/API/APITest.py
index b039a73b11..35df720683 100644
--- a/test/python/unittest/API/APITest.py
+++ b/test/python/unittest/API/APITest.py
@@ -398,6 +398,12 @@ def makeConstantRoundTrip(self, file_ending):
         ms["pybool"][SCALAR].reset_dataset(DS(DT.BOOL, extent))
         ms["pybool"][SCALAR].make_constant(False)
 
+        # just testing the data_order attribute
+        ms["char"].data_order = 'C'
+        ms["pyint"].data_order = 'F'
+        self.assertEqual(ms["char"].data_order, 'C')
+        self.assertEqual(ms["pyint"].data_order, 'F')
+
         # staggering meta data
         ms["pyint"][SCALAR].position = [0.25, 0.5]
         ms["pyfloat"][SCALAR].position = [0.5, 0.75]
@@ -457,6 +463,9 @@ def makeConstantRoundTrip(self, file_ending):
         o = [1, 2, 3]
         e = [1, 1, 1]
 
+        self.assertEqual(ms["char"].data_order, 'C')
+        self.assertEqual(ms["pyint"].data_order, 'F')
+
         self.assertTrue(ms["char"].scalar)
         self.assertTrue(ms["pyint"].scalar)
         self.assertTrue(ms["pyfloat"].scalar)
@@ -1267,11 +1276,6 @@ def testParticles(self):
             # self.assertIsInstance(ps, str)
             self.assertIsInstance(i.particles[ps], io.ParticleSpecies)
 
-    def testData_Order(self):
-        """ Test Data_Order. """
-        obj = io.Data_Order('C')
-        del obj
-
     def testDatatype(self):
         """ Test Datatype. """
         data_type = io.Datatype(1)
@@ -1320,6 +1324,7 @@ def testMesh(self):
         self.assertRaises(TypeError, io.Mesh)
         mesh = self.__series.iterations[100].meshes['E']
         copy_mesh = io.Mesh(mesh)
+        self.assertEqual(mesh.data_order, 'C')
 
         self.assertIsInstance(copy_mesh, io.Mesh)