diff --git a/README.md b/README.md
index 0ccdfaaa..d6c1b583 100644
--- a/README.md
+++ b/README.md
@@ -14,7 +14,6 @@ Python package for brain decoding analysis
 - scikit-learn
 - pandas
 - h5py
-- hdf5storage
 - pyyaml
 
 ### Optional requirements
diff --git a/bdpy/dataform/_mat_v73.py b/bdpy/dataform/_mat_v73.py
index 36f8d18c..4f00262b 100644
--- a/bdpy/dataform/_mat_v73.py
+++ b/bdpy/dataform/_mat_v73.py
@@ -1,10 +1,16 @@
-"""Read MATLAB v7.3 (HDF5) ``.mat`` files with h5py.
+"""h5py-based legacy reading support for MATLAB v7.3 / hdf5storage ``.mat`` files.
 
-bdpy historically relied on :func:`hdf5storage.loadmat` for the load path, but
-that broke under NumPy 2.0, which removed ``np.unicode_`` that older hdf5storage
-referenced (see issue #106). This module reimplements the *read* side on top of
-h5py for the data layouts bdpy actually uses (dense numeric arrays and the
-sparse-array struct). Saving still goes through hdf5storage.
+bdpy historically relied on a third-party MATLAB-v7.3 library (hdf5storage) to
+read ``.mat`` files, but it broke under NumPy 2.0 (it referenced the removed
+``np.unicode_``; see issue #106) and added an extra dependency. This module
+replaces that *read* path with h5py, understanding the MATLAB v7.3 / hdf5storage
+on-disk conventions (column-major / reversed dimension order, ``MATLAB_class``
+and ``Python.*`` attributes), so existing files still load correctly.
+
+This module is **read-only**. bdpy no longer writes MATLAB-compatible ``.mat``
+files; new files are saved as bdpy-native plain HDF5 directly at the save sites
+(see ``bdpy/dataform/sparse.py`` and ``bdpy/dataform/features.py``). The only
+compatibility promise is reading existing hdf5storage / MATLAB v7.3 / bdpy files.
 
 This file is a part of BdPy.
 """
@@ -13,16 +19,21 @@
 import numpy as np
 import scipy.io as sio
 
-__all__ = ["load_array", "loadmat_key", "read_cell", "read_dataset"]
+__all__ = [
+    "load_array",
+    "loadmat_key",
+    "read_cell",
+    "read_dataset",
+]
 
 
 def read_dataset(dset: h5py.Dataset) -> np.ndarray:
-    """Read an h5py dataset, undoing MATLAB v7.3 / hdf5storage conventions.
+    """Read an h5py dataset, undoing MATLAB v7.3 conventions.
 
     MATLAB stores arrays in Fortran (column-major) order, so multi-dimensional
-    datasets are written transposed relative to NumPy's C order. hdf5storage
-    additionally records the original Python shape and empty-array flags as
-    ``Python.*`` attributes, which we honor to reproduce ``hdf5storage.loadmat``.
+    datasets are written transposed relative to NumPy's C order. Older bdpy
+    files additionally record the original Python shape and empty-array flags as
+    ``Python.*`` attributes, which we honor to reproduce the original array.
     Only ``Python.Empty`` is special-cased; a bare ``MATLAB_empty`` (set by
     MATLAB without ``Python.Shape``) is read through the normal path so that
     empty non-scalar arrays such as ``(0, 3)`` keep their shape.
@@ -39,8 +50,8 @@ def read_dataset(dset: h5py.Dataset) -> np.ndarray:
     """
     attrs = dset.attrs
     if "Python.Empty" in attrs:
-        # Only Python.Empty (written by hdf5storage) implies a Python.Shape we
-        # can trust; fall back to the stored dataset shape if it is missing. A
+        # Only Python.Empty (written by the legacy writer) implies a Python.Shape
+        # we can trust; fall back to the stored dataset shape if it is missing. A
         # bare MATLAB_empty (written by MATLAB without Python.Shape) must NOT be
         # treated this way -- np.empty(()) would collapse e.g. (0, 3) to 0-d --
         # so it falls through to the normal read/transpose path below.
@@ -57,7 +68,7 @@ def read_dataset(dset: h5py.Dataset) -> np.ndarray:
 def read_cell(f: h5py.File, dset: h5py.Dataset) -> list:
     """Read a MATLAB cell array (or a plain matrix) into a list of arrays.
 
-    hdf5storage stores Python tuples/lists as MATLAB cell arrays, i.e. an object
+    Python tuples/lists are stored as MATLAB cell arrays, i.e. an object
     dataset of HDF5 references to the individual elements. Files written by other
     tools (e.g. MATLAB or Julia) may instead store the same information as a plain
     2-D matrix whose rows are the elements; both layouts are handled here.
@@ -77,7 +88,9 @@ def read_cell(f: h5py.File, dset: h5py.Dataset) -> list:
     data = dset[()]
     if isinstance(data, np.ndarray) and data.dtype == object:
         return [read_dataset(f[ref]) for ref in data.ravel()]
-    arr = np.asarray(data)
+    # Plain-matrix layout: route through read_dataset so MATLAB-style transposed
+    # matrices carrying MATLAB_class are de-transposed before we split rows.
+    arr = np.asarray(read_dataset(dset))
     return [arr[i] for i in range(arr.shape[0])]
 
 
@@ -105,8 +118,8 @@ def loadmat_key(path: str, key: str) -> np.ndarray:
 
     MATLAB v5 (and earlier) files are read with :func:`scipy.io.loadmat`; v7.3
     (HDF5) files, which scipy cannot read, fall back to the h5py reader. This
-    preserves v5 support while avoiding hdf5storage on the load path (which
-    breaks under NumPy 2.0).
+    preserves v5 support while avoiding the legacy MATLAB-v7.3 library on the
+    load path (which breaks under NumPy 2.0).
 
     Parameters
     ----------
diff --git a/bdpy/dataform/datastore.py b/bdpy/dataform/datastore.py
index 4098b157..3031d5e4 100644
--- a/bdpy/dataform/datastore.py
+++ b/bdpy/dataform/datastore.py
@@ -251,7 +251,7 @@ def get(self, **kargs):
         return dat
 
     def __load_feature(self, fpath):
-        # v5 .mat via scipy, v7.3 (HDF5) via h5py (avoids hdf5storage under NumPy 2.0).
+        # v5 .mat via scipy, v7.3 (HDF5) via h5py (avoids the legacy MAT-v7.3 library under NumPy 2.0).
         r = _mat_v73.loadmat_key(fpath, self.__variable)
         if self.__squeeze:
             r = np.squeeze(r)
diff --git a/bdpy/dataform/features.py b/bdpy/dataform/features.py
index 193537f0..18bb518c 100644
--- a/bdpy/dataform/features.py
+++ b/bdpy/dataform/features.py
@@ -17,15 +17,15 @@
 from multiprocessing import Pool
 from typing import Any, Dict, List, Optional, Union
 
-import hdf5storage
+import h5py
 import numpy as np
 
 from . import _mat_v73
 
 
 def _load_array_with_key(key: str, path: str) -> np.ndarray:
-    # v5 .mat via scipy, v7.3 (HDF5) via h5py; avoids hdf5storage on the load
-    # path, which breaks under NumPy 2.0 (see bdpy/dataform/_mat_v73.py).
+    # v5 .mat via scipy, v7.3 (HDF5) via h5py; avoids the legacy MAT-v7.3
+    # library on the load path, which breaks under NumPy 2.0 (see _mat_v73.py).
     return _mat_v73.loadmat_key(path, key)
 
 
@@ -535,7 +535,10 @@ def save_feature(feature: np.ndarray, base_dir: str, layer: str, label: str, ver
             print(f'{save_file} already exists. Skipped.')
         return None
 
-    hdf5storage.savemat(save_file, {'feat': feature})
+    # Write a bdpy-native plain HDF5 file (not MATLAB-``load`` compatible);
+    # it is read back via _mat_v73.loadmat_key / read_dataset.
+    with h5py.File(save_file, 'w') as f:
+        f.create_dataset('feat', data=feature)
     if verbose:
         print(f'Saved {save_file}.')
 
diff --git a/bdpy/dataform/sparse.py b/bdpy/dataform/sparse.py
index e22e0b3d..174b717a 100644
--- a/bdpy/dataform/sparse.py
+++ b/bdpy/dataform/sparse.py
@@ -9,17 +9,10 @@
 import os
 
 import h5py
-import hdf5storage
 import numpy as np
 
 from . import _mat_v73
 
-# hdf5storage.savemat can fail for SparseArray.save under NumPy 2.x, especially
-# when overwriting an existing sparse struct. SparseArray.save therefore uses a
-# direct h5py writer only under NumPy >= 2. Dense-array save paths are unchanged
-# and still rely on hdf5storage for MATLAB-v7.3 metadata/compatibility behavior.
-_NUMPY2 = int(np.__version__.split('.')[0]) >= 2
-
 
 def load_array(fname, key='data'):
     """Load an array (dense or sparse)."""
@@ -33,34 +26,35 @@ def load_array(fname, key='data'):
                 return s_ary.dense
             raise RuntimeError('Unsupported group: %s' % key)
         if isinstance(obj, h5py.Dataset):
-            # Dense array (read with h5py; hdf5storage breaks under NumPy 2.0)
+            # Dense array (read with h5py; the legacy MAT-v7.3 library breaks under NumPy 2.0)
             return _mat_v73.read_dataset(obj)
         raise RuntimeError('Unsupported data type: %s' % type(obj))
 
 
 def save_array(fname, array, key='data', dtype=np.float64, sparse=False):
-    """Save an array (dense or sparse)."""
+    """Save an array (dense or sparse).
+
+    Dense arrays are written as bdpy-native plain HDF5 datasets. This is
+    intentionally not MATLAB-``load`` compatible; bdpy reads it back via
+    ``load_array`` / ``_mat_v73``.
+    """
     if sparse:
         # Save as a SparseArray
         s_ary = SparseArray(array.astype(dtype))
         s_ary.save(fname, key=key, dtype=dtype)
     else:
-        # Save as a dense array
-        hdf5storage.savemat(fname,
-                            {key: array.astype(dtype)},
-                            format='7.3', oned_as='column',
-                            store_python_metadata=True)
+        # Save as a dense array (bdpy-native plain HDF5)
+        with h5py.File(fname, 'w') as f:
+            f.create_dataset(key, data=array.astype(dtype))
 
     return None
 
 
 def save_multiarrays(fname, arrays):
-    """Save arrays (dense)."""
-    save_dict = {k: v for k, v in arrays.items()}
-    hdf5storage.savemat(fname,
-                        save_dict,
-                        format='7.3', oned_as='column',
-                        store_python_metadata=True)
+    """Save arrays (dense) as bdpy-native plain HDF5 datasets."""
+    with h5py.File(fname, 'w') as f:
+        for k, v in arrays.items():
+            f.create_dataset(k, data=np.asarray(v))
 
     return None
 
@@ -85,28 +79,12 @@ def dense(self):
         return self.__make_dense()
 
     def save(self, fname, key='data', dtype=np.float64):
-        if _NUMPY2:
-            # Avoid hdf5storage.savemat here (it can fail when overwriting an
-            # existing sparse struct under NumPy 2.x) and write the struct with
-            # h5py instead. We replace only ``key`` in the target file, leaving
-            # any other variables already stored there untouched.
-            self.__save_h5py(fname, key=key, dtype=dtype)
-        else:
-            payload = {key: {u'__bdpy_sparse_arrray': True,
-                             u'index': self.__index,
-                             u'value': self.__value.astype(dtype),
-                             u'shape': self.__shape,
-                             u'background': self.__background}}
-            hdf5storage.savemat(fname, payload, format='7.3',
-                                oned_as='column', store_python_metadata=True)
-        return None
-
-    def __save_h5py(self, fname, key='data', dtype=np.float64):
-        # NumPy-2-only writer for the sparse-array struct. The layout mirrors
-        # what ``__load`` expects: ``index``/``shape`` as plain matrices (which
+        # Write the sparse-array struct as bdpy-native plain HDF5 (not
+        # MATLAB-``load`` compatible). The layout mirrors what ``__load``
+        # expects: ``index``/``shape`` as plain matrices (which
         # _mat_v73.read_cell restores row-by-row) and ``value``/``background``
-        # as plain datasets. Opening in append mode and deleting only ``key``
-        # avoids the unconditional full-file rewrite of the previous fallback.
+        # as plain datasets. Open in append mode when the file exists and
+        # replace only ``key`` so any other top-level variables are preserved.
         index = np.vstack([np.asarray(i, dtype=np.int64).ravel()
                            for i in self.__index])
         mode = 'a' if os.path.exists(fname) else 'w'
@@ -133,7 +111,7 @@ def __make_dense(self):
         return dense
 
     def __load(self, fname, key='data'):
-        # Read with h5py instead of hdf5storage, which breaks under NumPy 2.0.
+        # Read with h5py instead of the legacy MAT-v7.3 library (NumPy 2.0).
         # The struct stores ``index``/``shape`` as cell arrays (object refs) when
         # written by bdpy, or as plain matrices when written by other tools.
         with h5py.File(fname, 'r') as f:
diff --git a/mypy.ini b/mypy.ini
index bdec4597..a8f9d9b9 100644
--- a/mypy.ini
+++ b/mypy.ini
@@ -18,9 +18,6 @@ ignore_missing_imports = True
 [mypy-h5py]
 ignore_missing_imports = True
 
-[mypy-hdf5storage]
-ignore_missing_imports = True
-
 [mypy-nibabel]
 ignore_missing_imports = True
 
diff --git a/pyproject.toml b/pyproject.toml
index d6a27a36..6833c81b 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -20,8 +20,6 @@ dependencies = [
     "scipy",
     "scikit-learn",
     "h5py>=3.0",
-    "hdf5storage; python_version < '3.10'",
-    "hdf5storage>=0.2.0; python_version >= '3.10'",
     "pyyaml",
     "pandas",
     "tqdm",
diff --git a/tests/dataform/test_features.py b/tests/dataform/test_features.py
index d5d3a104..0cebfef9 100644
--- a/tests/dataform/test_features.py
+++ b/tests/dataform/test_features.py
@@ -5,11 +5,12 @@
 import os
 import tempfile
 
+import h5py
 import numpy as np
 from numpy.testing import assert_array_equal
-import hdf5storage
 
-from bdpy.dataform.features import Features
+from bdpy.dataform import _mat_v73
+from bdpy.dataform.features import Features, save_feature
 
 
 def _prepare_mock_data(
@@ -20,9 +21,9 @@ def _prepare_mock_data(
     ) -> dict:
     """Prepare mock data for testing.
 
-    Files are written as MATLAB v7.3 (HDF5) so the loader's h5py path (used for
-    NumPy 2.0 compatibility) is exercised. The stacked arrays are returned so
-    tests can compare against them without re-reading the files.
+    Files are written as bdpy-native plain HDF5 (the same format ``save_feature``
+    now produces) so the loader's h5py path is exercised. The stacked arrays are
+    returned so tests can compare against them without re-reading the files.
     """
     stacked = {}
     for layer_name, shape in zip(mock_layer_names, mock_shapes):
@@ -32,11 +33,10 @@ def _prepare_mock_data(
         # sorts the feature files when collecting labels.
         for image_name in sorted(mock_image_names):
             data = np.random.rand(*shape)
-            hdf5storage.savemat(
-                os.path.join(tmpdir, layer_name, image_name + '.mat'),
-                {'feat': data},
-                format='7.3',
-                store_python_metadata=True)
+            with h5py.File(
+                    os.path.join(tmpdir, layer_name, image_name + '.mat'),
+                    'w') as f:
+                f.create_dataset('feat', data=data)
             arrays.append(data)
         stacked[layer_name] = np.vstack(arrays)
     return stacked
@@ -120,5 +120,31 @@ def test_features_get_label(self):
             self.alexnet_conv5_all[index, :]
         )
 
+
+class TestSaveFeature(unittest.TestCase):
+    def test_save_feature_writes_plain_hdf5_readable_by_loader(self):
+        feature = np.random.rand(1, 100)
+        with tempfile.TemporaryDirectory() as tmpdir:
+            save_feature(feature, tmpdir, 'fc8', 'image0001')
+            save_file = os.path.join(tmpdir, 'fc8', 'image0001.mat')
+            self.assertTrue(os.path.exists(save_file))
+
+            # Read back through the loader bdpy uses for feature files.
+            out = _mat_v73.loadmat_key(save_file, 'feat')
+            assert_array_equal(out, feature)
+
+            # The file is bdpy-native plain HDF5: no MATLAB metadata.
+            with h5py.File(save_file, 'r') as f:
+                self.assertNotIn('MATLAB_class', f['feat'].attrs)
+                self.assertNotIn('Python.Shape', f['feat'].attrs)
+
+    def test_save_feature_readable_via_features(self):
+        feature = np.random.rand(1, 100)
+        with tempfile.TemporaryDirectory() as tmpdir:
+            save_feature(feature, tmpdir, 'fc8', 'image0001')
+            feat = Features(tmpdir)
+            assert_array_equal(feat.get_features('fc8'), feature)
+
+
 if __name__ == "__main__":
     unittest.main()
diff --git a/tests/dataform/test_mat_v73.py b/tests/dataform/test_mat_v73.py
index a58f37d6..f2bd82d2 100644
--- a/tests/dataform/test_mat_v73.py
+++ b/tests/dataform/test_mat_v73.py
@@ -55,5 +55,44 @@ def test_matlab_empty_without_matlab_class_keeps_on_disk_shape(self):
             self.assertEqual(out.size, 0)
 
 
+class TestReadCell(unittest.TestCase):
+
+    def test_plain_matrix_with_matlab_class_is_detransposed(self):
+        # A plain (non-object) matrix written MATLAB-style is stored transposed
+        # and tagged with MATLAB_class. read_cell must route it through
+        # read_dataset so it is de-transposed before being split into rows.
+        original = np.array([[1, 2, 3], [4, 5, 6]], dtype=np.float64)  # (2, 3)
+        with tempfile.TemporaryDirectory() as tmpdir:
+            fname = os.path.join(tmpdir, 'cell_matrix.mat')
+            with h5py.File(fname, 'w') as f:
+                # On-disk column-major layout: store the transpose (3, 2).
+                dset = f.create_dataset('a', data=np.ascontiguousarray(original.T))
+                dset.attrs['MATLAB_class'] = np.bytes_(b'double')
+
+            with h5py.File(fname, 'r') as f:
+                rows = _mat_v73.read_cell(f, f['a'])
+
+        # Expect the original (2, 3) orientation split into 2 rows.
+        self.assertEqual(len(rows), 2)
+        np.testing.assert_array_equal(rows[0], original[0])
+        np.testing.assert_array_equal(rows[1], original[1])
+
+    def test_plain_matrix_without_matlab_class_splits_rows_as_stored(self):
+        # No MATLAB_class -> no transpose; rows are returned as stored. This is
+        # how bdpy's own (plain) index/shape datasets are read.
+        stored = np.array([[0, 1, 2], [3, 4, 5]], dtype=np.int64)
+        with tempfile.TemporaryDirectory() as tmpdir:
+            fname = os.path.join(tmpdir, 'plain_matrix.mat')
+            with h5py.File(fname, 'w') as f:
+                f.create_dataset('a', data=stored)
+
+            with h5py.File(fname, 'r') as f:
+                rows = _mat_v73.read_cell(f, f['a'])
+
+        self.assertEqual(len(rows), 2)
+        np.testing.assert_array_equal(rows[0], stored[0])
+        np.testing.assert_array_equal(rows[1], stored[1])
+
+
 if __name__ == '__main__':
     unittest.main()
diff --git a/tests/dataform/test_sparse.py b/tests/dataform/test_sparse.py
index 3047f6ba..9ab7db01 100644
--- a/tests/dataform/test_sparse.py
+++ b/tests/dataform/test_sparse.py
@@ -7,7 +7,7 @@
 import h5py
 import numpy as np
 
-from bdpy.dataform.sparse import load_array, save_array
+from bdpy.dataform.sparse import load_array, save_array, save_multiarrays
 
 
 class TestSparse(unittest.TestCase):
@@ -27,6 +27,32 @@ def test_load_save_dense_array(self):
 
                 np.testing.assert_array_equal(original_data, from_file)
 
+    def test_dense_save_is_plain_hdf5(self):
+        # New dense saves are bdpy-native plain HDF5: no MATLAB metadata.
+        with tempfile.TemporaryDirectory() as tmpdir:
+            fname = os.path.join(tmpdir, 'dense_plain.mat')
+            save_array(fname, np.random.rand(3, 4), key='testdata')
+            with h5py.File(fname, 'r') as f:
+                self.assertIsInstance(f['testdata'], h5py.Dataset)
+                self.assertNotIn('MATLAB_class', f['testdata'].attrs)
+                self.assertNotIn('Python.Shape', f['testdata'].attrs)
+
+    def test_save_multiarrays(self):
+        arrays = {
+            'a': np.random.rand(2, 3),
+            'b': np.arange(5),
+            'c': np.random.rand(4, 1, 2),
+        }
+        with tempfile.TemporaryDirectory() as tmpdir:
+            fname = os.path.join(tmpdir, 'multi.mat')
+            save_multiarrays(fname, arrays)
+
+            # Each array is a top-level plain HDF5 dataset, reloadable individually.
+            with h5py.File(fname, 'r') as f:
+                self.assertEqual(set(f.keys()), set(arrays.keys()))
+            for key, expected in arrays.items():
+                np.testing.assert_array_equal(load_array(fname, key=key), expected)
+
     def test_load_save_sparse_array(self):
         payloads = [
             [(10,), 'test_array_sparse_ndim1.mat'],  # ndim = 1
@@ -44,10 +70,7 @@ def test_load_save_sparse_array(self):
 
                 np.testing.assert_array_equal(original_data, from_file)
 
-    def test_sparse_save_preserves_other_variables_under_numpy2(self):
-        if int(np.__version__.split('.')[0]) < 2:
-            self.skipTest('NumPy-2-only h5py writer')
-
+    def test_sparse_save_preserves_other_variables(self):
         with tempfile.TemporaryDirectory() as tmpdir:
             fname = os.path.join(tmpdir, 'test_sparse_preserve.mat')
 
diff --git a/tests/env/py27/Pipfile b/tests/env/py27/Pipfile
index 0ecd9be9..e174f0b8 100644
--- a/tests/env/py27/Pipfile
+++ b/tests/env/py27/Pipfile
@@ -10,7 +10,6 @@ numpy = "*"
 scipy = "*"
 scikit-learn = "*"
 h5py = "*"
-hdf5storage = "*"
 pyyaml = "*"
 
 [requires]
diff --git a/tests/env/py38/Pipfile b/tests/env/py38/Pipfile
index 026795c5..758bbec5 100644
--- a/tests/env/py38/Pipfile
+++ b/tests/env/py38/Pipfile
@@ -10,7 +10,6 @@ numpy = "*"
 scipy = "*"
 scikit-learn = "*"
 h5py = "*"
-hdf5storage = "*"
 pyyaml = "*"
 
 [requires]