diff --git a/.github/workflows/python-package-conda.yml b/.github/workflows/python-package-conda.yml
index f448047a..0a035c29 100644
--- a/.github/workflows/python-package-conda.yml
+++ b/.github/workflows/python-package-conda.yml
@@ -10,7 +10,7 @@ on:
 
 env:
   latest_python: "3.13"
-  supported_pythons: '["3.9", "3.10", "3.11", "3.12", "3.13"]'
+  supported_pythons: '["3.10", "3.11", "3.12", "3.13"]'
   miniforge_version: "23.11.0-0"
   miniforge_variant: "Mambaforge"
 
diff --git a/.github/workflows/release.yml b/.github/workflows/release.yml
index f43dcb9f..3dc3f67e 100644
--- a/.github/workflows/release.yml
+++ b/.github/workflows/release.yml
@@ -6,7 +6,7 @@ on:
       - '*'
 
 env:
-  earliest_python: "3.9"
+  earliest_python: "3.10"
   latest_python: "3.13"
   miniforge_version: "23.11.0-0"
   miniforge_variant: "Mambaforge"
@@ -16,10 +16,10 @@ jobs:
     runs-on: ubuntu-latest
     steps:
       - uses: actions/checkout@v2
-      - name: Set up Python 3.9
+      - name: Set up Python 3.10
         uses: actions/setup-python@v2
         with:
-          python-version: 3.9
+          python-version: 3.10
       - name: Build distribution
         run: |
           # set version from '${{ github.ref_name }}'
diff --git a/.gitignore b/.gitignore
index d1f9eafa..b781e587 100644
--- a/.gitignore
+++ b/.gitignore
@@ -3,6 +3,10 @@
 *.pyc
 biom_format.egg-info/
 
+# Cython converted files
+########################
+*.c
+
 # Cython compiled files
 #######################
 *.so
diff --git a/ChangeLog.md b/ChangeLog.md
index b2f048e6..72e9a6aa 100644
--- a/ChangeLog.md
+++ b/ChangeLog.md
@@ -4,6 +4,10 @@ BIOM-Format ChangeLog
 biom-2.1.16-dev
 ---------------
 
