diff --git a/include/openPMD/CustomHierarchy.hpp b/include/openPMD/CustomHierarchy.hpp
index 2aadadccf2..a83af8f614 100644
--- a/include/openPMD/CustomHierarchy.hpp
+++ b/include/openPMD/CustomHierarchy.hpp
@@ -39,8 +39,17 @@ namespace internal
 {
     struct MeshesParticlesPath
     {
-        std::set<std::string> paths;
+        std::optional<std::string> meshesPath;
+        std::optional<std::string> particlesPath;
         [[nodiscard]] bool ignore(std::string const &name) const;
+        [[nodiscard]] inline bool hasMeshes() const noexcept
+        {
+            return meshesPath.has_value();
+        }
+        [[nodiscard]] inline bool hasParticles() const noexcept
+        {
+            return particlesPath.has_value();
+        }
     };
 
     struct CustomHierarchyData : ContainerData<CustomHierarchy>
@@ -74,6 +83,9 @@ class CustomHierarchy : public Container<CustomHierarchy>
         return *m_customHierarchyData;
     }
 
+    void readMeshes(std::string const &meshesPath);
+    void readParticles(std::string const &particlesPath);
+
 protected:
     CustomHierarchy();
     CustomHierarchy(NoInit);
diff --git a/include/openPMD/Iteration.hpp b/include/openPMD/Iteration.hpp
index eeb9d791e5..73b596b88b 100644
--- a/include/openPMD/Iteration.hpp
+++ b/include/openPMD/Iteration.hpp
@@ -305,8 +305,6 @@ class Iteration : public CustomHierarchy
         std::string filePath, std::string const &groupPath, bool beginStep);
     void readGorVBased(std::string const &groupPath, bool beginStep);
     void read_impl(std::string const &groupPath);
-    void readMeshes(std::string const &meshesPath);
-    void readParticles(std::string const &particlesPath);
 
     /**
      * Status after beginning an IO step. Currently includes:
diff --git a/include/openPMD/Mesh.hpp b/include/openPMD/Mesh.hpp
index 17ce9373de..38db4df2e9 100644
--- a/include/openPMD/Mesh.hpp
+++ b/include/openPMD/Mesh.hpp
@@ -41,6 +41,7 @@ class Mesh : public BaseRecord<MeshRecordComponent>
 {
     friend class Container<Mesh>;
     friend class Iteration;
+    friend class CustomHierarchy;
 
 public:
     Mesh(Mesh const &) = default;
diff --git a/include/openPMD/backend/Attributable.hpp b/include/openPMD/backend/Attributable.hpp
index ba64ad1b60..214e7f538a 100644
--- a/include/openPMD/backend/Attributable.hpp
+++ b/include/openPMD/backend/Attributable.hpp
@@ -229,6 +229,7 @@ class Attributable
 
         /** Reconstructs a path that can be passed to a Series constructor */
         std::string filePath() const;
