diff --git a/.github/workflows/gh-ci.yaml b/.github/workflows/gh-ci.yaml
index 88c6c11a825..c9f84aa0f7c 100644
--- a/.github/workflows/gh-ci.yaml
+++ b/.github/workflows/gh-ci.yaml
@@ -61,7 +61,7 @@ jobs:
               full-deps: false
               install_hole: false
               codecov: false
-              numpy: numpy=1.20.0
+              numpy: numpy=1.21.0
             - name: asv_check
               os: ubuntu-latest
               python-version: 3.8
diff --git a/.github/workflows/linters.yaml b/.github/workflows/linters.yaml
index 835277c84da..256f33b7eb3 100644
--- a/.github/workflows/linters.yaml
+++ b/.github/workflows/linters.yaml
@@ -102,7 +102,7 @@ jobs:
         with:
           mamba: true
           full-deps: true
-          numpy: numpy=1.20.0
+          numpy: numpy=1.21.0
 
       - name: install
         run: |
diff --git a/.travis.yml b/.travis.yml
index e75ecde9232..927f80cd60d 100644
--- a/.travis.yml
+++ b/.travis.yml
@@ -44,7 +44,7 @@ matrix:
       group: edge
       if: type = cron
       before_install:
-        - python -m pip install cython "numpy>=1.20.0" scipy
+        - python -m pip install cython "numpy>=1.21.0" scipy
         - python -m pip install --no-build-isolation hypothesis matplotlib packaging pytest pytest-cov pytest-xdist tqdm threadpoolctl fasteners
       install:
         - cd package
diff --git a/azure-pipelines.yml b/azure-pipelines.yml
index 5233f7a96e4..d8ddc709dc4 100644
--- a/azure-pipelines.yml
+++ b/azure-pipelines.yml
@@ -25,7 +25,7 @@ jobs:
           PYTHON_VERSION: '3.8'
           PYTHON_ARCH: 'x86'
           BUILD_TYPE: 'normal'
-          NUMPY_MIN: '1.20.0'
+          NUMPY_MIN: '1.21.0'
           imageName: 'windows-2019'
         Win-Python38-64bit-full:
           PYTHON_VERSION: '3.8'
@@ -47,7 +47,7 @@ jobs:
           PYTHON_VERSION: '3.8'
           PYTHON_ARCH: 'x64'
           BUILD_TYPE: 'wheel'
-          NUMPY_MIN: '1.20.0'
+          NUMPY_MIN: '1.21.0'
           imageName: 'windows-2019'
         Linux-Python310-64bit-full-wheel:
           PYTHON_VERSION: '3.10'
@@ -59,7 +59,7 @@ jobs:
           PYTHON_VERSION: '3.8'
           PYTHON_ARCH: 'x64'
           BUILD_TYPE: 'wheel'
-          NUMPY_MIN: '1.20.0'
+          NUMPY_MIN: '1.21.0'
           imageName: 'ubuntu-latest'
   pool:
     vmImage: $(imageName)
diff --git a/maintainer/conda/environment.yml b/maintainer/conda/environment.yml
index 3b3e5fcb226..ed3e4ed1d2b 100644
--- a/maintainer/conda/environment.yml
+++ b/maintainer/conda/environment.yml
@@ -16,7 +16,7 @@ dependencies:
   - mmtf-python
   - mock
   - networkx
-  - numpy>=1.20
+  - numpy>=1.21
   - pytest
   - python==3.8
   - pytng>=0.2.3
