diff --git a/CMakeLists.txt b/CMakeLists.txt
index 6990d1bcb..1eabdfda1 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -15,9 +15,9 @@ endif()
 project(python-blosc2)
 
 # Update to the latest minimum C-Blosc2 version, and to the actual bundled one
-set(BLOSC2_MIN_VERSION 3.5.0)
-set(BLOSC2_BUNDLED_VERSION v3.5.0)
-#set(BLOSC2_BUNDLED_VERSION a54e259aeb1c0e62d6bbea28934369d4831e4046)
+set(BLOSC2_MIN_VERSION 3.5.1)
+#set(BLOSC2_BUNDLED_VERSION v3.5.0)
+set(BLOSC2_BUNDLED_VERSION e3ff64f261522803c6b9af7be913ef1c530cbc04)
 
 if(WIN32 AND NOT CMAKE_C_COMPILER_ID STREQUAL "Clang")
     message(FATAL_ERROR "Windows builds require clang-cl. Set CC/CXX to clang-cl or configure CMake with -T ClangCL.")
diff --git a/doc/reference/ndarray.rst b/doc/reference/ndarray.rst
index 8c492fce7..01821040e 100644
--- a/doc/reference/ndarray.rst
+++ b/doc/reference/ndarray.rst
@@ -25,6 +25,7 @@ It is a direct alias for ``array.vlmeta`` and uses the same persistent storage.
         __len__
         __getitem__
         __setitem__
+        read_into
         take
 
     Utility Methods
@@ -34,6 +35,7 @@ It is a direct alias for ``array.vlmeta`` and uses the same persistent storage.
     .. automethod:: __len__
     .. automethod:: __getitem__
     .. automethod:: __setitem__
+    .. automethod:: read_into
     .. automethod:: take
 
 Constructors
diff --git a/src/blosc2/blosc2_ext.pyx b/src/blosc2/blosc2_ext.pyx
index f59628236..2910c3455 100644
--- a/src/blosc2/blosc2_ext.pyx
+++ b/src/blosc2/blosc2_ext.pyx
@@ -648,6 +648,9 @@ cdef extern from "b2nd.h":
     int b2nd_get_slice_cbuffer(b2nd_array_t *array,
                                int64_t *start, int64_t *stop,
                                void *buffer, int64_t *buffershape, int64_t buffersize) nogil