+Important:
+
+* The underlying data structure has been migrated from SciPy sparse matrix (spmatrix) to sparse array (sparray). See PR[#993](https://github.com/biocore/biom-format/pull/993). This change does not affect any public-facing API of biom-format. However, be minded that the downstream analysis of data extracted from a BIOM table may be impacted by the API change in sparray vs. spmatrix. See SciPy's documentation: https://docs.scipy.org/doc/scipy/reference/sparse.html for details. This change requires SciPy >= 1.8.0. See its release note: https://docs.scipy.org/doc/scipy/release/1.8.0-notes.html for details.
+
 Performance improvements:
 
 * Decreased execution time of `import biom` by half with lazy imports. See PR[#987](https://github.com/biocore/biom-format/pull/987)
diff --git a/biom/_filter.pyx b/biom/_filter.pyx
index 879db01d..b840d8bc 100644
--- a/biom/_filter.pyx
+++ b/biom/_filter.pyx
@@ -116,7 +116,7 @@ def _filter(arr, ids, metadata, index, ids_to_keep, axis, invert):
         arr = arr.tocsr()
     elif axis == 1:
         arr = arr.tocsc()
-    fmt = arr.getformat()
+    fmt = arr.format
 
     cdef cnp.ndarray[cnp.uint8_t, ndim=1] bools
 
diff --git a/biom/_subsample.pyx b/biom/_subsample.pyx
index 4261ac92..fe79c55b 100644
--- a/biom/_subsample.pyx
+++ b/biom/_subsample.pyx
@@ -21,9 +21,9 @@ cdef _subsample_with_replacement(cnp.ndarray[cnp.float64_t, ndim=1] data,
 
     Parameters
     ----------
-    data : {csr_matrix, csc_matrix}.data
+    data : {csr_array, csc_array}.data
         A 1xM sparse vector data
-    indptr : {csr_matrix, csc_matrix}.indptr
+    indptr : {csr_array, csc_array}.indptr
         A 1xM sparse vector indptr
     n : int
         Number of items to subsample from `arr`
@@ -66,9 +66,9 @@ cdef _subsample_without_replacement(cnp.ndarray[cnp.float64_t, ndim=1] data,
 
     Parameters
     ----------
-    data : {csr_matrix, csc_matrix}.data
+    data : {csr_array, csc_array}.data
         A 1xM sparse vector data
-    indptr : {csr_matrix, csc_matrix}.indptr
+    indptr : {csr_array, csc_array}.indptr
         A 1xM sparse vector indptr
     n : int
         Number of items to subsample from `arr`
@@ -153,7 +153,7 @@ def subsample(arr, n, with_replacement, rng):
 
     Parameters
     ----------
-    arr : {csr_matrix, csc_matrix}
+    arr : {csr_array, csc_array}
         A 1xM sparse vector
     n : int
         Number of items to subsample from `arr`
diff --git a/biom/_transform.pyx b/biom/_transform.pyx
index 0d500042..78ab3a4f 100644
--- a/biom/_transform.pyx
+++ b/biom/_transform.pyx
@@ -20,7 +20,7 @@ def _transform(arr, ids, metadata, function, axis):
 
     Parameters
     ----------
-    arr : csr_matrix or csc_matrix
+    arr : csr_array or csc_array
         Matrix whose rows or columns (respectively) are to be
         transformed.
     ids : 1D array_like
diff --git a/biom/cli/installation_informer.py b/biom/cli/installation_informer.py
index 3fb9e4c4..2b98a851 100644
--- a/biom/cli/installation_informer.py
+++ b/biom/cli/installation_informer.py
@@ -8,6 +8,7 @@
 
 
 import sys
+from importlib.metadata import version, PackageNotFoundError
 
 import click
 
@@ -63,25 +64,25 @@ def _get_dependency_version_info():
     not_installed_msg = "Not installed"
 
     try:
-        from click import __version__ as click_lib_version
-    except ImportError:
+        click_lib_version = version("click")
+    except PackageNotFoundError:
         click_lib_version = not_installed_msg
 
     try:
-        from numpy import __version__ as numpy_lib_version
-    except ImportError:
+        numpy_lib_version = version("numpy")
+    except PackageNotFoundError:
         numpy_lib_version = ("ERROR: Not installed - this is required! "
                              "(This will also cause the BIOM library to "
                              "not be importable.)")
 
     try:
-        from scipy import __version__ as scipy_lib_version
-    except ImportError:
+        scipy_lib_version = version("scipy")
+    except PackageNotFoundError:
         scipy_lib_version = not_installed_msg
 
     try:
-        from h5py import __version__ as h5py_lib_version
-    except ImportError:
+        h5py_lib_version = version("h5py")
+    except PackageNotFoundError:
         h5py_lib_version = ("WARNING: Not installed - this is an optional "
                             "dependency. It is strongly recommended for "
                             "large datasets.")
@@ -97,8 +98,8 @@ def _get_package_info():
                         "numpy) - is it installed and in your "
                         "$PYTHONPATH?")
     try:
-        from biom import __version__ as biom_lib_version
-    except ImportError:
+        biom_lib_version = version("biom-format")
+    except PackageNotFoundError:
         biom_lib_version = import_error_msg
 
     return (("biom-format version", biom_lib_version),)
diff --git a/biom/table.py b/biom/table.py
index 2cd40967..35f161d7 100644
--- a/biom/table.py
+++ b/biom/table.py
@@ -180,8 +180,8 @@
 from collections import defaultdict
 from collections.abc import Hashable, Iterable
 from numpy import ndarray, asarray, zeros, newaxis
-from scipy.sparse import (coo_matrix, csc_matrix, csr_matrix, isspmatrix,
-                          vstack, hstack, dok_matrix)
+from scipy.sparse import (coo_array, csc_array, csr_array, issparse,
+                          vstack, hstack, dok_array)
 import re
 from biom.exception import (TableException, UnknownAxisError, UnknownIDError,
                             DisjointIDError)
@@ -479,7 +479,7 @@ def __init__(self, data, observation_ids, sample_ids,
         self.generated_by = generated_by
         self.format_version = __format_version__
 
-        if not isspmatrix(data):
+        if not issparse(data):
             shape = (len(observation_ids), len(sample_ids))
             input_is_dense = kwargs.get('input_is_dense', False)
             self._data = Table._to_sparse(data, input_is_dense=input_is_dense,
@@ -570,7 +570,7 @@ def _conv_to_self_type(self, vals, transpose=False, dtype=None):
         if dtype is None:
             dtype = self.dtype
 
-        if isspmatrix(vals):
+        if issparse(vals):
             return vals
         else:
             return Table._to_sparse(vals, transpose, dtype)
@@ -615,7 +615,7 @@ def _to_sparse(values, transpose=False, dtype=float, input_is_dense=False,
             return mat
         # the empty list
         elif isinstance(values, list) and len(values) == 0:
-            return coo_matrix((0, 0))
+            return coo_array((0, 0))
         # list of np vectors
         elif isinstance(values, list) and isinstance(values[0], ndarray):
             mat = list_nparray_to_sparse(values, dtype)
@@ -629,7 +629,7 @@ def _to_sparse(values, transpose=False, dtype=float, input_is_dense=False,
                 mat = mat.T
             return mat
         # list of scipy.sparse matrices, each representing a row in row order
-        elif isinstance(values, list) and isspmatrix(values[0]):
+        elif isinstance(values, list) and issparse(values[0]):
             mat = list_sparse_to_sparse(values, dtype)
             if transpose:
                 mat = mat.T
@@ -641,13 +641,13 @@ def _to_sparse(values, transpose=False, dtype=float, input_is_dense=False,
             return mat
         elif isinstance(values, list) and isinstance(values[0], list):
             if input_is_dense:
-                d = coo_matrix(values)
+                d = coo_array(values)
                 mat = coo_arrays_to_sparse((d.data, (d.row, d.col)),
                                            dtype=dtype, shape=shape)
             else:
                 mat = list_list_to_sparse(values, dtype, shape=shape)
             return mat
-        elif isspmatrix(values):
+        elif issparse(values):
             mat = values
             if transpose:
                 mat = mat.transpose()
@@ -879,9 +879,9 @@ def __getitem__(self, args):
 
         Returns
         -------
-        float or spmatrix
+        float or sparray
             A float is return if a specific element is specified, otherwise a
-            spmatrix object representing a vector of sparse data is returned.
+            sparray object representing a vector of sparse data is returned.
 
         Raises
         ------
@@ -924,7 +924,7 @@ def __getitem__(self, args):
             else:
                 raise IndexError("Can only handle full : slices per axis.")
         else:
-            if self._data.getformat() == 'coo':
+            if self._data.format == 'coo':
                 self._data = self._data.tocsr()
 
             return self._data[row, col]
@@ -944,7 +944,7 @@ def _get_row(self, row_idx):
 
         """
         self._data = self._data.tocsr()
-        return self._data.getrow(row_idx)
+        return self._data[row_idx, None, :]  # this makes a 2D row
 
     def _get_col(self, col_idx):
         """Return the column at ``col_idx``.
@@ -962,7 +962,7 @@ def _get_col(self, col_idx):
 
         """
         self._data = self._data.tocsc()
-        return self._data.getcol(col_idx)
+        return self._data[:, None, col_idx]  # this makes a 2D column
 
     def align_to_dataframe(self, metadata, axis='sample'):
         """ Aligns dataframe against biom table, only keeping common ids.
@@ -1202,7 +1202,7 @@ def transpose(self):
         sample_md_copy = deepcopy(self.metadata())
         obs_md_copy = deepcopy(self.metadata(axis='observation'))
 
-        if self._data.getformat() == 'lil':
+        if self._data.format == 'lil':
             # lil's transpose method doesn't have the copy kwarg, but all of
             # the others do.
             self._data = self._data.tocsr()
@@ -1910,8 +1910,8 @@ def data(self, id, axis='sample', dense=True):
 
         Returns
         -------
-        np.ndarray or scipy.sparse.spmatrix
-            np.ndarray if ``dense``, otherwise scipy.sparse.spmatrix
+        np.ndarray or scipy.sparse.sparray
+            np.ndarray if ``dense``, otherwise scipy.sparse.sparray
 
         Raises
         ------
@@ -2736,12 +2736,11 @@ def axis_update(offaxis, onaxis):
             dtype = np.float64 if one_to_many_mode == 'divide' else self.dtype
 
             if axis == 'observation':
-                new_data = dok_matrix((len(self.ids(axis='sample')),
-                                       len(new_md)),
-                                      dtype=dtype)
+                new_data = dok_array((len(self.ids(axis='sample')),
+                                      len(new_md)), dtype=dtype)
             else:
-                new_data = dok_matrix((len(self.ids(axis='observation')),
-                                       len(new_md)), dtype=dtype)
+                new_data = dok_array((len(self.ids(axis='observation')),
+                                      len(new_md)), dtype=dtype)
 
             # for each vector
             # for each bin in the metadata
@@ -2794,9 +2793,9 @@ def axis_update(offaxis, onaxis):
 
             # convert back to self type
             if axis == 'observation':
-                new_data = csr_matrix(new_data.T)
+                new_data = csr_array(new_data.T)
             else:
-                new_data = csc_matrix(new_data)
+                new_data = csc_array(new_data)
 
             data = self._conv_to_self_type(new_data)
         else:
@@ -3617,7 +3616,7 @@ def concat(self, others, axis='sample'):
                     shape = (n_axis, n_invaxis)
 
                 # create the padded matrix
-                zerod = csr_matrix(shape)
+                zerod = csr_array(shape)
                 tmp_mat = invstack([table.matrix_data, zerod])
 
                 # resolve invert axis ids and metadata
@@ -3732,8 +3731,8 @@ def _fast_merge(self, others):
             data[offset:offset + t_nnz] = coo.data
             offset += t_nnz
 
-        coo = coo_matrix((data, (rows, cols)),
-                         shape=(len(feature_order), len(sample_order)))
+        coo = coo_array((data, (rows, cols)),
+                        shape=(len(feature_order), len(sample_order)))
 
         return self.__class__(coo.tocsr(), feature_order, sample_order)
 
@@ -4139,12 +4138,12 @@ def from_hdf5(cls, h5grp, ids=None, axis='sample', parse_fs=None,
                 obs_ids = h5grp['observation/ids'][:]
                 samp_ids = axis_ids[to_keep]
                 shape = (len(obs_ids), len(to_keep))
-                mat = csc_matrix((data, indices, indptr), shape=shape)
+                mat = csc_array((data, indices, indptr), shape=shape)
             else:
                 samp_ids = h5grp['sample/ids'][:]
                 obs_ids = axis_ids[to_keep]
                 shape = (len(to_keep), len(samp_ids))
-                mat = csr_matrix((data, indices, indptr), shape=shape)
+                mat = csr_array((data, indices, indptr), shape=shape)
 
             # use a fixed width dtype
             obs_ids_dtype = 'U%d' % max([len(v) for v in obs_ids])
@@ -4234,7 +4233,7 @@ def _get_ids(source_ids, desired_ids):
                 else:
                     desired_ids = np.asarray(desired_ids)
                     # Get the index of the source ids to include
-                    idx = np.in1d(source_ids, desired_ids)
+                    idx = np.isin(source_ids, desired_ids)
                     # Retrieve only the ids that we are interested on
                     ids = source_ids[idx]
                     # Check that all desired ids have been found on source ids
@@ -4290,9 +4289,9 @@ def _subset_metadata(md, idx):
         cs = (data, indices, indptr)
 
         if axis == 'sample':
-            matrix = csc_matrix(cs, shape=shape)
+            matrix = csc_array(cs, shape=shape)
         else:
-            matrix = csr_matrix(cs, shape=shape)
+            matrix = csr_array(cs, shape=shape)
 
         t = Table(matrix, obs_ids, samp_ids, obs_md or None,
                   samp_md or None, type=type_, create_date=create_date,
@@ -4587,9 +4586,9 @@ def to_hdf5(self, h5grp, generated_by, compress=True, format_fs=None,
         References
         ----------
         .. [1] http://docs.scipy.org/doc/scipy-0.13.0/reference/generated/sci\
-py.sparse.csr_matrix.html
+py.sparse.csr_array.html
         .. [2] http://docs.scipy.org/doc/scipy-0.13.0/reference/generated/sci\
-py.sparse.csc_matrix.html
+py.sparse.csc_array.html
         .. [3] http://biom-format.org/documentation/format_versions/biom-2.1.\
 html
 
@@ -5066,10 +5065,11 @@ def is_num(item):
         samp_index = {s: i for i, s in enumerate(samp_order)}
 
         # fill the matrix
-        row = np.array([obs_index[obs] for obs in observations], dtype=int)
-        col = np.array([samp_index[samp] for samp in samples], dtype=int)
+        row = np.array([obs_index[obs] for obs in observations],
+                       dtype=np.int32)
+        col = np.array([samp_index[samp] for samp in samples], dtype=np.int32)
         data = np.asarray(values)
-        mat = coo_matrix((data, (row, col)))
+        mat = coo_array((data, (row, col)))
 
         return Table(mat, obs_order, samp_order)
 
@@ -5333,7 +5333,7 @@ def to_tsv(self, header_key=None, header_value=None,
 
 
 def coo_arrays_to_sparse(data, dtype=np.float64, shape=None):
-    """Map directly on to the coo_matrix constructor
+    """Map directly on to the coo_array constructor
 
     Parameters
     ----------
@@ -5352,14 +5352,14 @@ def coo_arrays_to_sparse(data, dtype=np.float64, shape=None):
     else:
         n_rows, n_cols = shape
 
-    # coo_matrix allows zeros to be added as data, and this affects
+    # coo_array allows zeros to be added as data, and this affects
     # nnz, items, and iteritems. Clean them out here, as this is
     # the only time these zeros can creep in.
-    # Note: coo_matrix allows duplicate entries; the entries will
+    # Note: coo_array allows duplicate entries; the entries will
     # be summed when converted. Not really sure how we want to
     # handle this generally within BIOM- I'm okay with leaving it
     # as undefined behavior for now.
-    matrix = coo_matrix(data, shape=(n_rows, n_cols), dtype=dtype)
+    matrix = coo_array(data, shape=(n_rows, n_cols), dtype=dtype)
     matrix = matrix.tocsr()
     matrix.eliminate_zeros()
     return matrix
@@ -5380,10 +5380,12 @@ def list_list_to_sparse(data, dtype=float, shape=None):
 
     Returns
     -------
-    scipy.csr_matrix
+    scipy.csr_array
         The newly generated matrix
     """
     rows, cols, values = zip(*data)
+    rows = np.asarray(rows, dtype=np.int32)
+    cols = np.asarray(cols, dtype=np.int32)
 
     if shape is None:
         n_rows = max(rows) + 1
@@ -5391,8 +5393,8 @@ def list_list_to_sparse(data, dtype=float, shape=None):
     else:
         n_rows, n_cols = shape
 
-    matrix = coo_matrix((values, (rows, cols)), shape=(n_rows, n_cols),
-                        dtype=dtype)
+    matrix = coo_array((values, (rows, cols)), shape=(n_rows, n_cols),
+                       dtype=dtype)
     matrix = matrix.tocsr()
     matrix.eliminate_zeros()
     return matrix
@@ -5410,27 +5412,30 @@ def nparray_to_sparse(data, dtype=float):
 
     Returns
     -------
-    scipy.csr_matrix
+    scipy.csr_array
         The newly generated matrix
     """
     if data.shape == (0,):
         # an empty vector. Note, this short circuit is necessary as calling
-        # csr_matrix([], shape=(0, 0), dtype=dtype) will result in a matrix
+        # csr_array([], shape=(0, 0), dtype=dtype) will result in a matrix
         # has a shape of (1, 0).
-        return csr_matrix((0, 0), dtype=dtype)
+        return csr_array((0, 0), dtype=dtype)
     elif data.shape in ((1, 0), (0, 1)) and data.size == 0:
         # an empty matrix. This short circuit is necessary for the same reason
         # as the empty vector. While a (1, 0) matrix is _empty_, this does
         # confound code that assumes that (1, 0) means there might be metadata
         # or IDs associated with that singular row
-        return csr_matrix((0, 0), dtype=dtype)
+        return csr_array((0, 0), dtype=dtype)
     elif len(data.shape) == 1:
         # a vector
-        shape = (1, data.shape[0])
+        n = data.size
+        row = np.zeros(n, dtype=np.int32)
+        col = np.arange(n, dtype=np.int32)
+        matrix = coo_array((data, (row, col)), shape=(1, n), dtype=dtype)
     else:
-        shape = data.shape
+        # a 2D array
+        matrix = coo_array(data, dtype=dtype)
 
-    matrix = coo_matrix(data, shape=shape, dtype=dtype)
     matrix = matrix.tocsr()
     matrix.eliminate_zeros()
     return matrix
@@ -5448,10 +5453,10 @@ def list_nparray_to_sparse(data, dtype=float):
 
     Returns
     -------
-    scipy.csr_matrix
+    scipy.csr_array
         The newly generated matrix
     """
-    matrix = coo_matrix(data, shape=(len(data), len(data[0])), dtype=dtype)
+    matrix = coo_array(data, shape=(len(data), len(data[0])), dtype=dtype)
     matrix = matrix.tocsr()
     matrix.eliminate_zeros()
     return matrix
@@ -5469,10 +5474,10 @@ def list_sparse_to_sparse(data, dtype=float):
 
     Returns
     -------
-    scipy.csr_matrix
+    scipy.csr_array
         The newly generated matrix
     """
-    if isspmatrix(data[0]):
+    if issparse(data[0]):
         if data[0].shape[0] > data[0].shape[1]:
             n_cols = len(data)
             n_rows = data[0].shape[0]
@@ -5489,8 +5494,7 @@ def list_sparse_to_sparse(data, dtype=float):
             n_rows = len(data)
 
     data = vstack(data)
-    matrix = coo_matrix(data, shape=(n_rows, n_cols),
-                        dtype=dtype)
+    matrix = coo_array(data, shape=(n_rows, n_cols), dtype=dtype)
     matrix = matrix.tocsr()
     matrix.eliminate_zeros()
     return matrix
@@ -5508,10 +5512,10 @@ def list_dict_to_sparse(data, dtype=float):
 
     Returns
     -------
-    scipy.csr_matrix
+    scipy.csr_array
         The newly generated matrix
     """
-    if isspmatrix(data[0]):
+    if issparse(data[0]):
         if data[0].shape[0] > data[0].shape[1]:
             is_col = True
             n_cols = len(data)
@@ -5546,8 +5550,10 @@ def list_dict_to_sparse(data, dtype=float):
                 cols.append(col_idx)
                 vals.append(val)
 
-    matrix = coo_matrix((vals, (rows, cols)), shape=(n_rows, n_cols),
-                        dtype=dtype)
+    rows = np.asarray(rows, dtype=np.int32)
+    cols = np.asarray(cols, dtype=np.int32)
+    matrix = coo_array((vals, (rows, cols)), shape=(n_rows, n_cols),
+                       dtype=dtype)
     matrix = matrix.tocsr()
     matrix.eliminate_zeros()
     return matrix
@@ -5565,7 +5571,7 @@ def dict_to_sparse(data, dtype=float, shape=None):
 
     Returns
     -------
-    scipy.csr_matrix
+    scipy.csr_array
         The newly generated matrix
     """
     if shape is None:
@@ -5582,5 +5588,7 @@ def dict_to_sparse(data, dtype=float, shape=None):
         cols.append(c)
         vals.append(v)
 
+    rows = np.asarray(rows, dtype=np.int32)
+    cols = np.asarray(cols, dtype=np.int32)
     return coo_arrays_to_sparse((vals, (rows, cols)),
                                 shape=(n_rows, n_cols), dtype=dtype)
diff --git a/biom/tests/test_table.py b/biom/tests/test_table.py
index f22283f6..feaa0700 100644
--- a/biom/tests/test_table.py
+++ b/biom/tests/test_table.py
@@ -16,7 +16,7 @@
 
 import numpy.testing as npt
 import numpy as np
-from scipy.sparse import lil_matrix, csr_matrix, csc_matrix
+from scipy.sparse import lil_array, csr_array, csc_array
 import scipy.sparse
 import pandas.testing as pdt
 import pandas as pd
@@ -115,7 +115,7 @@ def test_head_zero_or_neg(self):
             example_table.head(5, 0)
 
     def test_remove_empty_sample(self):
-        wrn = "Changing the sparsity structure of a csr_matrix is expensive. lil_matrix is more efficient."  # noqa
+        wrn = "Changing the sparsity structure of a csr_array is expensive. lil_array is more efficient."  # noqa
         with warnings.catch_warnings():
             warnings.filterwarnings("ignore", message=wrn)
             t = example_table.copy()
@@ -125,7 +125,7 @@ def test_remove_empty_sample(self):
             self.assertEqual(t, exp)
 
     def test_remove_empty_obs(self):
-        wrn = "Changing the sparsity structure of a csr_matrix is expensive. lil_matrix is more efficient."  # noqa
+        wrn = "Changing the sparsity structure of a csr_array is expensive. lil_array is more efficient."  # noqa
         with warnings.catch_warnings():
             warnings.filterwarnings("ignore", message=wrn)
             t = example_table.copy()
@@ -136,7 +136,7 @@ def test_remove_empty_obs(self):
             self.assertEqual(t, exp)
 
     def test_remove_empty_both(self):
-        wrn = "Changing the sparsity structure of a csr_matrix is expensive. lil_matrix is more efficient."  # noqa
+        wrn = "Changing the sparsity structure of a csr_array is expensive. lil_array is more efficient."  # noqa
         with warnings.catch_warnings():
             warnings.filterwarnings("ignore", message=wrn)
             t = example_table.copy()
@@ -446,7 +446,7 @@ def test_table_sparse_list_list(self):
         # list list test
         samp_ids = range(3)
         obs_ids = range(2)
-        exp_data = lil_matrix((2, 3))
+        exp_data = lil_array((2, 3))
         exp_data[0, 1] = 5
         exp_data[1, 2] = 10
         exp = Table(exp_data, obs_ids, samp_ids)
@@ -1594,8 +1594,8 @@ def test_add_group_metadata_w_existing_metadata(self):
              'tree': ('newick', '((4:0.1,5:0.1):0.2,(6:0.1,7:0.1):0.2):0.3;')})
 
     def test_to_dataframe(self):
-        mat = csr_matrix(np.array([[0.0, 1.0, 2.0],
-                                   [3.0, 4.0, 5.0]]))
+        mat = csr_array(np.array([[0.0, 1.0, 2.0],
+                                  [3.0, 4.0, 5.0]]))
         exp = pd.DataFrame.sparse.from_spmatrix(mat,
                                                 index=['O1', 'O2'],
                                                 columns=['S1', 'S2', 'S3'])
@@ -1609,7 +1609,7 @@ def test_to_dataframe(self):
 
     def test_to_dataframe_is_sparse(self):
         df = example_table.to_dataframe()
-        density = (float(example_table.matrix_data.getnnz()) /
+        density = (float(example_table.matrix_data.nnz) /
                    np.prod(example_table.shape))
         df_density = (df.values > 0).sum().sum() / np.prod(df.shape)
         assert np.allclose(df_density, density)
@@ -2128,21 +2128,21 @@ def test_is_empty(self):
 
     def test_convert_vector_to_dense(self):
         """Properly converts ScipySparseMat vectors to dense numpy repr."""
-        input_row = lil_matrix((1, 3))
+        input_row = lil_array((1, 3))
         input_row[(0, 0)] = 1
         input_row[(0, 2)] = 3
         exp = np.array([1, 0, 3])
         obs = self.row_vec._to_dense(input_row)
         npt.assert_array_equal(obs, exp)
 
-        input_row = lil_matrix((3, 1))
+        input_row = lil_array((3, 1))
         input_row[(0, 0)] = 1
         input_row[(2, 0)] = 3
         exp = np.array([1, 0, 3])
         obs = self.row_vec._to_dense(input_row)
         npt.assert_array_equal(obs, exp)
 
-        input_row = lil_matrix((1, 1))
+        input_row = lil_array((1, 1))
         input_row[(0, 0)] = 42
         exp = np.array([42])
         obs = self.single_ele._to_dense(input_row)
@@ -2190,7 +2190,7 @@ def test_nnz(self):
         self.assertEqual(self.explicit_zeros.nnz, 4)
 
     def test_nnz_issue_727(self):
-        wrn = "Changing the sparsity structure of a csr_matrix is expensive. lil_matrix is more efficient."  # noqa
+        wrn = "Changing the sparsity structure of a csr_array is expensive. lil_array is more efficient."  # noqa
         with warnings.catch_warnings():
             warnings.filterwarnings("ignore", message=wrn)
             tab = Table(np.array([[0, 1], [0, 0]]), ['a', 'b'], ['1', '2'])
@@ -2206,7 +2206,7 @@ def test_get_row(self):
             with self.assertRaises(IndexError):
                 self.nulls[i]._get_row(0)
 
-        exp = lil_matrix((1, 3))
+        exp = lil_array((1, 3))
         exp[(0, 0)] = 1
         exp[(0, 2)] = 2
 
@@ -2222,7 +2222,7 @@ def test_get_col(self):
             with self.assertRaises(IndexError):
                 self.nulls[i]._get_col(0)
 
-        exp = lil_matrix((2, 1))
+        exp = lil_array((2, 1))
         exp[(0, 0)] = 1
         exp[(1, 0)] = 3
 
@@ -2303,7 +2303,7 @@ def test_getitem_2(self):
         with self.assertRaises(IndexError):
             self.empty[0, 0:1]
 
-        exp = lil_matrix((2, 1))
+        exp = lil_array((2, 1))
         obs = self.empty[:, 0]
         self.assertEqual((obs != exp).sum(), 0)
 
@@ -2509,14 +2509,14 @@ def test_update_ids_cache_bug(self):
         exp_index = {'x': 0, 'y': 1}
         self.assertEqual(obs._sample_index, exp_index)
 
-    def test_other_spmatrix_type(self):
+    def test_other_sparray_type(self):
         ss = scipy.sparse
-        for c in [ss.lil_matrix, ss.bsr_matrix, ss.coo_matrix, ss.dia_matrix,
-                  ss.dok_matrix, ss.csc_matrix, ss.csr_matrix]:
+        for c in [ss.lil_array, ss.bsr_array, ss.coo_array, ss.dia_array,
+                  ss.dok_array, ss.csc_array, ss.csr_array]:
             mat = c((2, 2))
             t = Table(mat, ['a', 'b'], [1, 2])
             self.assertTrue(isinstance(t.matrix_data,
-                                       (csr_matrix, csc_matrix)))
+                                       (csr_array, csc_array)))
 
     def test_sort_order(self):
         """sorts tables by arbitrary order"""
@@ -2754,14 +2754,14 @@ def test_data(self):
 
     def test_data_sparse(self):
         # Returns observations for a given sample
-        exp = csc_matrix(np.array([[5], [7]]))
+        exp = csc_array(np.array([[5], [7]]))
         obs = self.st1.data('a', 'sample', dense=False)
         self.assertEqual((obs != exp).nnz, 0)
         with self.assertRaises(UnknownIDError):
             self.st1.data('asdasd', 'sample')
 
         # Returns samples for a given observation
-        exp = csr_matrix(np.array([5, 6]))
+        exp = csr_array(np.array([[5, 6]]))
         obs = self.st1.data('1', 'observation', dense=False)
         self.assertEqual((obs != exp).nnz, 0)
         with self.assertRaises(UnknownIDError):
@@ -2788,7 +2788,7 @@ def test_delimited_self(self):
 
     def test_conv_to_self_type(self):
         """Should convert other to sparse type"""
-        exp = lil_matrix((2, 2))
+        exp = lil_array((2, 2))
         exp[(0, 0)] = 5
         exp[(0, 1)] = 6
         exp[(1, 0)] = 7
@@ -2796,7 +2796,7 @@ def test_conv_to_self_type(self):
         obs = self.st1._conv_to_self_type(self.vals)
         self.assertEqual((obs != exp).sum(), 0)
 
-        exp = lil_matrix((2, 2))
+        exp = lil_array((2, 2))
         exp[(0, 0)] = 5
         exp[(0, 1)] = 7
         exp[(1, 0)] = 6
@@ -2805,7 +2805,7 @@ def test_conv_to_self_type(self):
         self.assertEqual((obs != exp).sum(), 0)
 
         # passing a single vector
-        exp = lil_matrix((1, 3))
+        exp = lil_array((1, 3))
         exp[(0, 0)] = 2
         exp[(0, 1)] = 0
         exp[(0, 2)] = 3
@@ -2813,7 +2813,7 @@ def test_conv_to_self_type(self):
         self.assertEqual((obs != exp).sum(), 0)
 
         # passing a list of dicts
-        exp = lil_matrix((2, 3))
+        exp = lil_array((2, 3))
         exp[(0, 0)] = 5
         exp[(0, 1)] = 6
         exp[(0, 2)] = 7
@@ -2826,20 +2826,20 @@ def test_conv_to_self_type(self):
 
     def test_to_dense(self):
         """Should convert a self styled vector to numpy type"""
-        input_row = lil_matrix((1, 3))
+        input_row = lil_array((1, 3))
         input_row[(0, 0)] = 10
         exp = np.array([10.0, 0, 0])
         obs = self.st1._to_dense(input_row)
         npt.assert_equal(obs, exp)
 
-        input_col = lil_matrix((3, 1))
+        input_col = lil_array((3, 1))
         input_col[(0, 0)] = 12
         exp = np.array([12.0, 0, 0])
         obs = self.st1._to_dense(input_col)
         npt.assert_equal(obs, exp)
 
         # 1x1
-        input_vec = lil_matrix((1, 1))
+        input_vec = lil_array((1, 1))
         input_vec[(0, 0)] = 42
         exp = np.array([42.0])
         obs = self.st1._to_dense(input_vec)
@@ -2851,8 +2851,8 @@ def test_iter_data_dense(self):
         npt.assert_equal(obs, exp)
 
     def test_iter_data_sparse(self):
-        exp = [csr_matrix(np.array([5, 7])),
-               csr_matrix(np.array([6, 8]))]
+        exp = [csr_array(np.array([[5, 7]])),
+               csr_array(np.array([[6, 8]]))]
         obs = list(self.st1.iter_data(dense=False))
         for o, e in zip(obs, exp):
             self.assertTrue((o != e).nnz == 0)
@@ -2895,8 +2895,8 @@ def test_iter(self):
 
     def test_iter_obs(self):
         """Iterate over observations of sparse matrix"""
-        r1 = lil_matrix((1, 2))
-        r2 = lil_matrix((1, 2))
+        r1 = lil_array((1, 2))
+        r2 = lil_array((1, 2))
         r1[(0, 0)] = 5
         r1[(0, 1)] = 6
         r2[(0, 0)] = 7
@@ -2910,8 +2910,8 @@ def test_iter_obs(self):
 
     def test_iter_samp(self):
         """Iterate over samples of sparse matrix"""
-        c1 = lil_matrix((1, 2))
-        c2 = lil_matrix((1, 2))
+        c1 = lil_array((1, 2))
+        c2 = lil_array((1, 2))
         c1[(0, 0)] = 5
         c1[(0, 1)] = 7
         c2[(0, 0)] = 6
@@ -3088,8 +3088,8 @@ def test_filter_general_sample(self):
         def f(vals, id_, md):
             return id_ == 'a'
 
-        values = csr_matrix(np.array([[5.],
-                                      [7.]]))
+        values = csr_array(np.array([[5.],
+                                     [7.]]))
         exp_table = Table(values, ['1', '2'], ['a'],
                           [{'taxonomy': ['k__a', 'p__b']},
                            {'taxonomy': ['k__a', 'p__c']}],
@@ -3109,7 +3109,7 @@ def test_filter_general_observation(self):
         def f(vals, id_, md):
             return md['taxonomy'][1] == 'p__c'
 
-        values = csr_matrix(np.array([[7., 8.]]))
+        values = csr_array(np.array([[7., 8.]]))
         exp_table = Table(values, ['2'], ['a', 'b'],
                           [{'taxonomy': ['k__a', 'p__c']}],
                           [{'barcode': 'aatt'}, {'barcode': 'ttgg'}])
@@ -3126,8 +3126,8 @@ def test_filter_sample_id(self):
         def f(vals, id_, md):
             return id_ == 'a'
 
-        values = csr_matrix(np.array([[5.],
-                                      [7.]]))
+        values = csr_array(np.array([[5.],
+                                     [7.]]))
         exp_table = Table(values, ['1', '2'], ['a'],
                           [{'taxonomy': ['k__a', 'p__b']},
                            {'taxonomy': ['k__a', 'p__c']}],
@@ -3140,8 +3140,8 @@ def f(vals, id_, md):
     def test_filter_sample_metadata(self):
         def f(vals, id_, md):
             return md['barcode'] == 'ttgg'
-        values = csr_matrix(np.array([[6.],
-                                      [8.]]))
+        values = csr_array(np.array([[6.],
+                                     [8.]]))
         exp_table = Table(values, ['1', '2'], ['b'],
                           [{'taxonomy': ['k__a', 'p__b']},
                            {'taxonomy': ['k__a', 'p__c']}],
@@ -3153,8 +3153,8 @@ def f(vals, id_, md):
     def test_filter_sample_invert(self):
         def f(vals, id_, md):
             return md['barcode'] == 'aatt'
-        values = csr_matrix(np.array([[6.],
-                                      [8.]]))
+        values = csr_array(np.array([[6.],
+                                     [8.]]))
         exp_table = Table(values, ['1', '2'], ['b'],
                           [{'taxonomy': ['k__a', 'p__b']},
                            {'taxonomy': ['k__a', 'p__c']}],
@@ -3174,7 +3174,7 @@ def test_filter_observations_id(self):
         def f(vals, id_, md):
             return id_ == '1'
 
-        values = csr_matrix(np.array([[5., 6.]]))
+        values = csr_array(np.array([[5., 6.]]))
         exp_table = Table(values, ['1'], ['a', 'b'],
                           [{'taxonomy': ['k__a', 'p__b']}],
                           [{'barcode': 'aatt'}, {'barcode': 'ttgg'}])
@@ -3186,7 +3186,7 @@ def test_filter_observations_metadata(self):
         def f(vals, id_, md):
             return md['taxonomy'][1] == 'p__c'
 
-        values = csr_matrix(np.array([[7., 8.]]))
+        values = csr_array(np.array([[7., 8.]]))
         exp_table = Table(values, ['2'], ['a', 'b'],
                           [{'taxonomy': ['k__a', 'p__c']}],
                           [{'barcode': 'aatt'}, {'barcode': 'ttgg'}])
@@ -3198,7 +3198,7 @@ def test_filter_observations_invert(self):
         def f(vals, id_, md):
             return md['taxonomy'][1] == 'p__c'
 
-        values = csr_matrix(np.array([[5., 6.]]))
+        values = csr_array(np.array([[5., 6.]]))
         exp_table = Table(values, ['1'], ['a', 'b'],
                           [{'taxonomy': ['k__a', 'p__b']}],
                           [{'barcode': 'aatt'}, {'barcode': 'ttgg'}])
@@ -4525,12 +4525,12 @@ class SupportTests2(TestCase):
     def test_coo_arrays_to_sparse(self):
         """convert (values, (row, col)) to scipy"""
         n_rows, n_cols = 3, 4
-        exp_d = lil_matrix((n_rows, n_cols))
+        exp_d = lil_array((n_rows, n_cols))
         exp_d[(0, 0)] = 10
         exp_d[(1, 3)] = 5
         exp_d[(2, 1)] = 2
         exp_d = exp_d.tocoo()
-        exp = lil_matrix((n_rows, n_cols))
+        exp = lil_array((n_rows, n_cols))
         exp[(0, 0)] = 10
         exp[(1, 3)] = 5
         exp[(2, 1)] = 2
@@ -4543,7 +4543,7 @@ def test_coo_arrays_to_sparse(self):
     def test_list_list_to_sparse(self):
         """convert [[row,col,value], ...] to scipy"""
         input = [[0, 0, 1], [1, 1, 5.0], [0, 2, 6]]
-        exp = lil_matrix((2, 3))
+        exp = lil_array((2, 3))
         exp[(0, 0)] = 1.0
         exp[(1, 1)] = 5.0
         exp[(0, 2)] = 6
@@ -4553,7 +4553,7 @@ def test_list_list_to_sparse(self):
     def test_nparray_to_sparse(self):
         """Convert nparray to sparse"""
         input = np.array([[1, 2, 3, 4], [-1, 6, 7, 8], [9, 10, 11, 12]])
-        exp = lil_matrix((3, 4))
+        exp = lil_array((3, 4))
         exp[(0, 0)] = 1
         exp[(0, 1)] = 2
         exp[(0, 2)] = 3
@@ -4572,7 +4572,7 @@ def test_nparray_to_sparse(self):
     def test_list_dict_to_sparse(self):
         """Take a list of dicts and condense down to a single dict"""
         input = [{(0, 0): 10, (0, 1): 2}, {(1, 2): 15}, {(0, 3): 7}]
-        exp = lil_matrix((3, 4))
+        exp = lil_array((3, 4))
         exp[(0, 0)] = 10
         exp[(0, 1)] = 2
         exp[(1, 2)] = 15
@@ -4583,7 +4583,7 @@ def test_list_dict_to_sparse(self):
     def test_dict_to_sparse(self):
         """Take a dict and convert to sparse"""
         input = {(0, 1): 5, (1, 0): 2, (2, 1): 6}
-        exp = lil_matrix((3, 2))
+        exp = lil_array((3, 2))
         exp[(0, 1)] = 5
         exp[(1, 0)] = 2
         exp[(2, 1)] = 6
@@ -4594,7 +4594,7 @@ def test_to_sparse(self):
         """Convert to expected sparse types"""
         vals = {(0, 0): 5, (0, 1): 6, (1, 0): 7, (1, 1): 8}
         obs = Table._to_sparse(vals)
-        exp = lil_matrix((2, 2))
+        exp = lil_array((2, 2))
         exp[(0, 0)] = 5
         exp[(0, 1)] = 6
         exp[(1, 0)] = 7
@@ -4604,21 +4604,21 @@ def test_to_sparse(self):
         input = {(0, 1): 5, (10, 8): -1.23}
         input_transpose = {(1, 0): 5, (8, 10): -1.23}
 
-        exp = lil_matrix((11, 9))
+        exp = lil_array((11, 9))
         exp[(0, 1)] = 5
         exp[(10, 8)] = -1.23
         obs = Table._to_sparse(input)
         self.assertEqual((obs != exp).sum(), 0)
 
         # test transpose
-        exp = lil_matrix((9, 11))
+        exp = lil_array((9, 11))
         exp[(1, 0)] = 5
         exp[(8, 10)] = -1.23
         obs = Table._to_sparse(input_transpose)
         self.assertEqual((obs != exp).sum(), 0)
 
         # passing a list of dicts, transpose
-        exp = lil_matrix((3, 2))
+        exp = lil_array((3, 2))
         exp[(0, 0)] = 5.0
         exp[(1, 0)] = 6.0
         exp[(2, 0)] = 7.0
@@ -4629,19 +4629,19 @@ def test_to_sparse(self):
                                 {(0, 1): 8, (1, 1): 9, (2, 1): 10}])
         self.assertEqual((obs != exp).sum(), 0)
 
-        # passing a list of lil_matrix
-        exp = lil_matrix((2, 3))
+        # passing a list of lil_array
+        exp = lil_array((2, 3))
         exp[(0, 0)] = 5
         exp[(0, 1)] = 6
         exp[(0, 2)] = 7
         exp[(1, 0)] = 8
         exp[(1, 1)] = 9
         exp[(1, 2)] = 10
-        row1 = lil_matrix((1, 3))
+        row1 = lil_array((1, 3))
         row1[(0, 0)] = 5
         row1[(0, 1)] = 6
         row1[(0, 2)] = 7
-        row2 = lil_matrix((1, 3))
+        row2 = lil_array((1, 3))
         row2[(0, 0)] = 8
         row2[(0, 1)] = 9
         row2[(0, 2)] = 10
@@ -4649,14 +4649,14 @@ def test_to_sparse(self):
         self.assertEqual((obs != exp).sum(), 0)
 
         # test empty set
-        exp = lil_matrix((0, 0))
+        exp = lil_array((0, 0))
         obs = Table._to_sparse([])
         self.assertEqual((obs != exp).sum(), 0)
 
     def test_list_nparray_to_sparse(self):
         """lists of nparrays to sparse"""
         ins = [np.array([0, 2, 1, 0]), np.array([1, 0, 0, 1])]
-        exp = lil_matrix((2, 4))
+        exp = lil_array((2, 4))
         exp[(0, 1)] = 2
         exp[(0, 2)] = 1
         exp[(1, 0)] = 1
@@ -4665,13 +4665,13 @@ def test_list_nparray_to_sparse(self):
         self.assertEqual((obs != exp).sum(), 0)
 
     def test_list_sparse_to_sparse(self):
-        """list of lil_matrix to sparse"""
-        ins = [lil_matrix((1, 4)), lil_matrix((1, 4))]
+        """list of lil_array to sparse"""
+        ins = [lil_array((1, 4)), lil_array((1, 4))]
         ins[0][0, 0] = 5
         ins[0][0, 1] = 10
         ins[1][0, 2] = 1
         ins[1][0, 3] = 2
-        exp = lil_matrix((2, 4))
+        exp = lil_array((2, 4))
         exp[0, 0] = 5
         exp[0, 1] = 10
         exp[1, 2] = 1
diff --git a/ci/conda_requirements.txt b/ci/conda_requirements.txt
index c07742dc..8f954770 100644
--- a/ci/conda_requirements.txt
+++ b/ci/conda_requirements.txt
@@ -1,6 +1,6 @@
 numpy >= 1.9.2
 pandas >= 0.20.0
-scipy >= 1.3.1
+scipy >= 1.8.0
 h5py >= 2.2.0
 anndata
 click
diff --git a/setup.cfg b/setup.cfg
index cea11b22..50eea740 100644
--- a/setup.cfg
+++ b/setup.cfg
@@ -5,4 +5,4 @@ test=pytest
 exclude=biom/tests/long_lines.py
 
 [options]
-python_requires = >=3.6
+python_requires = >=3.10
diff --git a/setup.py b/setup.py
index 0121adae..79329064 100644
--- a/setup.py
+++ b/setup.py
@@ -56,7 +56,6 @@
     Programming Language :: Python
     Programming Language :: Python :: 3
     Programming Language :: Python :: 3 :: Only
-    Programming Language :: Python :: 3.9
     Programming Language :: Python :: 3.10
     Programming Language :: Python :: 3.11
     Programming Language :: Python :: 3.12
@@ -85,7 +84,7 @@
 install_requires = [
     "click",
     "numpy >= 1.9.2",
-    "scipy >= 1.3.1",
+    "scipy >= 1.8.0",
     'pandas >= 0.20.0',
     "h5py",
 ]