diff --git a/package/CHANGELOG b/package/CHANGELOG
index fb23cb3823f..dadf1c64459 100644
--- a/package/CHANGELOG
+++ b/package/CHANGELOG
@@ -28,6 +28,9 @@ Enhancements
     and SegmentGroup. (PR #3953)
 
 Changes
+  * As per NEP29 the minimum supported NumPy version has been raised to 1.21
+    (note: in practice later versions of NumPy may be used depending on your
+     architecture, operating system, or Python version) (PR #3983)
   * The deprecated direct indexing and `times` from the `results` attribute of
     analysis.nucleicacids' NucPairDist and WatsonCrickDist classes has been
     removed. Please use the `results.pair_distances` and `times` attributes
diff --git a/package/MDAnalysis/lib/distances.py b/package/MDAnalysis/lib/distances.py
index 785a9712fcc..a36b1e81617 100644
--- a/package/MDAnalysis/lib/distances.py
+++ b/package/MDAnalysis/lib/distances.py
@@ -72,7 +72,7 @@
 .. autofunction:: minimize_vectors(vectors, box)
 """
 import numpy as np
-from numpy.lib.utils import deprecate
+import numpy.typing as npt
 
 from typing import Union, Optional, Callable
 from typing import TYPE_CHECKING
@@ -139,9 +139,8 @@ def _run(funcname: str, args: Optional[tuple] = None,
 from .c_distances_openmp import OPENMP_ENABLED as USED_OPENMP
 
 
-# typing: numpy
-def _check_result_array(result: Optional[np.ndarray],
-                        shape: tuple) -> np.ndarray:
+def _check_result_array(result: Optional[npt.NDArray],
+                        shape: tuple) -> npt.NDArray:
     """Check if the result array is ok to use.
 
     The `result` array must meet the following requirements:
@@ -182,14 +181,13 @@ def _check_result_array(result: Optional[np.ndarray],
     return result
 
 
-# typing: numpy
 @check_coords('reference', 'configuration', reduce_result_if_single=False,
               check_lengths_match=False, allow_atomgroup=True)
-def distance_array(reference: Union[np.ndarray, 'AtomGroup'],
-                   configuration: Union[np.ndarray, 'AtomGroup'],
-                   box: Optional[np.ndarray] = None,
-                   result: Optional[np.ndarray] = None,
-                   backend: str = "serial") -> np.ndarray:
+def distance_array(reference: Union[npt.NDArray, 'AtomGroup'],
+                   configuration: Union[npt.NDArray, 'AtomGroup'],
+                   box: Optional[npt.NDArray] = None,
+                   result: Optional[npt.NDArray] = None,
+                   backend: str = "serial") -> npt.NDArray:
     """Calculate all possible distances between a reference set and another
     configuration.
 
@@ -272,12 +270,11 @@ def distance_array(reference: Union[np.ndarray, 'AtomGroup'],
     return distances
 
 
-# typing: numpy
 @check_coords('reference', reduce_result_if_single=False, allow_atomgroup=True)
-def self_distance_array(reference: Union[np.ndarray, 'AtomGroup'],
-                        box: Optional[np.ndarray] = None,
-                        result: Optional[np.ndarray] = None,
-                        backend: str = "serial") -> np.ndarray:
+def self_distance_array(reference: Union[npt.NDArray, 'AtomGroup'],
+                        box: Optional[npt.NDArray] = None,
+                        result: Optional[npt.NDArray] = None,
+                        backend: str = "serial") -> npt.NDArray:
     """Calculate all possible distances within a configuration `reference`.
 
     If the optional argument `box` is supplied, the minimum image convention is
@@ -359,10 +356,10 @@ def self_distance_array(reference: Union[np.ndarray, 'AtomGroup'],
 @check_coords('reference', 'configuration', enforce_copy=False,
               reduce_result_if_single=False, check_lengths_match=False,
               allow_atomgroup=True)
-def capped_distance(reference: Union[np.ndarray, 'AtomGroup'],
-                    configuration: Union[np.ndarray, 'AtomGroup'],
+def capped_distance(reference: Union[npt.NDArray, 'AtomGroup'],
+                    configuration: Union[npt.NDArray, 'AtomGroup'],
                     max_cutoff: float, min_cutoff: Optional[float] = None,
-                    box: Optional[np.ndarray] = None,
+                    box: Optional[npt.NDArray] = None,
                     method: Optional[str] = None,
                     return_distances: Optional[bool] = True):
     """Calculates pairs of indices corresponding to entries in the `reference`
@@ -455,8 +452,8 @@ def capped_distance(reference: Union[np.ndarray, 'AtomGroup'],
     # The check_coords decorator made sure that reference and configuration
     # are arrays of positions. Mypy does not know about that so we have to
     # tell it.
-    reference_positions: np.ndarray = reference  # type: ignore
-    configuration_positions: np.ndarray = configuration  # type: ignore
+    reference_positions: npt.NDArray = reference  # type: ignore
+    configuration_positions: npt.NDArray = configuration  # type: ignore
     function = _determine_method(reference_positions, configuration_positions,
                                  max_cutoff, min_cutoff=min_cutoff,
                                  box=box, method=method)
@@ -465,9 +462,9 @@ def capped_distance(reference: Union[np.ndarray, 'AtomGroup'],
                     box=box, return_distances=return_distances)
 
 
-def _determine_method(reference: np.ndarray, configuration: np.ndarray,
+def _determine_method(reference: npt.NDArray, configuration: npt.NDArray,
                       max_cutoff: float, min_cutoff: Optional[float] = None,
-                      box: Optional[np.ndarray] = None,
+                      box: Optional[npt.NDArray] = None,
                       method: Optional[str] = None) -> Callable:
     """Guesses the fastest method for capped distance calculations based on the
     size of the coordinate sets and the relative size of the target volume.
@@ -541,10 +538,10 @@ def _determine_method(reference: np.ndarray, configuration: np.ndarray,
 @check_coords('reference', 'configuration', enforce_copy=False,
               reduce_result_if_single=False, check_lengths_match=False,
               allow_atomgroup=True)
-def _bruteforce_capped(reference: Union[np.ndarray, 'AtomGroup'],
-                       configuration: Union[np.ndarray, 'AtomGroup'],
+def _bruteforce_capped(reference: Union[npt.NDArray, 'AtomGroup'],
+                       configuration: Union[npt.NDArray, 'AtomGroup'],
                        max_cutoff: float, min_cutoff: Optional[float] = None,
-                       box: Optional[np.ndarray] = None,
+                       box: Optional[npt.NDArray] = None,
                        return_distances: Optional[bool] = True):
     """Capped distance evaluations using a brute force method.
 
@@ -628,10 +625,10 @@ def _bruteforce_capped(reference: Union[np.ndarray, 'AtomGroup'],
 @check_coords('reference', 'configuration', enforce_copy=False,
               reduce_result_if_single=False, check_lengths_match=False,
               allow_atomgroup=True)
-def _pkdtree_capped(reference: Union[np.ndarray, 'AtomGroup'],
-                    configuration: Union[np.ndarray, 'AtomGroup'],
+def _pkdtree_capped(reference: Union[npt.NDArray, 'AtomGroup'],
+                    configuration: Union[npt.NDArray, 'AtomGroup'],
                     max_cutoff: float, min_cutoff: Optional[float] = None,
-                    box: Optional[np.ndarray] = None,
+                    box: Optional[npt.NDArray] = None,
                     return_distances: Optional[bool] = True):
     """Capped distance evaluations using a KDtree method.
 
@@ -720,10 +717,10 @@ def _pkdtree_capped(reference: Union[np.ndarray, 'AtomGroup'],
 @check_coords('reference', 'configuration', enforce_copy=False,
               reduce_result_if_single=False, check_lengths_match=False,
               allow_atomgroup=True)
-def _nsgrid_capped(reference: Union[np.ndarray, 'AtomGroup'],
-                   configuration: Union[np.ndarray, 'AtomGroup'],
+def _nsgrid_capped(reference: Union[npt.NDArray, 'AtomGroup'],
+                   configuration: Union[npt.NDArray, 'AtomGroup'],
                    max_cutoff: float, min_cutoff: Optional[float] = None,
-                   box: Optional[np.ndarray] = None,
+                   box: Optional[npt.NDArray] = None,
                    return_distances: Optional[bool] = True):
     """Capped distance evaluations using a grid-based search method.
 
@@ -829,10 +826,10 @@ def _nsgrid_capped(reference: Union[np.ndarray, 'AtomGroup'],
 @check_coords('reference', enforce_copy=False,
               reduce_result_if_single=False, check_lengths_match=False,
               allow_atomgroup=True)
-def self_capped_distance(reference: Union[np.ndarray, 'AtomGroup'],
+def self_capped_distance(reference: Union[npt.NDArray, 'AtomGroup'],
                          max_cutoff: float,
                          min_cutoff: Optional[float] = None,
-                         box: Optional[np.ndarray] = None,
+                         box: Optional[npt.NDArray] = None,
                          method: Optional[str] = None,
                          return_distances: Optional[bool] = True):
     """Calculates pairs of indices corresponding to entries in the `reference`
@@ -926,7 +923,7 @@ def self_capped_distance(reference: Union[np.ndarray, 'AtomGroup'],
     # The check_coords decorator made sure that reference is an
     # array of positions. Mypy does not know about that so we have to
     # tell it.
-    reference_positions: np.ndarray = reference  # type: ignore
+    reference_positions: npt.NDArray = reference  # type: ignore
     function = _determine_method_self(reference_positions,
                                       max_cutoff, min_cutoff=min_cutoff,
                                       box=box, method=method)
@@ -934,9 +931,9 @@ def self_capped_distance(reference: Union[np.ndarray, 'AtomGroup'],
                     return_distances=return_distances)
 
 
-def _determine_method_self(reference: np.ndarray, max_cutoff: float,
+def _determine_method_self(reference: npt.NDArray, max_cutoff: float,
                            min_cutoff: Optional[float] = None,
-                           box: Optional[np.ndarray] = None,
+                           box: Optional[npt.NDArray] = None,
                            method: Optional[str] = None):
     """Guesses the fastest method for capped distance calculations based on the
     size of the `reference` coordinate set and the relative size of the target
@@ -1001,10 +998,10 @@ def _determine_method_self(reference: np.ndarray, max_cutoff: float,
 
 @check_coords('reference', enforce_copy=False, reduce_result_if_single=False,
               allow_atomgroup=True)
-def _bruteforce_capped_self(reference: Union[np.ndarray, 'AtomGroup'],
+def _bruteforce_capped_self(reference: Union[npt.NDArray, 'AtomGroup'],
                             max_cutoff: float,
                             min_cutoff: Optional[float] = None,
-                            box: Optional[np.ndarray] = None,
+                            box: Optional[npt.NDArray] = None,
                             return_distances: Optional[bool] = True):
     """Capped distance evaluations using a brute force method.
 
@@ -1084,10 +1081,10 @@ def _bruteforce_capped_self(reference: Union[np.ndarray, 'AtomGroup'],
 
 @check_coords('reference', enforce_copy=False, reduce_result_if_single=False,
               allow_atomgroup=True)
-def _pkdtree_capped_self(reference: Union[np.ndarray, 'AtomGroup'],
+def _pkdtree_capped_self(reference: Union[npt.NDArray, 'AtomGroup'],
                          max_cutoff: float,
                          min_cutoff: Optional[float] = None,
-                         box: Optional[np.ndarray] = None,
+                         box: Optional[npt.NDArray] = None,
                          return_distances: Optional[bool] = True):
     """Capped distance evaluations using a KDtree method.
 
@@ -1168,10 +1165,10 @@ def _pkdtree_capped_self(reference: Union[np.ndarray, 'AtomGroup'],
 
 @check_coords('reference', enforce_copy=False, reduce_result_if_single=False,
               allow_atomgroup=True)
-def _nsgrid_capped_self(reference: Union[np.ndarray, 'AtomGroup'],
+def _nsgrid_capped_self(reference: Union[npt.NDArray, 'AtomGroup'],
                         max_cutoff: float,
                         min_cutoff: Optional[float] = None,
-                        box: Optional[np.ndarray] = None,
+                        box: Optional[npt.NDArray] = None,
                         return_distances: Optional[bool] = True):
     """Capped distance evaluations using a grid-based search method.
 
@@ -1363,13 +1360,12 @@ def transform_StoR(coords, box, backend="serial"):
     return coords
 
 
-# typing: numpy
 @check_coords('coords1', 'coords2', allow_atomgroup=True)
-def calc_bonds(coords1: Union[np.ndarray, 'AtomGroup'],
-               coords2: Union[np.ndarray, 'AtomGroup'],
-               box: Optional[np.ndarray] = None,
-               result: Optional[np.ndarray] = None,
-               backend: str = "serial") -> np.ndarray:
+def calc_bonds(coords1: Union[npt.NDArray, 'AtomGroup'],
+               coords2: Union[npt.NDArray, 'AtomGroup'],
+               box: Optional[npt.NDArray] = None,
+               result: Optional[npt.NDArray] = None,
+               backend: str = "serial") -> npt.NDArray:
     """Calculates the bond lengths between pairs of atom positions from the two
     coordinate arrays `coords1` and `coords2`, which must contain the same
     number of coordinates. ``coords1[i]`` and ``coords2[i]`` represent the
@@ -1455,14 +1451,13 @@ def calc_bonds(coords1: Union[np.ndarray, 'AtomGroup'],
     return bondlengths
 
 
-# typing: numpy
 @check_coords('coords1', 'coords2', 'coords3', allow_atomgroup=True)
-def calc_angles(coords1: Union[np.ndarray, 'AtomGroup'],
-                coords2: Union[np.ndarray, 'AtomGroup'],
-                coords3: Union[np.ndarray, 'AtomGroup'],
-                box: Optional[np.ndarray] = None,
-                result: Optional[np.ndarray] = None,
-                backend: str = "serial") -> np.ndarray:
+def calc_angles(coords1: Union[npt.NDArray, 'AtomGroup'],
+                coords2: Union[npt.NDArray, 'AtomGroup'],
+                coords3: Union[npt.NDArray, 'AtomGroup'],
+                box: Optional[npt.NDArray] = None,
+                result: Optional[npt.NDArray] = None,
+                backend: str = "serial") -> npt.NDArray:
     """Calculates the angles formed between triplets of atom positions from the
     three coordinate arrays `coords1`, `coords2`, and `coords3`. All coordinate
     arrays must contain the same number of coordinates.
@@ -1558,15 +1553,14 @@ def calc_angles(coords1: Union[np.ndarray, 'AtomGroup'],
     return angles
 
 
-# typing: numpy
 @check_coords('coords1', 'coords2', 'coords3', 'coords4', allow_atomgroup=True)
-def calc_dihedrals(coords1: Union[np.ndarray, 'AtomGroup'],
-                   coords2: Union[np.ndarray, 'AtomGroup'],
-                   coords3: Union[np.ndarray, 'AtomGroup'],
-                   coords4: Union[np.ndarray, 'AtomGroup'],
-                   box: Optional[np.ndarray] = None,
-                   result: Optional[np.ndarray] = None,
-                   backend: str = "serial") -> np.ndarray:
+def calc_dihedrals(coords1: Union[npt.NDArray, 'AtomGroup'],
+                   coords2: Union[npt.NDArray, 'AtomGroup'],
+                   coords3: Union[npt.NDArray, 'AtomGroup'],
+                   coords4: Union[npt.NDArray, 'AtomGroup'],
+                   box: Optional[npt.NDArray] = None,
+                   result: Optional[npt.NDArray] = None,
+                   backend: str = "serial") -> npt.NDArray:
     r"""Calculates the dihedral angles formed between quadruplets of positions
     from the four coordinate arrays `coords1`, `coords2`, `coords3`, and
     `coords4`, which must contain the same number of coordinates.
@@ -1677,11 +1671,10 @@ def calc_dihedrals(coords1: Union[np.ndarray, 'AtomGroup'],
     return dihedrals
 
 
-# typing: numpy
 @check_coords('coords', allow_atomgroup=True)
-def apply_PBC(coords: Union[np.ndarray, 'AtomGroup'],
-              box: Optional[np.ndarray] = None,
-              backend: str = "serial") -> np.ndarray:
+def apply_PBC(coords: Union[npt.NDArray, 'AtomGroup'],
+              box: Optional[npt.NDArray] = None,
+              backend: str = "serial") -> npt.NDArray:
     """Moves coordinates into the primary unit cell.
 
     Parameters
@@ -1717,7 +1710,7 @@ def apply_PBC(coords: Union[np.ndarray, 'AtomGroup'],
     """
     # coords is an array, the check_coords decorator made sure of that.
     # Mypy, however, is not aware of that so we have to tell it explicitly.
-    coords_array: np.ndarray = coords  # type: ignore
+    coords_array: npt.NDArray = coords  # type: ignore
 
     if len(coords_array) == 0:
         return coords_array
@@ -1730,9 +1723,8 @@ def apply_PBC(coords: Union[np.ndarray, 'AtomGroup'],
     return coords_array
 
 
-# typing: numpy
 @check_coords('vectors', enforce_copy=False, enforce_dtype=False)
-def minimize_vectors(vectors: np.ndarray, box: np.ndarray) -> np.ndarray:
+def minimize_vectors(vectors: npt.NDArray, box: npt.NDArray) -> npt.NDArray:
     """Apply minimum image convention to an array of vectors
 
     This function is required for calculating the correct vectors between two
diff --git a/package/MDAnalysis/lib/mdamath.py b/package/MDAnalysis/lib/mdamath.py
index fff63d73bd5..260bf0b87d5 100644
--- a/package/MDAnalysis/lib/mdamath.py
+++ b/package/MDAnalysis/lib/mdamath.py
@@ -92,8 +92,7 @@ def norm(v: npt.ArrayLike) -> float:
     return np.sqrt(np.dot(v, v))
 
 
-# typing: numpy
-def normal(vec1: npt.ArrayLike, vec2: npt.ArrayLike) -> np.ndarray:
+def normal(vec1: npt.ArrayLike, vec2: npt.ArrayLike) -> npt.NDArray:
     r"""Returns the unit vector normal to two vectors.
 
     .. math::
@@ -105,7 +104,7 @@ def normal(vec1: npt.ArrayLike, vec2: npt.ArrayLike) -> np.ndarray:
     .. versionchanged:: 0.11.0
        Moved into lib.mdamath
     """
-    normal = np.cross(vec1, vec2)
+    normal: npt.NDArray = np.cross(vec1, vec2)
     n = norm(normal)
     if n == 0.0:
         return normal  # returns [0,0,0] instead of [nan,nan,nan]
@@ -113,8 +112,7 @@ def normal(vec1: npt.ArrayLike, vec2: npt.ArrayLike) -> np.ndarray:
     return normal / n
 
 
-# typing: numpy
-def pdot(a: npt.ArrayLike, b: npt.ArrayLike) -> np.ndarray:
+def pdot(a: npt.NDArray, b: npt.NDArray) -> npt.NDArray:
     """Pairwise dot product.
 
     ``a`` must be the same shape as ``b``.
@@ -131,8 +129,7 @@ def pdot(a: npt.ArrayLike, b: npt.ArrayLike) -> np.ndarray:
     return np.einsum('ij,ij->i', a, b)
 
 
-# typing: numpy
-def pnorm(a: npt.ArrayLike) -> np.ndarray:
+def pnorm(a: npt.NDArray) -> npt.NDArray:
     """Euclidean norm of each vector in a matrix
 
     Parameters
@@ -199,8 +196,7 @@ def dihedral(ab: npt.ArrayLike, bc: npt.ArrayLike, cd: npt.ArrayLike) -> float:
     return (x if stp(ab, bc, cd) <= 0.0 else -x)
 
 
-# typing: numpy
-def sarrus_det(matrix: np.ndarray) -> Union[float, np.ndarray]:
+def sarrus_det(matrix: npt.NDArray) -> Union[float, npt.NDArray]:
     """Computes the determinant of a 3x3 matrix according to the
     `rule of Sarrus`_.
 
@@ -245,8 +241,7 @@ def sarrus_det(matrix: np.ndarray) -> Union[float, np.ndarray]:
     return _sarrus_det_multiple(m.reshape((-1, 3, 3))).reshape(shape[:-2])
 
 
-# typing: numpy
-def triclinic_box(x: npt.ArrayLike, y: npt.ArrayLike, z: npt.ArrayLike) -> np.ndarray:
+def triclinic_box(x: npt.ArrayLike, y: npt.ArrayLike, z: npt.ArrayLike) -> npt.NDArray:
     """Convert the three triclinic box vectors to
     ``[lx, ly, lz, alpha, beta, gamma]``.
 
@@ -308,8 +303,8 @@ def triclinic_box(x: npt.ArrayLike, y: npt.ArrayLike, z: npt.ArrayLike) -> np.nd
     return np.zeros(6, dtype=np.float32)
 
 
-# typing: numpy
-def triclinic_vectors(dimensions: npt.ArrayLike, dtype: npt.DTypeLike = np.float32) -> np.ndarray:
+def triclinic_vectors(dimensions: npt.ArrayLike,
+                      dtype: npt.DTypeLike = np.float32) -> npt.NDArray:
     """Convert ``[lx, ly, lz, alpha, beta, gamma]`` to a triclinic matrix
     representation.
 
diff --git a/package/MDAnalysis/lib/pkdtree.py b/package/MDAnalysis/lib/pkdtree.py
index fa4f8ace0d5..33d48fa1ddf 100644
--- a/package/MDAnalysis/lib/pkdtree.py
+++ b/package/MDAnalysis/lib/pkdtree.py
@@ -159,8 +159,7 @@ def set_coords(self, coords: npt.ArrayLike, cutoff: Optional[float] = None) -> N
             self.ckdt = cKDTree(self.coords, self.leafsize)
         self._built = True
 
-    #  typing: numpy
-    def search(self, centers: npt.ArrayLike, radius: float) -> np.ndarray:
+    def search(self, centers: npt.ArrayLike, radius: float) -> npt.NDArray:
         """Search all points within radius from centers and their periodic images.
 
         All the centers coordinates are wrapped around the central cell
@@ -210,8 +209,7 @@ def search(self, centers: npt.ArrayLike, radius: float) -> np.ndarray:
         self._indices = np.asarray(unique_int_1d(self._indices))
         return self._indices
 
-    # typing: numpy
-    def get_indices(self) -> np.ndarray:
+    def get_indices(self) -> npt.NDArray:
         """Return the neighbors from the last query.
 
         Returns
@@ -221,8 +219,7 @@ def get_indices(self) -> np.ndarray:
         """
         return self._indices
 
-    # typing: numpy
-    def search_pairs(self, radius: float) -> np.ndarray:
+    def search_pairs(self, radius: float) -> npt.NDArray:
         """Search all the pairs within a specified radius
 
         Parameters
diff --git a/package/pyproject.toml b/package/pyproject.toml
index ad168894fd4..e4336b7ee81 100644
--- a/package/pyproject.toml
+++ b/package/pyproject.toml
@@ -6,10 +6,10 @@ requires = [
     # lowest NumPy we can use for a given Python,
     # In part adapted from: https://github.com/scipy/oldest-supported-numpy/blob/main/setup.cfg
     # except for more exotic platform (Mac Arm flavors)
-    # aarch64, AIX, s390x all support < 1.20 so we can safely pin to this
+    # aarch64, AIX, s390x all support < 1.21 so we can safely pin to this
     # Note: MDA does not build with PyPy so we do not support it in the build system
-    "numpy==1.20.0; python_version=='3.8' and (platform_machine!='arm64' or platform_system!='Darwin') and platform_python_implementation != 'PyPy'",
-    "numpy==1.20.0; python_version=='3.9' and (platform_machine!='arm64' or platform_system!='Darwin') and platform_python_implementation != 'PyPy'",
+    "numpy==1.21.0; python_version=='3.8' and (platform_machine!='arm64' or platform_system!='Darwin') and platform_python_implementation != 'PyPy'",
+    "numpy==1.21.0; python_version=='3.9' and (platform_machine!='arm64' or platform_system!='Darwin') and platform_python_implementation != 'PyPy'",
     # arm64 on darwin for py3.8+ requires numpy >=1.21.0
     "numpy==1.21.0; python_version=='3.8' and platform_machine=='arm64' and platform_system=='Darwin' and platform_python_implementation != 'PyPy'",
     "numpy==1.21.0; python_version=='3.9' and platform_machine=='arm64' and platform_system=='Darwin' and platform_python_implementation != 'PyPy'",
@@ -39,7 +39,7 @@ maintainers = [
 ]
 requires-python = ">=3.8"
 dependencies = [
-    'numpy>=1.20.0',
+    'numpy>=1.21.0',
     'biopython>=1.80',
     'networkx>=2.0',
     'GridDataFormats>=0.4.0',
diff --git a/package/requirements.txt b/package/requirements.txt
index 44c34303a9f..a9daf1375ff 100644
--- a/package/requirements.txt
+++ b/package/requirements.txt
@@ -11,7 +11,7 @@ mmtf-python
 msmb_theme==1.2.0
 netcdf4
 networkx
-numpy>=1.20.0
+numpy>=1.21.0
 packaging
 parmed
 pytest
diff --git a/package/setup.py b/package/setup.py
index 0c6a5446a86..186fe94161f 100755
--- a/package/setup.py
+++ b/package/setup.py
@@ -183,7 +183,7 @@ def get_numpy_include():
         import numpy as np
     except ImportError:
         print('*** package "numpy" not found ***')
-        print('MDAnalysis requires a version of NumPy (>=1.20.0), even for setup.')
+        print('MDAnalysis requires a version of NumPy (>=1.21.0), even for setup.')
         print('Please get it from http://numpy.scipy.org/ or install it through '
               'your package manager.')
         sys.exit(-1)
@@ -593,7 +593,7 @@ def long_description(readme):
     exts, cythonfiles = extensions(config)
 
     install_requires = [
-          'numpy>=1.20.0',
+          'numpy>=1.21.0',
           'biopython>=1.80',
           'networkx>=2.0',
           'GridDataFormats>=0.4.0',
@@ -641,7 +641,7 @@ def long_description(readme):
           # all standard requirements are available through PyPi and
           # typically can be installed without difficulties through setuptools
           setup_requires=[
-              'numpy>=1.20.0',
+              'numpy>=1.21.0',
               'packaging',
           ],
           install_requires=install_requires,