diff --git a/Tests/test_image_resize.py b/Tests/test_image_resize.py
index d9ddf5009..8548fb5da 100644
--- a/Tests/test_image_resize.py
+++ b/Tests/test_image_resize.py
@@ -44,9 +44,19 @@ class TestImagingCoreResize:
             self.resize(hopper("1"), (15, 12), Image.Resampling.BILINEAR)
         with pytest.raises(ValueError):
             self.resize(hopper("P"), (15, 12), Image.Resampling.BILINEAR)
-        with pytest.raises(ValueError):
-            self.resize(hopper("I;16"), (15, 12), Image.Resampling.BILINEAR)
-        for mode in ["L", "I", "F", "RGB", "RGBA", "CMYK", "YCbCr"]:
+        for mode in [
+            "L",
+            "I",
+            "I;16",
+            "I;16L",
+            "I;16B",
+            "I;16N",
+            "F",
+            "RGB",
+            "RGBA",
+            "CMYK",
+            "YCbCr",
+        ]:
             im = hopper(mode)
             r = self.resize(im, (15, 12), Image.Resampling.BILINEAR)
             assert r.mode == mode
@@ -305,14 +315,14 @@ class TestImageResize:
             im = im.resize((64, 64))
             assert im.size == (64, 64)
 
-    @pytest.mark.parametrize("mode", ("L", "RGB", "I", "F"))
+    @pytest.mark.parametrize(
+        "mode", ("L", "RGB", "I", "I;16", "I;16L", "I;16B", "I;16N", "F")
+    )
     def test_default_filter_bicubic(self, mode: str) -> None:
         im = hopper(mode)
         assert im.resize((20, 20), Image.Resampling.BICUBIC) == im.resize((20, 20))
 
-    @pytest.mark.parametrize(
-        "mode", ("1", "P", "I;16", "I;16L", "I;16B", "BGR;15", "BGR;16")
-    )
+    @pytest.mark.parametrize("mode", ("1", "P", "BGR;15", "BGR;16"))
     def test_default_filter_nearest(self, mode: str) -> None:
         im = hopper(mode)
         assert im.resize((20, 20), Image.Resampling.NEAREST) == im.resize((20, 20))
diff --git a/docs/releasenotes/11.0.0.rst b/docs/releasenotes/11.0.0.rst
index 36334e39f..1218c5014 100644
--- a/docs/releasenotes/11.0.0.rst
+++ b/docs/releasenotes/11.0.0.rst
@@ -119,10 +119,11 @@ Specific WebP Feature Checks
 API Changes
 ===========
 
-TODO
-^^^^
+Default resampling filter for I;16* image modes
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 
-TODO
+The default resampling filter for I;16, I;16L, I;16B and I;16N has been changed from
+``Image.NEAREST`` to ``Image.BICUBIC``, to match the majority of modes.
 
 API Additions
 =============
diff --git a/src/PIL/Image.py b/src/PIL/Image.py
index 3f94cef38..bf9079f8b 100644
--- a/src/PIL/Image.py
+++ b/src/PIL/Image.py
@@ -2278,8 +2278,8 @@ class Image:
            :py:data:`Resampling.BILINEAR`, :py:data:`Resampling.HAMMING`,
            :py:data:`Resampling.BICUBIC` or :py:data:`Resampling.LANCZOS`.
            If the image has mode "1" or "P", it is always set to
-           :py:data:`Resampling.NEAREST`. If the image mode specifies a number
-           of bits, such as "I;16", then the default filter is
+           :py:data:`Resampling.NEAREST`. If the image mode is "BGR;15",
+           "BGR;16" or "BGR;24", then the default filter is
            :py:data:`Resampling.NEAREST`. Otherwise, the default filter is
            :py:data:`Resampling.BICUBIC`. See: :ref:`concept-filters`.
         :param box: An optional 4-tuple of floats providing