+        std::string printGroup() const;
     };
 
     /**
diff --git a/src/CustomHierarchy.cpp b/src/CustomHierarchy.cpp
index 3860fbf832..9b8571e879 100644
--- a/src/CustomHierarchy.cpp
+++ b/src/CustomHierarchy.cpp
@@ -25,8 +25,10 @@
 #include "openPMD/IO/IOTask.hpp"
 #include "openPMD/RecordComponent.hpp"
 #include "openPMD/Series.hpp"
+#include "openPMD/auxiliary/StringManip.hpp"
 #include "openPMD/backend/Attributable.hpp"
 
+#include <algorithm>
 #include <deque>
 #include <memory>
 
@@ -36,7 +38,13 @@ namespace internal
 {
     bool MeshesParticlesPath::ignore(const std::string &name) const
     {
-        return paths.find(name) != paths.end();
+        auto no_slashes = [](std::string const &str) {
+            return auxiliary::trim(str, [](char const &c) { return c == '/'; });
+        };
+        return (meshesPath.has_value() &&
+                name == no_slashes(meshesPath.value())) ||
+            (particlesPath.has_value() &&
+             name == no_slashes(particlesPath.value()));
     }
 
     CustomHierarchyData::CustomHierarchyData()
@@ -62,6 +70,126 @@ CustomHierarchy::CustomHierarchy()
 CustomHierarchy::CustomHierarchy(NoInit) : Container_t(NoInit())
 {}
 
+void CustomHierarchy::readMeshes(std::string const &meshesPath)
+{
+    Parameter<Operation::OPEN_PATH> pOpen;
+    Parameter<Operation::LIST_PATHS> pList;
+
+    pOpen.path = meshesPath;
+    IOHandler()->enqueue(IOTask(&meshes, pOpen));
+
+    meshes.readAttributes(ReadMode::FullyReread);
+
+    internal::EraseStaleEntries<decltype(meshes)> map{meshes};
+
+    /* obtain all non-scalar meshes */
+    IOHandler()->enqueue(IOTask(&meshes, pList));
+    IOHandler()->flush(internal::defaultFlushParams);
+
+    Parameter<Operation::LIST_ATTS> aList;
+    for (auto const &mesh_name : *pList.paths)
+    {
+        Mesh &m = map[mesh_name];
+        pOpen.path = mesh_name;
+        aList.attributes->clear();
+        IOHandler()->enqueue(IOTask(&m, pOpen));
+        IOHandler()->enqueue(IOTask(&m, aList));
+        IOHandler()->flush(internal::defaultFlushParams);
+
+        auto att_begin = aList.attributes->begin();
+        auto att_end = aList.attributes->end();
+        auto value = std::find(att_begin, att_end, "value");
+        auto shape = std::find(att_begin, att_end, "shape");
+        if (value != att_end && shape != att_end)
+        {
+            MeshRecordComponent &mrc = m;
+            IOHandler()->enqueue(IOTask(&mrc, pOpen));
+            IOHandler()->flush(internal::defaultFlushParams);
+            mrc.get().m_isConstant = true;
+        }
+        m.read();
+        try
+        {
+            m.read();
+        }
+        catch (error::ReadError const &err)
+        {
+            std::cerr << "Cannot read mesh with name '" << mesh_name
+                      << "' and will skip it due to read error:\n"
+                      << err.what() << std::endl;
+            map.forget(mesh_name);
+        }
+    }
+
+    /* obtain all scalar meshes */
+    Parameter<Operation::LIST_DATASETS> dList;
+    IOHandler()->enqueue(IOTask(&meshes, dList));
+    IOHandler()->flush(internal::defaultFlushParams);
+
+    Parameter<Operation::OPEN_DATASET> dOpen;
+    for (auto const &mesh_name : *dList.datasets)
+    {
+        Mesh &m = map[mesh_name];
+        dOpen.name = mesh_name;
+        IOHandler()->enqueue(IOTask(&m, dOpen));
+        IOHandler()->flush(internal::defaultFlushParams);
+        MeshRecordComponent &mrc = m;
+        IOHandler()->enqueue(IOTask(&mrc, dOpen));
+        IOHandler()->flush(internal::defaultFlushParams);
+        mrc.written() = false;
+        mrc.resetDataset(Dataset(*dOpen.dtype, *dOpen.extent));
+        mrc.written() = true;
+        try
+        {
+            m.read();
+        }
+        catch (error::ReadError const &err)
+        {
+            std::cerr << "Cannot read mesh with name '" << mesh_name
+                      << "' and will skip it due to read error:\n"
+                      << err.what() << std::endl;
+            map.forget(mesh_name);
+        }
+    }
+}
+
+void CustomHierarchy::readParticles(std::string const &particlesPath)
+{
+    Parameter<Operation::OPEN_PATH> pOpen;
+    Parameter<Operation::LIST_PATHS> pList;
+
+    pOpen.path = particlesPath;
+    IOHandler()->enqueue(IOTask(&particles, pOpen));
+
+    particles.readAttributes(ReadMode::FullyReread);
+
+    /* obtain all particle species */
+    pList.paths->clear();
+    IOHandler()->enqueue(IOTask(&particles, pList));
+    IOHandler()->flush(internal::defaultFlushParams);
+
+    internal::EraseStaleEntries<decltype(particles)> map{particles};
+    for (auto const &species_name : *pList.paths)
+    {
+        ParticleSpecies &p = map[species_name];
+        pOpen.path = species_name;
+        IOHandler()->enqueue(IOTask(&p, pOpen));
+        IOHandler()->flush(internal::defaultFlushParams);
+        try
+        {
+            p.read();
+        }
+        catch (error::ReadError const &err)
+        {
+            std::cerr << "Cannot read particle species with name '"
+                      << species_name
+                      << "' and will skip it due to read error:\n"
+                      << err.what() << std::endl;
+            map.forget(species_name);
+        }
+    }
+}
+
 void CustomHierarchy::read(internal::MeshesParticlesPath const &mpp)
 {
     /*
@@ -69,9 +197,73 @@ void CustomHierarchy::read(internal::MeshesParticlesPath const &mpp)
      * Path is created/opened already at entry point of method, method needs
      * to create/open path for contained subpaths.
      */
-    Attributable::readAttributes(ReadMode::FullyReread);
+
     Parameter<Operation::LIST_PATHS> pList;
     IOHandler()->enqueue(IOTask(this, pList));
+    IOHandler()->flush(internal::defaultFlushParams);
+
+    auto thisGroupHasMeshesOrParticles =
+        [&pList](std::optional<std::string> meshesOrParticlesPath) -> bool {
+        if (!meshesOrParticlesPath.has_value())
+        {
+            return false;
+        }
+        auto no_slashes = [](std::string const &str) {
+            return auxiliary::trim(str, [](char const &c) { return c == '/'; });
+        };
+        std::string look_for = no_slashes(meshesOrParticlesPath.value());
+        auto const &paths = *pList.paths;
+        return std::find_if(
+                   paths.begin(),
+                   paths.end(),
+                   [&no_slashes, &look_for](std::string const &entry) {
+                       return no_slashes(entry) == look_for;
+                   }) != paths.end();
+    };
+
+    if (thisGroupHasMeshesOrParticles(mpp.meshesPath))
+    {
+        try
+        {
+            readMeshes(mpp.meshesPath.value());
+        }
+        catch (error::ReadError const &err)
+        {
+            std::cerr << "Cannot read meshes at location '"
+                      << myPath().printGroup()
+                      << "' and will skip them due to read error:\n"
+                      << err.what() << std::endl;
+            meshes = {};
+            meshes.dirty() = false;
+        }
+    }
+    else
+    {
+        meshes.dirty() = false;
+    }
+
+    if (thisGroupHasMeshesOrParticles(mpp.particlesPath))
+    {
+        try
+        {
+            readParticles(mpp.particlesPath.value());
+        }
+        catch (error::ReadError const &err)
+        {
+            std::cerr << "Cannot read particles at location '"
+                      << myPath().printGroup()
+                      << "' and will skip them due to read error:\n"
+                      << err.what() << std::endl;
+            particles = {};
+            particles.dirty() = false;
+        }
+    }
+    else
+    {
+        particles.dirty() = false;
+    }
+
+    Attributable::readAttributes(ReadMode::FullyReread);
     Parameter<Operation::LIST_DATASETS> dList;
     IOHandler()->enqueue(IOTask(this, dList));
     IOHandler()->flush(internal::defaultFlushParams);
@@ -147,7 +339,7 @@ void CustomHierarchy::flush(
          * meshesPath and particlesPath are stored there */
         Series s = retrieveSeries();
 
-        if (!meshes.empty() || s.containsAttribute("meshesPath"))
+        if (!meshes.empty())
         {
             if (!s.containsAttribute("meshesPath"))
             {
@@ -163,7 +355,7 @@ void CustomHierarchy::flush(
             meshes.dirty() = false;
         }
 
-        if (!particles.empty() || s.containsAttribute("particlesPath"))
+        if (!particles.empty())
         {
             if (!s.containsAttribute("particlesPath"))
             {
diff --git a/src/Iteration.cpp b/src/Iteration.cpp
index 5a595338da..a3c7f55d16 100644
--- a/src/Iteration.cpp
+++ b/src/Iteration.cpp
@@ -43,6 +43,8 @@ Iteration::Iteration() : CustomHierarchy(NoInit())
     setTime(static_cast<double>(0));
     setDt(static_cast<double>(1));
     setTimeUnitSI(1);
+    meshes.writable().ownKeyWithinParent = "meshes";
+    particles.writable().ownKeyWithinParent = "particles";
 }
 
 template <typename T>
@@ -418,17 +420,15 @@ void Iteration::read_impl(std::string const &groupPath)
             "found " +
                 datatypeToString(Attribute(*aRead.resource).dtype) + ")");
 
-    /* Find the root point [Series] of this file,
-     * meshesPath and particlesPath are stored there */
     Series s = retrieveSeries();
 
     Parameter<Operation::LIST_PATHS> pList;
+    IOHandler()->enqueue(IOTask(this, pList));
     std::string version = s.openPMD();
     bool hasMeshes = false;
     bool hasParticles = false;
     if (version == "1.0.0" || version == "1.0.1")
     {
-        IOHandler()->enqueue(IOTask(this, pList));
         IOHandler()->flush(internal::defaultFlushParams);
         hasMeshes = std::count(
                         pList.paths->begin(),
@@ -447,62 +447,17 @@ void Iteration::read_impl(std::string const &groupPath)
         hasParticles = s.containsAttribute("particlesPath");
     }
 
-    {
-        internal::MeshesParticlesPath mpp;
-        if (hasMeshes)
-        {
-            mpp.paths.emplace(auxiliary::trim(
-                s.meshesPath(), [](char const &c) { return c == '/'; }));
-        }
-        if (hasParticles)
-        {
-            mpp.paths.emplace(auxiliary::trim(
-                s.particlesPath(), [](char const &c) { return c == '/'; }));
-        }
-        CustomHierarchy::read(mpp);
-    }
-
+    internal::MeshesParticlesPath mpp;
     if (hasMeshes)
     {
-        try
-        {
-            readMeshes(s.meshesPath());
-        }
-        catch (error::ReadError const &err)
-        {
-            std::cerr << "Cannot read meshes in iteration " << groupPath
-                      << " and will skip them due to read error:\n"
-                      << err.what() << std::endl;
-            meshes = {};
-            meshes.dirty() = false;
-        }
-    }
-    else
-    {
-        meshes.dirty() = false;
+        mpp.meshesPath = s.meshesPath();
     }
-
     if (hasParticles)
     {
-        try
-        {
-            readParticles(s.particlesPath());
-        }
-        catch (error::ReadError const &err)
-        {
-            std::cerr << "Cannot read particles in iteration " << groupPath
-                      << " and will skip them due to read error:\n"
-                      << err.what() << std::endl;
-            particles = {};
-            particles.dirty() = false;
-        }
-    }
-    else
-    {
-        particles.dirty() = false;
+        mpp.particlesPath = s.particlesPath();
     }
+    CustomHierarchy::read(mpp);
 
-    readAttributes(ReadMode::FullyReread);
 #ifdef openPMD_USE_INVASIVE_TESTS
     if (containsAttribute("__openPMD_internal_fail"))
     {
@@ -515,126 +470,6 @@ void Iteration::read_impl(std::string const &groupPath)
 #endif
 }
 
-void Iteration::readMeshes(std::string const &meshesPath)
-{
-    Parameter<Operation::OPEN_PATH> pOpen;
-    Parameter<Operation::LIST_PATHS> pList;
-
-    pOpen.path = meshesPath;
-    IOHandler()->enqueue(IOTask(&meshes, pOpen));
-
-    meshes.readAttributes(ReadMode::FullyReread);
-
-    internal::EraseStaleEntries<decltype(meshes)> map{meshes};
-
-    /* obtain all non-scalar meshes */
-    IOHandler()->enqueue(IOTask(&meshes, pList));
-    IOHandler()->flush(internal::defaultFlushParams);
-
-    Parameter<Operation::LIST_ATTS> aList;
-    for (auto const &mesh_name : *pList.paths)
-    {
-        Mesh &m = map[mesh_name];
-        pOpen.path = mesh_name;
-        aList.attributes->clear();
-        IOHandler()->enqueue(IOTask(&m, pOpen));
-        IOHandler()->enqueue(IOTask(&m, aList));
-        IOHandler()->flush(internal::defaultFlushParams);
-
-        auto att_begin = aList.attributes->begin();
-        auto att_end = aList.attributes->end();
-        auto value = std::find(att_begin, att_end, "value");
-        auto shape = std::find(att_begin, att_end, "shape");
-        if (value != att_end && shape != att_end)
-        {
-            MeshRecordComponent &mrc = m;
-            IOHandler()->enqueue(IOTask(&mrc, pOpen));
-            IOHandler()->flush(internal::defaultFlushParams);
-            mrc.get().m_isConstant = true;
-        }
-        m.read();
-        try
-        {
-            m.read();
-        }
-        catch (error::ReadError const &err)
-        {
-            std::cerr << "Cannot read mesh with name '" << mesh_name
-                      << "' and will skip it due to read error:\n"
-                      << err.what() << std::endl;
-            map.forget(mesh_name);
-        }
-    }
-
-    /* obtain all scalar meshes */
-    Parameter<Operation::LIST_DATASETS> dList;
-    IOHandler()->enqueue(IOTask(&meshes, dList));
-    IOHandler()->flush(internal::defaultFlushParams);
-
-    Parameter<Operation::OPEN_DATASET> dOpen;
-    for (auto const &mesh_name : *dList.datasets)
-    {
-        Mesh &m = map[mesh_name];
-        dOpen.name = mesh_name;
-        IOHandler()->enqueue(IOTask(&m, dOpen));
-        IOHandler()->flush(internal::defaultFlushParams);
-        MeshRecordComponent &mrc = m;
-        IOHandler()->enqueue(IOTask(&mrc, dOpen));
-        IOHandler()->flush(internal::defaultFlushParams);
-        mrc.written() = false;
-        mrc.resetDataset(Dataset(*dOpen.dtype, *dOpen.extent));
-        mrc.written() = true;
-        try
-        {
-            m.read();
-        }
-        catch (error::ReadError const &err)
-        {
-            std::cerr << "Cannot read mesh with name '" << mesh_name
-                      << "' and will skip it due to read error:\n"
-                      << err.what() << std::endl;
-            map.forget(mesh_name);
-        }
-    }
-}
-
-void Iteration::readParticles(std::string const &particlesPath)
-{
-    Parameter<Operation::OPEN_PATH> pOpen;
-    Parameter<Operation::LIST_PATHS> pList;
-
-    pOpen.path = particlesPath;
-    IOHandler()->enqueue(IOTask(&particles, pOpen));
-
-    particles.readAttributes(ReadMode::FullyReread);
-
-    /* obtain all particle species */
-    pList.paths->clear();
-    IOHandler()->enqueue(IOTask(&particles, pList));
-    IOHandler()->flush(internal::defaultFlushParams);
-
-    internal::EraseStaleEntries<decltype(particles)> map{particles};
-    for (auto const &species_name : *pList.paths)
-    {
-        ParticleSpecies &p = map[species_name];
-        pOpen.path = species_name;
-        IOHandler()->enqueue(IOTask(&p, pOpen));
-        IOHandler()->flush(internal::defaultFlushParams);
-        try
-        {
-            p.read();
-        }
-        catch (error::ReadError const &err)
-        {
-            std::cerr << "Cannot read particle species with name '"
-                      << species_name
-                      << "' and will skip it due to read error:\n"
-                      << err.what() << std::endl;
-            map.forget(species_name);
-        }
-    }
-}
-
 auto Iteration::beginStep(bool reread) -> BeginStepStatus
 {
     BeginStepStatus res;
diff --git a/src/backend/Attributable.cpp b/src/backend/Attributable.cpp
index cfa9bdf15f..68d03bfc07 100644
--- a/src/backend/Attributable.cpp
+++ b/src/backend/Attributable.cpp
@@ -28,6 +28,7 @@
 #include <complex>
 #include <iostream>
 #include <set>
+#include <sstream>
 
 namespace openPMD
 {
@@ -187,6 +188,21 @@ std::string Attributable::MyPath::filePath() const
     return directory + seriesName + seriesExtension;
 }
 
+std::string Attributable::MyPath::printGroup() const
+{
+    std::stringstream res;
+    if (this->group.empty())
+    {
+        return res.str();
+    }
+    res << *this->group.begin();
+    for (auto it = ++this->group.begin(); it != this->group.end(); ++it)
+    {
+        res << '/' << *it;
+    }
+    return res.str();
+}
+
 auto Attributable::myPath() const -> MyPath
 {
     MyPath res;
diff --git a/test/CoreTest.cpp b/test/CoreTest.cpp
index d3a95b7f7a..7ac8b82c9e 100644
--- a/test/CoreTest.cpp
+++ b/test/CoreTest.cpp
@@ -242,6 +242,44 @@ TEST_CASE("custom_hierarchies", "[core]")
         REQUIRE(constant_dataset.getDatatype() == Datatype::FLOAT);
         REQUIRE(constant_dataset.getExtent() == Extent{0, 0, 0});
     }
+    read.close();
+
+    write = Series(filePath, Access::READ_WRITE);
+    {
+        std::vector<int> data(10, 3);
+
+        auto E_x = write.iterations[0]["custom_meshes"].meshes["E"]["x"];
+        E_x.resetDataset({Datatype::INT, {10}});
+        E_x.storeChunk(data, {0}, {10});
+
+        auto e_pos_x = write.iterations[0]["custom_particles"]
+                           .particles["e"]["position"]["x"];
+        e_pos_x.resetDataset({Datatype::INT, {10}});
+        e_pos_x.storeChunk(data, {0}, {10});
+        write.close();
+    }
+
+    read = Series(filePath, Access::READ_ONLY);
+    {
+        REQUIRE(read.iterations[0]["custom_meshes"].meshes.size() == 1);
+        REQUIRE(read.iterations[0]["custom_meshes"].meshes.count("E") == 1);
+        auto E_x_loaded = read.iterations[0]["custom_meshes"]
+                              .meshes["E"]["x"]
+                              .loadChunk<int>();
+        REQUIRE(read.iterations[0]["custom_particles"].particles.size() == 1);
+        REQUIRE(
+            read.iterations[0]["custom_particles"].particles.count("e") == 1);
+        auto e_pos_x_loaded = read.iterations[0]["custom_particles"]
+                                  .particles["e"]["position"]["x"]
+                                  .loadChunk<int>();
+        read.flush();
+
+        for (size_t i = 0; i < 10; ++i)
+        {
+            REQUIRE(E_x_loaded.get()[i] == 3);
+            REQUIRE(e_pos_x_loaded.get()[i] == 3);
+        }
+    }
 }
 
 TEST_CASE("myPath", "[core]")