+    int b2nd_get_slice_cbuffer_into(b2nd_array_t *array,
+                                    int64_t *start, int64_t *stop,
+                                    void *buffer, int64_t *buffershape, int64_t buffersize) nogil
     int b2nd_set_slice_cbuffer(void *buffer, int64_t *buffershape, int64_t buffersize,
                                int64_t *start, int64_t *stop, b2nd_array_t *array)
     int b2nd_get_slice(b2nd_context_t *ctx, b2nd_array_t **array, b2nd_array_t *src, const int64_t *start,
@@ -4386,6 +4389,73 @@ cdef class NDArray:
         _check_rc(rc, "Error while getting the buffer")
         return arr
 
+    def get_basic_into(self, np.ndarray out, key):
+        """NDArray.read_into() for a basic index: decompress ``self[key]`` into the NumPy array
+        `out` and return it, or return None if *key* is not a basic index.
+
+        The data goes straight into `out` if it has the array's dtype and is a box of a C-order
+        buffer (see _box_layout()); otherwise into a new C-contiguous array, which is returned
+        instead for the caller to copy into `out`.  Keys are as in get_basic_numpy(), but without
+        None; any other key, and any key that is out of range, gives None, so that the caller can
+        use the general indexing path, which also raises the appropriate exception.  A selection
+        of another shape than `out` raises ValueError.
+        """
+        cdef tuple items = key if type(key) is tuple else (key,)
+        cdef Py_ssize_t nkey = len(items), i
+        cdef char[_MAX_BASIC_KEY] kinds
+        cdef int64_t[_MAX_BASIC_KEY] ints
+        cdef Py_ssize_t[_MAX_BASIC_KEY] starts, stops
+        cdef Py_ssize_t nidx = self._classify_basic_key(items, kinds, ints, starts, stops)
+        if nidx < 0:
+            return None
+        for i in range(nkey):
+            if kinds[i] == _KEY_NEWAXIS:
+                return None
+
+        cdef int64_t[B2ND_MAX_DIM] start_, stop_, buffershape_
+        cdef c_bool[B2ND_MAX_DIM] int_dims
+        cdef int64_t tick, buffersize = 0
+        cdef int rc = 0, attempt
+        cdef c_bool direct = False
+        target = None
+        with self._handle_lock:
+            # As in get_basic_numpy(): a key resolved against the cached shape is right if the read
+            # does not find shape changes made by another handle (SWMR); otherwise follow them and
+            # try again.  Also follow them first if the selection does not have out's shape (they
+            # may be why) or is empty (there is no read then).
+            for attempt in range(_BASIC_READ_ATTEMPTS):
+                if attempt > 0:
+                    _check_rc(b2nd_refresh(self.array), "Error while refreshing the array")
+                tick = self.array.last_tick
+                out_shape = self._resolve_basic_key(nkey, kinds, ints, starts, stops, nidx, start_, stop_,
+                                                    int_dims)
+                if out_shape is None:
+                    if attempt > 0:
+                        return None
+                    continue
+                if tuple(out_shape) != (<object> out).shape:
+                    if attempt > 0:
+                        raise ValueError(f"out has shape {(<object> out).shape}, but the selection has "
+                                         f"shape {tuple(out_shape)}")
+                    continue
+                if np.PyArray_SIZE(out) == 0:
+                    if attempt > 0:
+                        return out
+                    continue
+                if target is None:
+                    direct = out.dtype == self.dtype and self._box_layout(out, int_dims, buffershape_,
+                                                                          &buffersize)
+                    target = out if direct else np.empty(out_shape, dtype=self.dtype)
+                if direct:
+                    rc = self._read_slice(start_, stop_, np.PyArray_DATA(out), buffershape_, buffersize, True)
+                else:
+                    rc = self._get_slice_cbuffer(target, start_, stop_)
+                if rc >= 0 and self.array.last_tick == tick:
+                    return target
+        # Out of attempts: as the general path, return what the last read gave
+        _check_rc(rc, "Error while getting the slice")
+        return target
+
     cdef _resolve_basic_key(self, Py_ssize_t nkey, const char *kinds, const int64_t *ints,
                             const Py_ssize_t *starts, const Py_ssize_t *stops, Py_ssize_t nidx,
                             int64_t *start_, int64_t *stop_, c_bool *int_dims=NULL):
@@ -4449,7 +4519,6 @@ cdef class NDArray:
 
     cdef int _get_slice_cbuffer(self, arr, int64_t *start_, int64_t *stop_) except? -9999:
         # The caller holds self._handle_lock; returns the rc of b2nd_get_slice_cbuffer()
-        global _nogil_reads
         cdef int64_t[B2ND_MAX_DIM] buffershape_
         cdef int i
         for i in range(self.array.ndim):
@@ -4459,22 +4528,73 @@ cdef class NDArray:
         cdef int rc = 0
         PyObject_GetBuffer(arr, &view, PyBUF_SIMPLE)
         try:
-            # Decompress without the GIL, so that threads reading other arrays (or other handles
-            # of the same file) run in parallel; everything else that uses this super-chunk
-            # holds the lock meanwhile.
-            if _can_release_gil(self.array.sc):
-                _nogil_reads += 1
-                with nogil:
-                    rc = b2nd_get_slice_cbuffer(self.array, start_, stop_,
-                                                <void *> view.buf, buffershape_, view.len)
-                _nogil_reads -= 1
-            else:
-                rc = b2nd_get_slice_cbuffer(self.array, start_, stop_,
-                                            <void *> view.buf, buffershape_, view.len)
+            rc = self._read_slice(start_, stop_, view.buf, buffershape_, view.len, False)
         finally:
             PyBuffer_Release(&view)
         return rc
 
+    cdef int _read_slice(self, int64_t *start_, int64_t *stop_, void *buf, int64_t *buffershape_,
+                         int64_t buffersize, c_bool into) except? -9999:
+        """b2nd_get_slice_cbuffer(), or b2nd_get_slice_cbuffer_into() if `into`; returns its rc.
+        The caller holds self._handle_lock. The decompression runs without the GIL when no Python
+        callback can run on the way, so that threads reading other arrays (or other handles of the
+        same file) run in parallel; everything else that uses this super-chunk holds the lock
+        meanwhile."""
+        global _nogil_reads
+        cdef int rc
+        if _can_release_gil(self.array.sc):
+            _nogil_reads += 1
+            with nogil:
+                if into:
+                    rc = b2nd_get_slice_cbuffer_into(self.array, start_, stop_, buf, buffershape_, buffersize)
+                else:
+                    rc = b2nd_get_slice_cbuffer(self.array, start_, stop_, buf, buffershape_, buffersize)
+            _nogil_reads -= 1
+        elif into:
+            rc = b2nd_get_slice_cbuffer_into(self.array, start_, stop_, buf, buffershape_, buffersize)
+        else:
+            rc = b2nd_get_slice_cbuffer(self.array, start_, stop_, buf, buffershape_, buffersize)
+        return rc
+
+    cdef c_bool _box_layout(self, np.ndarray out, const c_bool *int_dims, int64_t *buffershape_,
+                            int64_t *buffersize):
+        """Whether `out`, of the shape of a selection whose integer-indexed dimensions are
+        `int_dims`, is a box of a C-order buffer that starts at out's first item; if so, fill that
+        buffer's shape and its size in bytes from there (for b2nd_get_slice_cbuffer_into()).
+        Integer-indexed dimensions and the others of length 1 are size-1 axes of the box, whose
+        strides do not matter. Of the longer ones, the innermost must be contiguous, and every
+        other one's stride a whole multiple of the next one's that spans at least that one's
+        length."""
+        cdef np.npy_intp *shape = np.PyArray_DIMS(out)
+        cdef np.npy_intp *strides = np.PyArray_STRIDES(out)
+        cdef int64_t itemsize = np.PyArray_ITEMSIZE(out)
+        cdef int64_t[B2ND_MAX_DIM] lens, steps
+        cdef int64_t size = itemsize
+        cdef int ndim = self.array.ndim, d, j = 0, inner = -1
+        for d in range(ndim):
+            if int_dims[d]:
+                lens[d] = 1
+                steps[d] = 0
+            else:
+                lens[d] = shape[j]
+                steps[d] = strides[j]
+                j += 1
+            buffershape_[d] = lens[d]
+        for d in range(ndim - 1, -1, -1):  # the long dimensions, innermost first
+            if lens[d] <= 1:
+                continue
+            if inner < 0:
+                if steps[d] != itemsize or itemsize <= 0:
+                    return False
+            elif steps[d] % steps[inner] != 0 or steps[d] // steps[inner] < lens[inner]:
+                return False
+            else:
+                buffershape_[inner] = steps[d] // steps[inner]
+            size += (lens[d] - 1) * steps[d]
+            inner = d
+        buffersize[0] = size
+        return True
+
     @_handle_locked
     def get_1d_span_numpy(self, arr, int64_t nchunk, int32_t start, int32_t nitems):
         if self.ndim != 1:
diff --git a/src/blosc2/ndarray.py b/src/blosc2/ndarray.py
index 2b517307b..23bad401a 100644
--- a/src/blosc2/ndarray.py
+++ b/src/blosc2/ndarray.py
@@ -4967,6 +4967,78 @@ def _remember_read(self, key, nparr: np.ndarray) -> None:
         self._last_read.clear()
         self._last_read[make_key_hashable(key)] = nparr
 
+    def read_into(self, out: np.ndarray, key: Any = None, casting: str = "same_kind") -> np.ndarray:
+        """Decompress ``self[key]`` into the existing NumPy array ``out``.
+
+        This is ``out[...] = self[key]`` without the temporary array: the data is decompressed
+        straight to its place in ``out``, which can be a view into a larger array, e.g. the
+        interior of a padded array, one item of a batch or one channel of a multi-channel array.
+        Nothing of the larger array outside ``out`` is touched. Reusing ``out`` across calls also
+        saves the page faults of a fresh result array every time.
+
+        Parameters
+        ----------
+        out: np.ndarray
+            Writeable NumPy array with the shape of ``self[key]``. Contiguous or not: any layout
+            works. Data goes straight into ``out`` (no temporary) when ``out`` has ``self.dtype`` and
+            its last axis is contiguous, with any strides in the other axes that are whole
+            multiples of the next axis' stride. That covers C-contiguous arrays and every box of one,
+            including views with steps in all but the last axis. Other layouts (e.g. Fortran order,
+            transposed views, steps along the last axis) and other dtypes go through one temporary.
+            For a CPU ``torch.Tensor`` ``t``, pass ``t.numpy()``.
+        key: int, slice, Ellipsis or tuple of them, or any key that :meth:`__getitem__` accepts
+            The selection to read, as in ``self[key]``; None (default) reads the whole array.
+            Integers, slices with step 1 and an Ellipsis decompress straight into ``out``; other
+            keys (steps, integer or boolean arrays, ``None``) go through ``self[key]``.
+        casting: str
+            How ``self.dtype`` may be converted to ``out.dtype``, as in :func:`numpy.copyto`.
+            Default ``"same_kind"``.
+
+        Returns
+        -------
+        out: np.ndarray
+            The ``out`` array.
+
+        Examples
+        --------
+        >>> import numpy as np
+        >>> import blosc2
+        >>> a = blosc2.arange(12, dtype=np.int32, shape=(3, 4))
+        >>> padded = np.full((5, 6), -1, dtype=np.int32)
+        >>> _ = a.read_into(padded[1:4, 1:5])  # into the interior of a padded array
+        >>> padded
+        array([[-1, -1, -1, -1, -1, -1],
+               [-1,  0,  1,  2,  3, -1],
+               [-1,  4,  5,  6,  7, -1],
+               [-1,  8,  9, 10, 11, -1],
+               [-1, -1, -1, -1, -1, -1]], dtype=int32)
+        >>> batch = np.empty((2, 2, 4), dtype=np.int32)
+        >>> _ = a.read_into(batch[1], (slice(1, 3), slice(None)))  # into one item of a batch
+        >>> batch[1]
+        array([[ 4,  5,  6,  7],
+               [ 8,  9, 10, 11]], dtype=int32)
+        """
+        if not isinstance(out, np.ndarray):
+            raise TypeError(f"out must be a NumPy array, not {type(out).__name__}")
+        if not out.flags.writeable:
+            raise ValueError("out is read-only")
+        if key is None:
+            key = ()
+        # Integers, slices with step 1 and an Ellipsis: straight into out, or into a temporary for
+        # another dtype or layout (None: another key, see get_basic_into)
+        value = self.get_basic_into(out, key)
+        if value is out:
+            return out
+        if value is None:  # steps, arrays, newaxis, masks, fields...: through __getitem__
+            value = self[key]
+            if isinstance(value, (blosc2.LazyExpr, NDField)):
+                value = value[:]
+            value = np.asarray(value)
+            if out.shape != value.shape:
+                raise ValueError(f"out has shape {out.shape}, but the selection has shape {value.shape}")
+        np.copyto(out, value, casting=casting)
+        return out
+
     def __setitem__(  # noqa: C901
         self,
         key: int | slice | Sequence[slice | int | np.bool_ | np.ndarray[int | np.bool_] | None] | None,
diff --git a/tests/ndarray/test_read_into.py b/tests/ndarray/test_read_into.py
new file mode 100644
index 000000000..2c7bcaaae
--- /dev/null
+++ b/tests/ndarray/test_read_into.py
@@ -0,0 +1,348 @@
+#######################################################################
+# Copyright (c) 2019-present, Blosc Development Team <blosc@blosc.org>
+# All rights reserved.
+#
+# SPDX-License-Identifier: BSD-3-Clause
+#######################################################################
+
+import threading
+
+import numpy as np
+import pytest
+
+import blosc2
+
+SENTINEL = 77
+
+
+@pytest.fixture(params=["memory", "contiguous", "sparse", "mmap"])
+def storage(request, tmp_path):
+    kind = request.param
+
+    def make(a, chunks, blocks):
+        if kind == "memory":
+            return blosc2.asarray(a, chunks=chunks, blocks=blocks)
+        urlpath = str(tmp_path / ("a.b2nd" if kind != "sparse" else "a.b2frame"))
+        blosc2.asarray(
+            a, chunks=chunks, blocks=blocks, urlpath=urlpath, contiguous=kind != "sparse", mode="w"
+        )
+        return blosc2.open(urlpath, mode="r", **({"mmap_mode": "r"} if kind == "mmap" else {}))
+
+    return make
+
+
+def source(shape, dtype, zero_half=True):
+    a = (np.arange(np.prod(shape)).reshape(shape) * 7 % 1000).astype(dtype)
+    if zero_half:  # some all-zero chunks
+        a[: shape[0] // 2] = 0
+    return a
+
+
+def no_getitem(monkeypatch):
+    """Make NDArray.__getitem__ fail, to check that read_into decompresses straight into out"""
+
+    def fail(self, key):
+        raise AssertionError("read_into went through __getitem__")
+
+    monkeypatch.setattr(blosc2.NDArray, "__getitem__", fail)
+
+
+KEYS = [
+    None,
+    (),
+    Ellipsis,
+    1,
+    -1,
+    (slice(2, 9), slice(1, 6)),
+    (slice(None), 3, slice(-5, None)),
+    (Ellipsis, slice(0, 4)),
+    (2, Ellipsis),
+    (slice(3, 100), slice(-100, 2)),  # clipped like NumPy
+    (np.int64(4), slice(None), np.int32(-2)),
+]
+
+
+def read_and_check(b, a, key, monkeypatch):
+    expected = a[key] if key is not None else a
+    out = np.full(expected.shape, SENTINEL, dtype=a.dtype)
+    no_getitem(monkeypatch)
+    assert b.read_into(out, key) is out
+    np.testing.assert_array_equal(out, expected)
+
+
+@pytest.mark.parametrize("key", KEYS)
+def test_keys(storage, key, monkeypatch):
+    a = source((12, 9, 7), np.float32)
+    read_and_check(storage(a, (5, 4, 4), (2, 2, 3)), a, key, monkeypatch)
+
+
+@pytest.mark.parametrize("dtype", [np.float32, np.int16, np.uint8, np.float64, np.complex64, np.bool_])
+def test_dtypes(storage, dtype, monkeypatch):
+    a = source((12, 9, 7), dtype)
+    read_and_check(storage(a, (5, 4, 4), (2, 2, 3)), a, (slice(1, 11), 4, slice(None)), monkeypatch)
+
+
+@pytest.mark.parametrize(
+    ("key", "pad"),
+    [
+        ((slice(0, 12), slice(0, 9), slice(0, 7)), ((2, 1), (0, 3), (1, 1))),
+        ((slice(5, 12), slice(2, 9), slice(3, 4)), ((0, 0), (4, 4), (2, 0))),
+        ((slice(1, 2), slice(0, 9), slice(0, 7)), ((3, 3), (1, 0), (0, 2))),
+    ],
+)
+def test_interior_of_padded_array(storage, key, pad, monkeypatch):
+    a = source((12, 9, 7), np.float32)
+    b = storage(a, (5, 4, 4), (2, 2, 3))
+    inner = tuple(s.stop - s.start for s in key)
+    padded = np.full(
+        [n + lo + hi for n, (lo, hi) in zip(inner, pad, strict=True)], SENTINEL, dtype=np.float32
+    )
+    view = padded[tuple(slice(lo, lo + n) for n, (lo, _) in zip(inner, pad, strict=True))]
+    expected = np.full_like(padded, SENTINEL)
+    expected[tuple(slice(lo, lo + n) for n, (lo, _) in zip(inner, pad, strict=True))] = a[key]
+    no_getitem(monkeypatch)
+    b.read_into(view, key)
+    np.testing.assert_array_equal(padded, expected)  # nothing outside the view changed
+
+
+def test_batch_items_and_channels(storage, monkeypatch):
+    a = source((4, 10, 11), np.int16)
+    b = storage(a, (1, 6, 5), (1, 3, 5))
+    batch = np.full((3, 2, 6, 11), SENTINEL, dtype=np.int16)
+    no_getitem(monkeypatch)
+    for j in range(3):
+        b.read_into(batch[j], (slice(j, j + 2), slice(2, 8)))  # one item of a batch
+    for j in range(3):
+        np.testing.assert_array_equal(batch[j], a[j : j + 2, 2:8])
+    channels = np.full((6, 4, 10, 11), SENTINEL, dtype=np.int16)
+    b.read_into(channels[3, 1], 2)  # one channel of a multi-channel buffer
+    np.testing.assert_array_equal(channels[3, 1], a[2])
+    assert (np.delete(channels.reshape(24, 10, 11), 3 * 4 + 1, axis=0) == SENTINEL).all()
+
+
+@pytest.mark.parametrize(
+    "make_view",
+    [
+        lambda p: p[::2, :6, :5],  # steps in the outer axes are fine
+        lambda p: p[1:9, ::3, :5],
+        lambda p: p[:8, :6, 2:][:, 1:2, :5],  # size-1 axis
+        lambda p: p[3:4, :, :5][:, ::2],
+    ],
+)
+def test_strided_views_without_temporary(storage, make_view, monkeypatch):
+    a = source((8, 6, 5), np.float64)
+    b = storage(a, (4, 3, 5), (2, 3, 5))
+    parent = np.full((16, 12, 9), SENTINEL, dtype=np.float64)
+    view = make_view(parent)
+    key = tuple(slice(0, n) for n in view.shape)
+    expected = parent.copy()
+    make_view(expected)[...] = a[key]
+    no_getitem(monkeypatch)
+    b.read_into(view, key)
+    np.testing.assert_array_equal(parent, expected)
+
+
+@pytest.mark.parametrize(
+    "make_out",
+    [
+        lambda shape: np.full(shape, SENTINEL, dtype=np.float32, order="F"),
+        lambda shape: np.full(shape[::-1], SENTINEL, dtype=np.float32).T,  # transposed
+        lambda shape: np.full((*shape[:-1], 2 * shape[-1]), SENTINEL, dtype=np.float32)[
+            ..., ::2
+        ],  # last-axis step
+        lambda shape: np.full(shape, SENTINEL, dtype=np.float32)[::-1],  # negative stride
+    ],
+)
+def test_other_layouts(storage, make_out):
+    a = source((9, 8, 7), np.float32)
+    b = storage(a, (4, 4, 4), (2, 2, 2))
+    out = make_out((5, 8, 6))
+    b.read_into(out, (slice(2, 7), slice(None), slice(1, 7)))
+    np.testing.assert_array_equal(out, a[2:7, :, 1:7])
+
+
+def test_other_keys_and_dtypes():
+    a = source((10, 8), np.int16, zero_half=False)
+    b = blosc2.asarray(a, chunks=(4, 4), blocks=(2, 2))
+    out = np.empty((5, 4), dtype=np.int16)
+    b.read_into(out, (slice(0, 10, 2), slice(1, 8, 2)))  # steps
+    np.testing.assert_array_equal(out, a[::2, 1::2])
+    out = np.empty((3, 8), dtype=np.int16)
+    b.read_into(out, [1, 4, 7])  # integer array
+    np.testing.assert_array_equal(out, a[[1, 4, 7]])
+    out = np.empty((2, 8), dtype=np.int16)
+    b.read_into(out, (blosc2.asarray(np.array([1, 3])), slice(None)))  # a blosc2 index array in a tuple
+    np.testing.assert_array_equal(out, a[[1, 3]])
+    out = np.empty(int((a > 10).sum()), dtype=np.int16)
+    # a mask (lazy expression); its values come in chunk order, as from b[b > 10][:]
+    b.read_into(out, b > 10)
+    np.testing.assert_array_equal(out, b[b > 10][:])
+    s = blosc2.asarray(np.array([(1.5, 2), (3.5, 4)], dtype=[("x", "f4"), ("y", "i4")]))
+    out = np.empty(2, dtype=np.float32)
+    s.read_into(out, "x")  # a field
+    np.testing.assert_array_equal(out, [1.5, 3.5])
+    out = np.empty((10, 8), dtype=np.float32)
+    b.read_into(out)  # int16 into float32 (same_kind)
+    np.testing.assert_array_equal(out, a.astype(np.float32))
+    out = np.empty((10, 8), dtype=np.uint8)
+    with pytest.raises(TypeError):
+        b.read_into(out)  # int16 into uint8 is not same_kind
+    b.read_into(out, casting="unsafe")
+    np.testing.assert_array_equal(out, a.astype(np.uint8))
+
+
+def test_empty_and_0d():
+    a = source((6, 5), np.float32)
+    b = blosc2.asarray(a, chunks=(4, 4), blocks=(2, 2))
+    out = np.full((0, 5), SENTINEL, dtype=np.float32)
+    assert b.read_into(out, slice(3, 3)) is out
+    out = np.full((), SENTINEL, dtype=np.float32)
+    b.read_into(out, (4, 2))
+    assert out == a[4, 2]
+    z = blosc2.asarray(np.float32(3.5))
+    out = np.zeros((), dtype=np.float32)
+    z.read_into(out)
+    assert out == 3.5
+    u = blosc2.uninit((), dtype=np.float32)  # reads as zero, as u[()] does
+    out = np.full((), SENTINEL, dtype=np.float32)
+    u.read_into(out)
+    assert out == u[()] == 0
+
+
+def test_errors():
+    a = source((6, 5), np.float32)
+    b = blosc2.asarray(a, chunks=(4, 4), blocks=(2, 2))
+    with pytest.raises(ValueError, match="shape"):
+        b.read_into(np.empty((6, 4), dtype=np.float32))
+    with pytest.raises(ValueError, match="shape"):
+        b.read_into(np.empty((3, 3), dtype=np.float32), (slice(0, 6, 2), slice(0, 4, 2)))
+    out = np.empty((6, 5), dtype=np.float32)
+    out.flags.writeable = False
+    with pytest.raises(ValueError, match="read-only"):
+        b.read_into(out)
+    with pytest.raises(TypeError, match="NumPy array"):
+        b.read_into([[0] * 5] * 6)
+    with pytest.raises(IndexError):
+        b.read_into(np.empty((), dtype=np.float32), (6, 0))
+    with pytest.raises(IndexError):
+        b.read_into(np.empty((), dtype=np.float32), (0, 0, 0))
+
+
+@pytest.mark.parametrize(
+    "key", [None, 15, -1, slice(None), slice(-12, None), (Ellipsis, 3), (np.int64(19), 2), (slice(3, 3), 0)]
+)
+@pytest.mark.parametrize("order", ["C", "F"])
+def test_follows_resize(tmp_path, key, order):
+    """Keys are resolved against the shape another handle just set (SWMR), into C- and F-order out"""
+    urlpath = str(tmp_path / "keys.b2nd")
+    w = blosc2.asarray(
+        np.arange(100).reshape(10, 10), chunks=(4, 4), blocks=(2, 2), urlpath=urlpath, mode="w"
+    )
+    r = blosc2.open(urlpath)
+    _ = r[0]
+    np_key = () if key is None else key
+
+    def check(expected):
+        try:
+            want = expected[np_key]
+        except IndexError:
+            with pytest.raises(IndexError):
+                r.read_into(np.empty((10,), dtype=w.dtype), key)
+        else:
+            out = np.full(want.shape, SENTINEL, dtype=w.dtype, order=order)
+            np.testing.assert_array_equal(r.read_into(out, key), want)
+
+    w.resize((20, 10))
+    w[10:20] = np.arange(100, 200).reshape(10, 10)
+    check(np.arange(200).reshape(20, 10))
+    w.resize((5, 10))
+    check(np.arange(50).reshape(5, 10))
+
+
+def test_stale_out_shape_after_resize(tmp_path):
+    """An out of the shape the selection had before another handle resized the array"""
+    urlpath = str(tmp_path / "stale.b2nd")
+    w = blosc2.asarray(
+        np.arange(100).reshape(10, 10), chunks=(4, 4), blocks=(2, 2), urlpath=urlpath, mode="w"
+    )
+    r = blosc2.open(urlpath)
+    _ = r[0]
+    w.resize((20, 10))
+    with pytest.raises(ValueError, match="shape"):
+        r.read_into(np.empty((10, 10), dtype=w.dtype))
+    out = np.empty((20, 10), dtype=w.dtype)
+    np.testing.assert_array_equal(r.read_into(out)[:10], np.arange(100).reshape(10, 10))
+
+
+@pytest.mark.skipif(blosc2.IS_WASM, reason="Pyodide cannot start threads")
+def test_threads(tmp_path):
+    a = source((64, 32, 32), np.float32, zero_half=False)
+    urlpath = str(tmp_path / "t.b2nd")
+    blosc2.asarray(a, chunks=(16, 16, 16), blocks=(8, 8, 8), urlpath=urlpath, mode="w")
+    out = np.full((8, 66, 32, 32), SENTINEL, dtype=np.float32)
+    errors = []
+
+    def work(t):
+        try:
+            h = blosc2.open(urlpath, mode="r")  # one handle per thread
+            for k in range(t, 64, 8):
+                h.read_into(out[t, 1 + k : 2 + k], (slice(k, k + 1),))
+        except Exception as e:  # pragma: no cover
+            errors.append(e)
+
+    ts = [threading.Thread(target=work, args=(t,)) for t in range(8)]
+    for t in ts:
+        t.start()
+    for t in ts:
+        t.join()
+    assert not errors
+    for t in range(8):
+        for k in range(t, 64, 8):
+            np.testing.assert_array_equal(out[t, 1 + k], a[k])
+
+
+@pytest.mark.parametrize("seed", range(4))
+def test_random_views(seed):
+    """Random keys into random views of larger arrays (steps, size-1 axes, integer axes,
+    reversed axes, Fortran order): read_into equals NumPy and changes nothing outside out"""
+    rng = np.random.default_rng(seed)
+    for it in range(60):
+        ndim = int(rng.integers(1, 5))
+        shape = tuple(int(rng.integers(1, 9)) for _ in range(ndim))
+        chunks = tuple(int(rng.integers(1, s + 1)) for s in shape)
+        blocks = tuple(int(rng.integers(1, c + 1)) for c in chunks)
+        dtype = rng.choice([np.uint8, np.int16, np.float32, np.float64])
+        a = (rng.random(shape) * 100).astype(dtype)  # in range for every dtype (no invalid casts)
+        a[tuple(slice(0, s // 2) for s in shape)] = 0
+        b = blosc2.asarray(a, chunks=chunks, blocks=blocks)
+        key = []
+        for s in shape:
+            if rng.random() < 0.25:
+                key.append(int(rng.integers(-s, s)))
+            else:
+                lo = int(rng.integers(0, s + 1))
+                key.append(slice(lo, int(rng.integers(lo, s + 1))))
+        key = tuple(key)
+        expected = np.asarray(a[key])
+        # a view with the result's shape into a parent with room around it
+        steps = [int(rng.integers(1, 3)) for _ in expected.shape]
+        pads = [int(rng.integers(0, 3)) for _ in expected.shape]
+        pshape = [
+            n * st + p + int(rng.integers(0, 2))
+            for n, st, p in zip(expected.shape, steps, pads, strict=True)
+        ]
+        order = "F" if rng.random() < 0.15 else "C"
+        reverse = expected.ndim > 0 and rng.random() < 0.15
+        sel = tuple(slice(p, p + n * st, st) for n, st, p in zip(expected.shape, steps, pads, strict=True))
+
+        def view_of(arr, sel=sel, reverse=reverse):
+            v = arr[(*sel, Ellipsis)]  # a view also for 0-d arrays
+            return v[::-1] if reverse else v
+
+        parent = np.full(pshape, SENTINEL, dtype=dtype, order=order)
+        ref = parent.copy(order="K")
+        view_of(ref)[...] = expected
+        b.read_into(view_of(parent), key)
+        np.testing.assert_array_equal(
+            parent, ref, err_msg=f"{it} {shape} {chunks} {blocks} {key} {sel} {reverse}"
+        )