@@ -2302,8 +2302,8 @@ class Image:
         """
 
         if resample is None:
-            type_special = ";" in self.mode
-            resample = Resampling.NEAREST if type_special else Resampling.BICUBIC
+            bgr = self.mode.startswith("BGR;")
+            resample = Resampling.NEAREST if bgr else Resampling.BICUBIC
         elif resample not in (
             Resampling.NEAREST,
             Resampling.BILINEAR,
diff --git a/src/_imaging.c b/src/_imaging.c
index 07d9a64cc..426036335 100644
--- a/src/_imaging.c
+++ b/src/_imaging.c
@@ -1579,16 +1579,12 @@ _putdata(ImagingObject *self, PyObject *args) {
                 int bigendian = 0;
                 if (image->type == IMAGING_TYPE_SPECIAL) {
                     // I;16*
-                    if (strcmp(image->mode, "I;16N") == 0) {
+                    if (strcmp(image->mode, "I;16B") == 0
 #ifdef WORDS_BIGENDIAN
-                        bigendian = 1;
-#else
-                        bigendian = 0;
+                        || strcmp(image->mode, "I;16N") == 0
 #endif
-                    } else if (strcmp(image->mode, "I;16B") == 0) {
+                    ) {
                         bigendian = 1;
-                    } else {
-                        bigendian = 0;
                     }
                 }
                 for (i = x = y = 0; i < n; i++) {
diff --git a/src/libImaging/Resample.c b/src/libImaging/Resample.c
index 222d6bca4..f5e386dc2 100644
--- a/src/libImaging/Resample.c
+++ b/src/libImaging/Resample.c
@@ -460,6 +460,83 @@ ImagingResampleVertical_8bpc(
     ImagingSectionLeave(&cookie);
 }
 
+void
+ImagingResampleHorizontal_16bpc(
+    Imaging imOut, Imaging imIn, int offset, int ksize, int *bounds, double *kk
+) {
+    ImagingSectionCookie cookie;
+    double ss;
+    int xx, yy, x, xmin, xmax, ss_int;
+    double *k;
+
+    int bigendian = 0;
+    if (strcmp(imIn->mode, "I;16N") == 0
+#ifdef WORDS_BIGENDIAN
+        || strcmp(imIn->mode, "I;16B") == 0
+#endif
+    ) {
+        bigendian = 1;
+    }
+
+    ImagingSectionEnter(&cookie);
+    for (yy = 0; yy < imOut->ysize; yy++) {
+        for (xx = 0; xx < imOut->xsize; xx++) {
+            xmin = bounds[xx * 2 + 0];
+            xmax = bounds[xx * 2 + 1];
+            k = &kk[xx * ksize];
+            ss = 0.0;
+            for (x = 0; x < xmax; x++) {
+                ss += (imIn->image8[yy + offset][(x + xmin) * 2 + (bigendian ? 1 : 0)] +
+                       (imIn->image8[yy + offset][(x + xmin) * 2 + (bigendian ? 0 : 1)]
+                        << 8)) *
+                      k[x];
+            }
+            ss_int = ROUND_UP(ss);
+            imOut->image8[yy][xx * 2 + (bigendian ? 1 : 0)] = CLIP8(ss_int % 256);
+            imOut->image8[yy][xx * 2 + (bigendian ? 0 : 1)] = CLIP8(ss_int >> 8);
+        }
+    }
+    ImagingSectionLeave(&cookie);
+}
+
+void
+ImagingResampleVertical_16bpc(
+    Imaging imOut, Imaging imIn, int offset, int ksize, int *bounds, double *kk
+) {
+    ImagingSectionCookie cookie;
+    double ss;
+    int xx, yy, y, ymin, ymax, ss_int;
+    double *k;
+
+    int bigendian = 0;
+    if (strcmp(imIn->mode, "I;16N") == 0
+#ifdef WORDS_BIGENDIAN
+        || strcmp(imIn->mode, "I;16B") == 0
+#endif
+    ) {
+        bigendian = 1;
+    }
+
+    ImagingSectionEnter(&cookie);
+    for (yy = 0; yy < imOut->ysize; yy++) {
+        ymin = bounds[yy * 2 + 0];
+        ymax = bounds[yy * 2 + 1];
+        k = &kk[yy * ksize];
+        for (xx = 0; xx < imOut->xsize; xx++) {
+            ss = 0.0;
+            for (y = 0; y < ymax; y++) {
+                ss += (imIn->image8[y + ymin][xx * 2 + (bigendian ? 1 : 0)] +
+                       (imIn->image8[y + ymin][xx * 2 + (bigendian ? 0 : 1)] << 8)) *
+                      k[y];
+            }
+            ss_int = ROUND_UP(ss);
+            imOut->image8[yy][xx * 2 + (bigendian ? 1 : 0)] = CLIP8(ss_int % 256);
+            imOut->image8[yy][xx * 2 + (bigendian ? 0 : 1)] = CLIP8(ss_int >> 8);
+        }
+    }
+    ImagingSectionLeave(&cookie);
+}
+
 void
 ImagingResampleHorizontal_32bpc(
     Imaging imOut, Imaging imIn, int offset, int ksize, int *bounds, double *kk
@@ -574,7 +651,12 @@ ImagingResample(Imaging imIn, int xsize, int ysize, int filter, float box[4]) {
     }
 
     if (imIn->type == IMAGING_TYPE_SPECIAL) {
-        return (Imaging)ImagingError_ModeError();
+        if (strncmp(imIn->mode, "I;16", 4) == 0) {
+            ResampleHorizontal = ImagingResampleHorizontal_16bpc;
+            ResampleVertical = ImagingResampleVertical_16bpc;
+        } else {
+            return (Imaging)ImagingError_ModeError();
+        }
     } else if (imIn->image8) {
         ResampleHorizontal = ImagingResampleHorizontal_8bpc;
         ResampleVertical = ImagingResampleVertical_8bpc;