diff --git a/include/pybind11/numpy.h b/include/pybind11/numpy.h
index 03e1ed61..ceb06fe3 100644
--- a/include/pybind11/numpy.h
+++ b/include/pybind11/numpy.h
@@ -419,6 +419,10 @@ class unchecked_mutable_reference : public unchecked_reference<T, Dims> {
     using ConstBase::ConstBase;
     using ConstBase::Dynamic;
 public:
+    // Bring in const-qualified versions from base class
+    using ConstBase::operator();
+    using ConstBase::operator[];
+
     /// Mutable, unchecked access to data at the given indices.
     template <typename... Ix> T& operator()(Ix... index) {
         static_assert(ssize_t{sizeof...(Ix)} == Dims || Dynamic,
diff --git a/tests/test_numpy_array.cpp b/tests/test_numpy_array.cpp
index 33f1d785..c6bc5e06 100644
--- a/tests/test_numpy_array.cpp
+++ b/tests/test_numpy_array.cpp
@@ -318,6 +318,18 @@ TEST_SUBMODULE(numpy_array, sm) {
         return auxiliaries(r, r2);
     });
 
+    sm.def("proxy_auxiliaries1_const_ref", [](py::array_t<double> a) {
+        const auto &r = a.unchecked<1>();
+        const auto &r2 = a.mutable_unchecked<1>();
+        return r(0) == r2(0) && r[0] == r2[0];
+    });
+
+    sm.def("proxy_auxiliaries2_const_ref", [](py::array_t<double> a) {
+        const auto &r = a.unchecked<2>();
+        const auto &r2 = a.mutable_unchecked<2>();
+        return r(0, 0) == r2(0, 0);
+    });
+
     // test_array_unchecked_dyn_dims
     // Same as the above, but without a compile-time dimensions specification:
     sm.def("proxy_add2_dyn", [](py::array_t<double> a, double v) {
diff --git a/tests/test_numpy_array.py b/tests/test_numpy_array.py
index a36e707c..e3e7c5f5 100644
--- a/tests/test_numpy_array.py
+++ b/tests/test_numpy_array.py
@@ -364,6 +364,9 @@ def test_array_unchecked_fixed_dims(msg):
     assert m.proxy_auxiliaries2(z1) == [11, 11, True, 2, 8, 2, 2, 4, 32]
     assert m.proxy_auxiliaries2(z1) == m.array_auxiliaries2(z1)
 
+    assert m.proxy_auxiliaries1_const_ref(z1[0, :])
+    assert m.proxy_auxiliaries2_const_ref(z1)
+
 
 def test_array_unchecked_dyn_dims(msg):
     z1 = np.array([[1, 2], [3, 4]], dtype='float64')