diff --git a/pandas-stubs/_typing.pyi b/pandas-stubs/_typing.pyi
index 1dc22c5d9..e2c8bedd7 100644
--- a/pandas-stubs/_typing.pyi
+++ b/pandas-stubs/_typing.pyi
@@ -15,6 +15,7 @@ from re import Pattern
 import sys
 from typing import (
     Any,
+    Generic,
     Literal,
     Protocol,
     SupportsIndex,
@@ -37,6 +38,7 @@ from typing_extensions import (
     ParamSpec,
     TypeAlias,
     TypeVar,
+    override,
 )
 
 from pandas._libs.interval import Interval
@@ -818,6 +820,8 @@ np_ndarray_float: TypeAlias = npt.NDArray[np.floating]
 np_ndarray_complex: TypeAlias = npt.NDArray[np.complexfloating]
 np_ndarray_bool: TypeAlias = npt.NDArray[np.bool_]
 np_ndarray_str: TypeAlias = npt.NDArray[np.str_]
+np_ndarray_dt: TypeAlias = npt.NDArray[np.datetime64]
+np_ndarray_td: TypeAlias = npt.NDArray[np.timedelta64]
 
 # Define shape and generic type variables with defaults similar to numpy
 GenericT = TypeVar("GenericT", bound=np.generic, default=Any)
@@ -1055,4 +1059,14 @@ DictConvertible: TypeAlias = FulldatetimeDict | DataFrame
 # where it is the only acceptable type.
 Incomplete: TypeAlias = Any
 
+# differentiating between bool and int/float/complex
+# https://github.com/pandas-dev/pandas-stubs/pull/1312#pullrequestreview-3126128971
+class Just(Protocol, Generic[T]):
+    @property  # type: ignore[override]
+    @override
+    def __class__(self, /) -> type[T]: ...
+    @__class__.setter
+    @override
+    def __class__(self, t: type[T], /) -> None: ...
+
 __all__ = ["npt", "type_t"]
diff --git a/pandas-stubs/core/groupby/groupby.pyi b/pandas-stubs/core/groupby/groupby.pyi
index 1490e7d17..9a620a450 100644
--- a/pandas-stubs/core/groupby/groupby.pyi
+++ b/pandas-stubs/core/groupby/groupby.pyi
@@ -74,9 +74,9 @@ from pandas._typing import (
 from pandas.plotting import PlotAccessor
 
 _ResamplerGroupBy: TypeAlias = (
-    DatetimeIndexResamplerGroupby[NDFrameT]
-    | PeriodIndexResamplerGroupby[NDFrameT]
-    | TimedeltaIndexResamplerGroupby[NDFrameT]
+    DatetimeIndexResamplerGroupby[NDFrameT]  # ty: ignore[invalid-argument-type]
+    | PeriodIndexResamplerGroupby[NDFrameT]  # ty: ignore[invalid-argument-type]
+    | TimedeltaIndexResamplerGroupby[NDFrameT]  # ty: ignore[invalid-argument-type]
 )
 
 class GroupBy(BaseGroupBy[NDFrameT]):
diff --git a/pandas-stubs/core/series.pyi b/pandas-stubs/core/series.pyi
index f0f59ffc0..eddb77f23 100644
--- a/pandas-stubs/core/series.pyi
+++ b/pandas-stubs/core/series.pyi
@@ -148,6 +148,7 @@ from pandas._typing import (
     JoinHow,
     JSONSerializable,
     JsonSeriesOrient,
+    Just,
     Label,
     Level,
     ListLike,
@@ -185,7 +186,9 @@ from pandas._typing import (
     np_ndarray_anyint,
     np_ndarray_bool,
     np_ndarray_complex,
+    np_ndarray_dt,
     np_ndarray_float,
+    np_ndarray_td,
     npt,
     num,
 )
@@ -2062,7 +2065,6 @@ class Series(IndexOpsMixin[S1], NDFrame):
     ) -> Series[bool]: ...
     @overload
     def __ror__(self, other: int | np_ndarray_anyint | Series[int]) -> Series[int]: ...
-    def __rsub__(self, other: num | _ListLike | Series[S1]) -> Series: ...
     # ignore needed for mypy as we want different results based on the arguments
     @overload  # type: ignore[override]
     def __rxor__(  # pyright: ignore[reportOverlappingOverload]
@@ -2071,21 +2073,417 @@ class Series(IndexOpsMixin[S1], NDFrame):
     @overload
     def __rxor__(self, other: int | np_ndarray_anyint | Series[int]) -> Series[int]: ...
     @overload
+    def __sub__(
+        self: Series[Never],
+        other: datetime | np.datetime64 | np_ndarray_dt | TimestampSeries,
+    ) -> TimedeltaSeries: ...
+    @overload
+    def __sub__(self: Series[Never], other: complex | _ListLike | Series) -> Series: ...
+    @overload
+    def __sub__(self, other: Series[Never]) -> Series: ...  # type: ignore[overload-overlap]
+    @overload
+    def __sub__(
+        self: Series[bool],
+        other: Just[int] | Sequence[Just[int]] | np_ndarray_anyint | Series[int],
+    ) -> Series[int]: ...
+    @overload
+    def __sub__(
+        self: Series[bool],
+        other: Just[float] | Sequence[Just[float]] | np_ndarray_float | Series[float],
+    ) -> Series[float]: ...
+    @overload
+    def __sub__(
+        self: Series[int],
+        other: (
+            int
+            | Sequence[int]
+            | np_ndarray_bool
+            | np_ndarray_anyint
+            | Series[bool]
+            | Series[int]
+        ),
+    ) -> Series[int]: ...
+    @overload
+    def __sub__(
+        self: Series[int],
+        other: Just[float] | Sequence[Just[float]] | np_ndarray_float | Series[float],
+    ) -> Series[float]: ...
+    @overload
+    def __sub__(
+        self: Series[float],
+        other: (
+            float
+            | Sequence[float]
+            | np_ndarray_bool
+            | np_ndarray_anyint
+            | np_ndarray_float
+            | Series[bool]
+            | Series[int]
+            | Series[float]
+        ),
+    ) -> Series[float]: ...
+    @overload
+    def __sub__(
+        self: Series[complex],
+        other: (
+            _T_COMPLEX
+            | Sequence[_T_COMPLEX]
+            | np_ndarray_bool
+            | np_ndarray_anyint
+            | np_ndarray_float
+            | Series[_T_COMPLEX]
+        ),
+    ) -> Series[complex]: ...
+    @overload
+    def __sub__(
+        self: Series[_T_COMPLEX],
+        other: (
+            Just[complex]
+            | Sequence[Just[complex]]
+            | np_ndarray_complex
+            | Series[complex]
+        ),
+    ) -> Series[complex]: ...
+    @overload
     def __sub__(
         self: Series[Timestamp],
-        other: Timedelta | TimedeltaSeries | TimedeltaIndex | np.timedelta64,
+        other: (
+            timedelta
+            | np.timedelta64
+            | np_ndarray_td
+            | TimedeltaSeries
+            | TimedeltaIndex
+        ),
     ) -> TimestampSeries: ...
     @overload
     def __sub__(
         self: Series[Timedelta],
-        other: Timedelta | TimedeltaSeries | TimedeltaIndex | np.timedelta64,
+        other: (
+            timedelta
+            | np.timedelta64
+            | np_ndarray_td
+            | TimedeltaSeries
+            | TimedeltaIndex
+        ),
     ) -> TimedeltaSeries: ...
     @overload
-    def __sub__(
-        self, other: Timestamp | datetime | TimestampSeries
+    def sub(
+        self: Series[Never],
+        other: datetime | np.datetime64 | np_ndarray_dt | TimestampSeries,
+        level: Level | None = None,
+        fill_value: float | None = None,
+        axis: int = 0,
+    ) -> TimedeltaSeries: ...
+    @overload
+    def sub(
+        self: Series[Never],
+        other: complex | _ListLike | Series,
+        level: Level | None = None,
+        fill_value: float | None = None,
+        axis: int = 0,
+    ) -> Series: ...
+    @overload
+    def sub(  # type: ignore[overload-overlap]
+        self,
+        other: Series[Never],
+        level: Level | None = None,
+        fill_value: float | None = None,
+        axis: int = 0,
+    ) -> Series: ...
+    @overload
+    def sub(
+        self: Series[bool],
+        other: Just[int] | Sequence[Just[int]] | np_ndarray_anyint | Series[int],
+        level: Level | None = None,
+        fill_value: float | None = None,
+        axis: int = 0,
+    ) -> Series[int]: ...
+    @overload
+    def sub(
+        self: Series[bool],
+        other: Just[float] | Sequence[Just[float]] | np_ndarray_float | Series[float],
+        level: Level | None = None,
+        fill_value: float | None = None,
+        axis: int = 0,
+    ) -> Series[float]: ...
+    @overload
+    def sub(
+        self: Series[int],
+        other: (
+            int
+            | Sequence[int]
+            | np_ndarray_bool
+            | np_ndarray_anyint
+            | Series[bool]
+            | Series[int]
+        ),
+        level: Level | None = None,
+        fill_value: float | None = None,
+        axis: int = 0,
+    ) -> Series[int]: ...
+    @overload
+    def sub(
+        self: Series[int],
+        other: Just[float] | Sequence[Just[float]] | np_ndarray_float | Series[float],
+        level: Level | None = None,
+        fill_value: float | None = None,
+        axis: int = 0,
+    ) -> Series[float]: ...
+    @overload
+    def sub(
+        self: Series[float],
+        other: (
+            float
+            | Sequence[float]
+            | np_ndarray_bool
+            | np_ndarray_anyint
+            | np_ndarray_float
+            | Series[bool]
+            | Series[int]
+            | Series[float]
+        ),
+        level: Level | None = None,
+        fill_value: float | None = None,
+        axis: int = 0,
+    ) -> Series[float]: ...
+    @overload
+    def sub(
+        self: Series[complex],
+        other: (
+            _T_COMPLEX
+            | Sequence[_T_COMPLEX]
+            | np_ndarray_bool
+            | np_ndarray_anyint
+            | np_ndarray_float
+            | Series[_T_COMPLEX]
+        ),
+        level: Level | None = None,
+        fill_value: float | None = None,
+        axis: int = 0,
+    ) -> Series[complex]: ...
+    @overload
+    def sub(
+        self: Series[_T_COMPLEX],
+        other: (
+            Just[complex]
+            | Sequence[Just[complex]]
+            | np_ndarray_complex
+            | Series[complex]
+        ),
+        level: Level | None = None,
+        fill_value: float | None = None,
+        axis: int = 0,
+    ) -> Series[complex]: ...
+    @overload
+    def sub(
+        self: Series[Timestamp],
+        other: (
+            timedelta
+            | np.timedelta64
+            | np_ndarray_td
+            | TimedeltaSeries
+            | TimedeltaIndex
+        ),
+        level: Level | None = None,
+        fill_value: float | None = None,
+        axis: int = 0,
+    ) -> TimestampSeries: ...
+    @overload
+    def sub(
+        self: Series[Timedelta],
+        other: (
+            timedelta
+            | np.timedelta64
+            | np_ndarray_td
+            | TimedeltaSeries
+            | TimedeltaIndex
+        ),
+        level: Level | None = None,
+        fill_value: float | None = None,
+        axis: int = 0,
+    ) -> TimedeltaSeries: ...
+    @overload
+    def __rsub__(  # type: ignore[misc]
+        self: Series[Never],
+        other: datetime | np.datetime64 | np_ndarray_dt | TimestampSeries,
+    ) -> TimedeltaSeries: ...
+    @overload
+    def __rsub__(
+        self: Series[Never], other: complex | _ListLike | Series
+    ) -> Series: ...
+    @overload
+    def __rsub__(self, other: Series[Never]) -> Series: ...
+    @overload
+    def __rsub__(
+        self: Series[bool],
+        other: Just[int] | Sequence[Just[int]] | np_ndarray_anyint | Series[int],
+    ) -> Series[int]: ...
+    @overload
+    def __rsub__(
+        self: Series[bool],
+        other: Just[float] | Sequence[Just[float]] | np_ndarray_float | Series[float],
+    ) -> Series[float]: ...
+    @overload
+    def __rsub__(
+        self: Series[int],
+        other: (
+            int
+            | Sequence[int]
+            | np_ndarray_bool
+            | np_ndarray_anyint
+            | Series[bool]
+            | Series[int]
+        ),
+    ) -> Series[int]: ...
+    @overload
+    def __rsub__(
+        self: Series[int],
+        other: Just[float] | Sequence[Just[float]] | np_ndarray_float | Series[float],
+    ) -> Series[float]: ...
+    @overload
+    def __rsub__(
+        self: Series[float],
+        other: (
+            float
+            | Sequence[float]
+            | np_ndarray_bool
+            | np_ndarray_anyint
+            | np_ndarray_float
+            | Series[bool]
+            | Series[int]
+            | Series[float]
+        ),
+    ) -> Series[float]: ...
+    @overload
+    def __rsub__(
+        self: Series[complex],
+        other: (
+            _T_COMPLEX
+            | Sequence[_T_COMPLEX]
+            | np_ndarray_bool
+            | np_ndarray_anyint
+            | np_ndarray_float
+            | Series[_T_COMPLEX]
+        ),
+    ) -> Series[complex]: ...
+    @overload
+    def __rsub__(
+        self: Series[_T_COMPLEX],
+        other: (
+            Just[complex]
+            | Sequence[Just[complex]]
+            | np_ndarray_complex
+            | Series[complex]
+        ),
+    ) -> Series[complex]: ...
+    @overload
+    def rsub(
+        self: Series[Never],
+        other: datetime | np.datetime64 | np_ndarray_dt | TimestampSeries,
+        level: Level | None = None,
+        fill_value: float | None = None,
+        axis: int = 0,
     ) -> TimedeltaSeries: ...
     @overload
-    def __sub__(self, other: num | _ListLike | Series) -> Series: ...
+    def rsub(
+        self: Series[Never],
+        other: complex | _ListLike | Series,
+        level: Level | None = None,
+        fill_value: float | None = None,
+        axis: int = 0,
+    ) -> Series: ...
+    @overload
+    def rsub(
+        self,
+        other: Series[Never],
+        level: Level | None = None,
+        fill_value: float | None = None,
+        axis: int = 0,
+    ) -> Series: ...
+    @overload
+    def rsub(
+        self: Series[bool],
+        other: Just[int] | Sequence[Just[int]] | np_ndarray_anyint | Series[int],
+        level: Level | None = None,
+        fill_value: float | None = None,
+        axis: int = 0,
+    ) -> Series[int]: ...
+    @overload
+    def rsub(
+        self: Series[bool],
+        other: Just[float] | Sequence[Just[float]] | np_ndarray_float | Series[float],
+        level: Level | None = None,
+        fill_value: float | None = None,
+        axis: int = 0,
+    ) -> Series[float]: ...
+    @overload
+    def rsub(
+        self: Series[int],
+        other: (
+            int
+            | Sequence[int]
+            | np_ndarray_bool
+            | np_ndarray_anyint
+            | Series[bool]
+            | Series[int]
+        ),
+        level: Level | None = None,
+        fill_value: float | None = None,
+        axis: int = 0,
+    ) -> Series[int]: ...
+    @overload
+    def rsub(
+        self: Series[int],
+        other: Just[float] | Sequence[Just[float]] | np_ndarray_float | Series[float],
+        level: Level | None = None,
+        fill_value: float | None = None,
+        axis: int = 0,
+    ) -> Series[float]: ...
+    @overload
+    def rsub(
+        self: Series[float],
+        other: (
+            float
+            | Sequence[float]
+            | np_ndarray_bool
+            | np_ndarray_anyint
+            | np_ndarray_float
+            | Series[bool]
+            | Series[int]
+            | Series[float]
+        ),
+        level: Level | None = None,
+        fill_value: float | None = None,
+        axis: int = 0,
+    ) -> Series[float]: ...
+    @overload
+    def rsub(
+        self: Series[complex],
+        other: (
+            _T_COMPLEX
+            | Sequence[_T_COMPLEX]
+            | np_ndarray_bool
+            | np_ndarray_anyint
+            | np_ndarray_float
+            | Series[_T_COMPLEX]
+        ),
+        level: Level | None = None,
+        fill_value: float | None = None,
+        axis: int = 0,
+    ) -> Series[complex]: ...
+    @overload
+    def rsub(
+        self: Series[_T_COMPLEX],
+        other: (
+            Just[complex]
+            | Sequence[Just[complex]]
+            | np_ndarray_complex
+            | Series[complex]
+        ),
+        level: Level | None = None,
+        fill_value: float | None = None,
+        axis: int = 0,
+    ) -> Series[complex]: ...
     @overload
     def __truediv__(
         self: Series[Never], other: complex | _ListLike | Series
@@ -2093,12 +2491,13 @@ class Series(IndexOpsMixin[S1], NDFrame):
     @overload
     def __truediv__(self, other: Series[Never]) -> Series: ...
     @overload
-    def __truediv__(self: Series[bool], other: bool | np_ndarray_bool) -> Never: ...
+    def __truediv__(self: Series[bool], other: np_ndarray_bool) -> Never: ...
     @overload
-    def __truediv__(  # pyright: ignore[reportOverlappingOverload]
+    def __truediv__(
         self: Series[bool],
         other: (
-            float
+            Just[int]
+            | Just[float]
             | Sequence[float]
             | np_ndarray_anyint
             | np_ndarray_float
@@ -2108,7 +2507,8 @@ class Series(IndexOpsMixin[S1], NDFrame):
     ) -> Series[float]: ...
     @overload
     def __truediv__(
-        self: Series[bool], other: complex | Sequence[complex] | Series[complex]
+        self: Series[bool],
+        other: Just[complex] | Sequence[Just[complex]] | Series[complex],
     ) -> Series[complex]: ...
     @overload
     def __truediv__(
@@ -2180,17 +2580,10 @@ class Series(IndexOpsMixin[S1], NDFrame):
     ) -> Series: ...
     @overload
     def truediv(
-        self: Series[bool],
-        other: bool | np_ndarray_bool,
-        level: Level | None = None,
-        fill_value: float | None = None,
-        axis: AxisIndex = 0,
-    ) -> Never: ...
-    @overload
-    def truediv(  # pyright: ignore[reportOverlappingOverload]
         self: Series[bool],
         other: (
-            float
+            Just[int]
+            | Just[float]
             | Sequence[float]
             | np_ndarray_anyint
             | np_ndarray_float
@@ -2204,7 +2597,7 @@ class Series(IndexOpsMixin[S1], NDFrame):
     @overload
     def truediv(
         self: Series[bool],
-        other: complex | Sequence[complex] | Series[complex],
+        other: Just[complex] | Sequence[Just[complex]] | Series[complex],
         level: Level | None = None,
         fill_value: float | None = None,
         axis: AxisIndex = 0,
@@ -2294,12 +2687,13 @@ class Series(IndexOpsMixin[S1], NDFrame):
     @overload
     def __rtruediv__(self, other: Series[Never]) -> Series: ...
     @overload
-    def __rtruediv__(self: Series[bool], other: bool | np_ndarray_bool) -> Never: ...
+    def __rtruediv__(self: Series[bool], other: np_ndarray_bool) -> Never: ...
     @overload
-    def __rtruediv__(  # pyright: ignore[reportOverlappingOverload]
+    def __rtruediv__(
         self: Series[bool],
         other: (
-            float
+            Just[int]
+            | Just[float]
             | Sequence[float]
             | np_ndarray_anyint
             | np_ndarray_float
@@ -2309,7 +2703,8 @@ class Series(IndexOpsMixin[S1], NDFrame):
     ) -> Series[float]: ...
     @overload
     def __rtruediv__(
-        self: Series[bool], other: complex | Sequence[complex] | Series[complex]
+        self: Series[bool],
+        other: Just[complex] | Sequence[Just[complex]] | Series[complex],
     ) -> Series[complex]: ...
     @overload
     def __rtruediv__(
@@ -2380,17 +2775,10 @@ class Series(IndexOpsMixin[S1], NDFrame):
     ) -> Series: ...
     @overload
     def rtruediv(
-        self: Series[bool],
-        other: bool | np_ndarray_bool,
-        level: Level | None = None,
-        fill_value: float | None = None,
-        axis: AxisIndex = 0,
-    ) -> Never: ...
-    @overload
-    def rtruediv(  # pyright: ignore[reportOverlappingOverload]
         self: Series[bool],
         other: (
-            float
+            Just[int]
+            | Just[float]
             | Sequence[float]
             | np_ndarray_anyint
             | np_ndarray_float
@@ -2404,7 +2792,7 @@ class Series(IndexOpsMixin[S1], NDFrame):
     @overload
     def rtruediv(
         self: Series[bool],
-        other: complex | Sequence[complex] | Series[complex],
+        other: Just[complex] | Sequence[Just[complex]] | Series[complex],
         level: Level | None = None,
         fill_value: float | None = None,
         axis: AxisIndex = 0,
@@ -2811,13 +3199,6 @@ class Series(IndexOpsMixin[S1], NDFrame):
         fill_value: float | None = None,
         axis: AxisIndex = ...,
     ) -> Series[S1]: ...
-    def rsub(
-        self,
-        other: Series[S1] | Scalar,
-        level: Level | None = ...,
-        fill_value: float | None = None,
-        axis: AxisIndex = ...,
-    ) -> Series[S1]: ...
     def sem(
         self,
         axis: AxisIndex | None = 0,
@@ -2841,20 +3222,6 @@ class Series(IndexOpsMixin[S1], NDFrame):
         numeric_only: _bool = False,
         **kwargs: Any,
     ) -> float: ...
-    def sub(
-        self,
-        other: num | _ListLike | Series[S1],
-        level: Level | None = ...,
-        fill_value: float | None = None,
-        axis: AxisIndex | None = 0,
-    ) -> Series[S1]: ...
-    def subtract(
-        self,
-        other: num | _ListLike | Series[S1],
-        level: Level | None = ...,
-        fill_value: float | None = None,
-        axis: AxisIndex | None = 0,
-    ) -> Series[S1]: ...
     @overload
     def sum(
         self: Series[Never],
diff --git a/tests/series/arithmetic/bool/test_sub.py b/tests/series/arithmetic/bool/test_sub.py
new file mode 100644
index 000000000..958cb8e3b
--- /dev/null
+++ b/tests/series/arithmetic/bool/test_sub.py
@@ -0,0 +1,152 @@
+import numpy as np
+from numpy import typing as npt  # noqa: F401
+import pandas as pd
+from typing_extensions import (
+    Never,
+    assert_type,
+)
+
+from tests import (
+    TYPE_CHECKING_INVALID_USAGE,
+    check,
+)
+
+left = pd.Series([True, True, False])  # left operand
+
+
+def test_sub_py_scalar() -> None:
+    """Test pd.Series[bool] - Python native scalars"""
+    b, i, f, c = True, 1, 1.0, 1j
+
+    if TYPE_CHECKING_INVALID_USAGE:
+        _0 = left - b  # type: ignore[operator] # pyright: ignore[reportOperatorIssue]
+    check(assert_type(left - i, "pd.Series[int]"), pd.Series, np.integer)
+    check(assert_type(left - f, "pd.Series[float]"), pd.Series, np.floating)
+    check(assert_type(left - c, "pd.Series[complex]"), pd.Series, np.complexfloating)
+
+    if TYPE_CHECKING_INVALID_USAGE:
+        _1 = b - left  # type: ignore[operator] # pyright: ignore[reportOperatorIssue]
+    check(assert_type(i - left, "pd.Series[int]"), pd.Series, np.integer)
+    check(assert_type(f - left, "pd.Series[float]"), pd.Series, np.floating)
+    check(assert_type(c - left, "pd.Series[complex]"), pd.Series, np.complexfloating)
+
+    if TYPE_CHECKING_INVALID_USAGE:
+        left.sub(b)  # type: ignore[arg-type] # pyright: ignore[reportArgumentType,reportCallIssue]
+    check(assert_type(left.sub(i), "pd.Series[int]"), pd.Series, np.integer)
+    check(assert_type(left.sub(f), "pd.Series[float]"), pd.Series, np.floating)
+    check(assert_type(left.sub(c), "pd.Series[complex]"), pd.Series, np.complexfloating)
+
+    if TYPE_CHECKING_INVALID_USAGE:
+        left.rsub(b)  # type: ignore[arg-type] # pyright: ignore[reportArgumentType,reportCallIssue]
+    check(assert_type(left.rsub(i), "pd.Series[int]"), pd.Series, np.integer)
+    check(assert_type(left.rsub(f), "pd.Series[float]"), pd.Series, np.floating)
+    check(
+        assert_type(left.rsub(c), "pd.Series[complex]"), pd.Series, np.complexfloating
+    )
+
+
+def test_sub_py_sequence() -> None:
+    """Test pd.Series[bool] - Python native sequence"""
+    b, i, f, c = [True, False, True], [2, 3, 5], [1.0, 2.0, 3.0], [1j, 1j, 4j]
+
+    if TYPE_CHECKING_INVALID_USAGE:
+        _0 = left - b  # type: ignore[operator] # pyright: ignore[reportOperatorIssue]
+    check(assert_type(left - i, "pd.Series[int]"), pd.Series, np.integer)
+    check(assert_type(left - f, "pd.Series[float]"), pd.Series, np.floating)
+    check(assert_type(left - c, "pd.Series[complex]"), pd.Series, np.complexfloating)
+
+    if TYPE_CHECKING_INVALID_USAGE:
+        _1 = b - left  # type: ignore[operator] # pyright: ignore[reportOperatorIssue]
+    check(assert_type(i - left, "pd.Series[int]"), pd.Series, np.integer)
+    check(assert_type(f - left, "pd.Series[float]"), pd.Series, np.floating)
+    check(assert_type(c - left, "pd.Series[complex]"), pd.Series, np.complexfloating)
+
+    if TYPE_CHECKING_INVALID_USAGE:
+        left.sub(b)  # type: ignore[arg-type] # pyright: ignore[reportArgumentType,reportCallIssue]
+    check(assert_type(left.sub(i), "pd.Series[int]"), pd.Series, np.integer)
+    check(assert_type(left.sub(f), "pd.Series[float]"), pd.Series, np.floating)
+    check(assert_type(left.sub(c), "pd.Series[complex]"), pd.Series, np.complexfloating)
+
+    if TYPE_CHECKING_INVALID_USAGE:
+        left.rsub(b)  # type: ignore[arg-type] # pyright: ignore[reportArgumentType,reportCallIssue]
+    check(assert_type(left.rsub(i), "pd.Series[int]"), pd.Series, np.integer)
+    check(assert_type(left.rsub(f), "pd.Series[float]"), pd.Series, np.floating)
+    check(
+        assert_type(left.rsub(c), "pd.Series[complex]"), pd.Series, np.complexfloating
+    )
+
+
+def test_sub_numpy_array() -> None:
+    """Test pd.Series[bool] - numpy array"""
+    b = np.array([True, False, True], np.bool_)
+    i = np.array([2, 3, 5], np.int64)
+    f = np.array([1.0, 2.0, 3.0], np.float64)
+    c = np.array([1.1j, 2.2j, 4.1j], np.complex128)
+
+    if TYPE_CHECKING_INVALID_USAGE:
+        assert_type(left - b, Never)
+    check(assert_type(left - i, "pd.Series[int]"), pd.Series, np.integer)
+    check(assert_type(left - f, "pd.Series[float]"), pd.Series, np.floating)
+    check(assert_type(left - c, "pd.Series[complex]"), pd.Series, np.complexfloating)
+
+    # `numpy` typing gives the corresponding `ndarray`s in the static type
+    # checking, where our `__rsub__` cannot override. At runtime, they return
+    # `Series`s with the correct element type.
+    if TYPE_CHECKING_INVALID_USAGE:
+        assert_type(b - left, Never)
+    check(assert_type(i - left, "npt.NDArray[np.int64]"), pd.Series, np.integer)
+    check(assert_type(f - left, "npt.NDArray[np.float64]"), pd.Series, np.floating)
+    check(
+        assert_type(c - left, "npt.NDArray[np.complex128]"),
+        pd.Series,
+        np.complexfloating,
+    )
+
+    if TYPE_CHECKING_INVALID_USAGE:
+        left.sub(b)  # type: ignore[arg-type] # pyright: ignore[reportArgumentType,reportCallIssue]
+    check(assert_type(left.sub(i), "pd.Series[int]"), pd.Series, np.integer)
+    check(assert_type(left.sub(f), "pd.Series[float]"), pd.Series, np.floating)
+    check(assert_type(left.sub(c), "pd.Series[complex]"), pd.Series, np.complexfloating)
+
+    if TYPE_CHECKING_INVALID_USAGE:
+        left.rsub(b)  # type: ignore[arg-type] # pyright: ignore[reportArgumentType,reportCallIssue]
+    check(assert_type(left.rsub(i), "pd.Series[int]"), pd.Series, np.integer)
+    check(assert_type(left.rsub(f), "pd.Series[float]"), pd.Series, np.floating)
+    check(
+        assert_type(left.rsub(c), "pd.Series[complex]"), pd.Series, np.complexfloating
+    )
+
+
+def test_sub_pd_series() -> None:
+    """Test pd.Series[bool] - pandas series"""
+    b = pd.Series([True, False, True])
+    i = pd.Series([2, 3, 5])
+    f = pd.Series([1.0, 2.0, 3.0])
+    c = pd.Series([1.1j, 2.2j, 4.1j])
+
+    # In the following two cases, mypy fails to recognise the second operand as pd.Series[bool]
+    if TYPE_CHECKING_INVALID_USAGE:
+        _0 = left - b  # type: ignore[operator] # pyright: ignore[reportOperatorIssue]
+    check(assert_type(left - i, "pd.Series[int]"), pd.Series, np.integer)
+    check(assert_type(left - f, "pd.Series[float]"), pd.Series, np.floating)
+    check(assert_type(left - c, "pd.Series[complex]"), pd.Series, np.complexfloating)
+
+    if TYPE_CHECKING_INVALID_USAGE:
+        _1 = b - left  # type: ignore[operator] # pyright: ignore[reportOperatorIssue]
+    check(assert_type(i - left, "pd.Series[int]"), pd.Series, np.integer)
+    check(assert_type(f - left, "pd.Series[float]"), pd.Series, np.floating)
+    check(assert_type(c - left, "pd.Series[complex]"), pd.Series, np.complexfloating)
+
+    if TYPE_CHECKING_INVALID_USAGE:
+        left.sub(b)  # type: ignore[arg-type] # pyright: ignore[reportArgumentType,reportCallIssue]
+    check(assert_type(left.sub(i), "pd.Series[int]"), pd.Series, np.integer)
+    check(assert_type(left.sub(f), "pd.Series[float]"), pd.Series, np.floating)
+    check(assert_type(left.sub(c), "pd.Series[complex]"), pd.Series, np.complexfloating)
+
+    if TYPE_CHECKING_INVALID_USAGE:
+        left.rsub(b)  # type: ignore[arg-type] # pyright: ignore[reportArgumentType,reportCallIssue]
+    check(assert_type(left.rsub(i), "pd.Series[int]"), pd.Series, np.integer)
+    check(assert_type(left.rsub(f), "pd.Series[float]"), pd.Series, np.floating)
+    check(
+        assert_type(left.rsub(c), "pd.Series[complex]"), pd.Series, np.complexfloating
+    )
diff --git a/tests/series/arithmetic/bool/test_truediv.py b/tests/series/arithmetic/bool/test_truediv.py
index caf61d0e3..a96e712b8 100644
--- a/tests/series/arithmetic/bool/test_truediv.py
+++ b/tests/series/arithmetic/bool/test_truediv.py
@@ -19,19 +19,19 @@ def test_truediv_py_scalar() -> None:
     b, i, f, c = True, 1, 1.0, 1j
 
     if TYPE_CHECKING_INVALID_USAGE:
-        assert_type(left / b, Never)
+        _0 = left / b  # type: ignore[operator] # pyright: ignore[reportOperatorIssue]
     check(assert_type(left / i, "pd.Series[float]"), pd.Series, np.floating)
     check(assert_type(left / f, "pd.Series[float]"), pd.Series, np.floating)
     check(assert_type(left / c, "pd.Series[complex]"), pd.Series, np.complexfloating)
 
     if TYPE_CHECKING_INVALID_USAGE:
-        assert_type(b / left, Never)
+        _1 = b / left  # type: ignore[operator] # pyright: ignore[reportOperatorIssue]
     check(assert_type(i / left, "pd.Series[float]"), pd.Series, np.floating)
     check(assert_type(f / left, "pd.Series[float]"), pd.Series, np.floating)
     check(assert_type(c / left, "pd.Series[complex]"), pd.Series, np.complexfloating)
 
     if TYPE_CHECKING_INVALID_USAGE:
-        assert_type(left.truediv(b), Never)
+        left.truediv(b)  # type: ignore[arg-type] # pyright: ignore[reportArgumentType,reportCallIssue]
     check(assert_type(left.truediv(i), "pd.Series[float]"), pd.Series, np.floating)
     check(assert_type(left.truediv(f), "pd.Series[float]"), pd.Series, np.floating)
     check(
@@ -41,13 +41,13 @@ def test_truediv_py_scalar() -> None:
     )
 
     if TYPE_CHECKING_INVALID_USAGE:
-        assert_type(left.div(b), Never)
+        left.div(b)  # type: ignore[arg-type] # pyright: ignore[reportArgumentType,reportCallIssue]
     check(assert_type(left.div(i), "pd.Series[float]"), pd.Series, np.floating)
     check(assert_type(left.div(f), "pd.Series[float]"), pd.Series, np.floating)
     check(assert_type(left.div(c), "pd.Series[complex]"), pd.Series, np.complexfloating)
 
     if TYPE_CHECKING_INVALID_USAGE:
-        assert_type(left.rtruediv(b), Never)
+        left.rtruediv(b)  # type: ignore[arg-type] # pyright: ignore[reportArgumentType,reportCallIssue]
     check(assert_type(left.rtruediv(i), "pd.Series[float]"), pd.Series, np.floating)
     check(assert_type(left.rtruediv(f), "pd.Series[float]"), pd.Series, np.floating)
     check(
@@ -57,7 +57,7 @@ def test_truediv_py_scalar() -> None:
     )
 
     if TYPE_CHECKING_INVALID_USAGE:
-        assert_type(left.rdiv(b), Never)
+        left.rdiv(b)  # type: ignore[arg-type] # pyright: ignore[reportArgumentType,reportCallIssue]
     check(assert_type(left.rdiv(i), "pd.Series[float]"), pd.Series, np.floating)
     check(assert_type(left.rdiv(f), "pd.Series[float]"), pd.Series, np.floating)
     check(
@@ -137,7 +137,7 @@ def test_truediv_numpy_array() -> None:
     )
 
     if TYPE_CHECKING_INVALID_USAGE:
-        assert_type(left.truediv(b), Never)
+        left.truediv(b)  # type: ignore[arg-type]  # pyright: ignore[reportArgumentType,reportCallIssue]
     check(assert_type(left.truediv(i), "pd.Series[float]"), pd.Series, np.floating)
     check(assert_type(left.truediv(f), "pd.Series[float]"), pd.Series, np.floating)
     check(
@@ -147,13 +147,13 @@ def test_truediv_numpy_array() -> None:
     )
 
     if TYPE_CHECKING_INVALID_USAGE:
-        assert_type(left.div(b), Never)
+        left.div(b)  # type: ignore[arg-type]  # pyright: ignore[reportArgumentType,reportCallIssue]
     check(assert_type(left.div(i), "pd.Series[float]"), pd.Series, np.floating)
     check(assert_type(left.div(f), "pd.Series[float]"), pd.Series, np.floating)
     check(assert_type(left.div(c), "pd.Series[complex]"), pd.Series, np.complexfloating)
 
     if TYPE_CHECKING_INVALID_USAGE:
-        assert_type(left.rtruediv(b), Never)
+        left.rtruediv(b)  # type: ignore[arg-type]  # pyright: ignore[reportArgumentType,reportCallIssue]
     check(assert_type(left.rtruediv(i), "pd.Series[float]"), pd.Series, np.floating)
     check(assert_type(left.rtruediv(f), "pd.Series[float]"), pd.Series, np.floating)
     check(
@@ -163,7 +163,7 @@ def test_truediv_numpy_array() -> None:
     )
 
     if TYPE_CHECKING_INVALID_USAGE:
-        assert_type(left.rdiv(b), Never)
+        left.rdiv(b)  # type: ignore[arg-type]  # pyright: ignore[reportArgumentType,reportCallIssue]
     check(assert_type(left.rdiv(i), "pd.Series[float]"), pd.Series, np.floating)
     check(assert_type(left.rdiv(f), "pd.Series[float]"), pd.Series, np.floating)
     check(
@@ -179,13 +179,13 @@ def test_truediv_pd_series() -> None:
     c = pd.Series([1.1j, 2.2j, 4.1j])
 
     if TYPE_CHECKING_INVALID_USAGE:
-        _ = left / b  # type: ignore[operator] # pyright: ignore[reportOperatorIssue]
+        _0 = left / b  # type: ignore[operator] # pyright: ignore[reportOperatorIssue]
     check(assert_type(left / i, "pd.Series[float]"), pd.Series, np.floating)
     check(assert_type(left / f, "pd.Series[float]"), pd.Series, np.floating)
     check(assert_type(left / c, "pd.Series[complex]"), pd.Series, np.complexfloating)
 
     if TYPE_CHECKING_INVALID_USAGE:
-        _ = b / left  # type: ignore[operator] # pyright: ignore[reportOperatorIssue]
+        _1 = b / left  # type: ignore[operator] # pyright: ignore[reportOperatorIssue]
     check(assert_type(i / left, "pd.Series[float]"), pd.Series, np.floating)
     check(assert_type(f / left, "pd.Series[float]"), pd.Series, np.floating)
     check(assert_type(c / left, "pd.Series[complex]"), pd.Series, np.complexfloating)
diff --git a/tests/series/arithmetic/complex/test_sub.py b/tests/series/arithmetic/complex/test_sub.py
new file mode 100644
index 000000000..5a7cc0de6
--- /dev/null
+++ b/tests/series/arithmetic/complex/test_sub.py
@@ -0,0 +1,157 @@
+from typing import NoReturn
+
+import numpy as np
+from numpy import typing as npt  # noqa: F401
+import pandas as pd
+from typing_extensions import assert_type
+
+from tests import check
+
+left = pd.Series([1j, 2j, 3j])  # left operand
+
+
+def test_sub_py_scalar() -> None:
+    """Test pd.Series[complex] - Python native scalars"""
+    b, i, f, c = True, 1, 1.0, 1j
+
+    check(assert_type(left - b, "pd.Series[complex]"), pd.Series, np.complexfloating)
+    check(assert_type(left - i, "pd.Series[complex]"), pd.Series, np.complexfloating)
+    check(assert_type(left - f, "pd.Series[complex]"), pd.Series, np.complexfloating)
+    check(assert_type(left - c, "pd.Series[complex]"), pd.Series, np.complexfloating)
+
+    check(assert_type(b - left, "pd.Series[complex]"), pd.Series, np.complexfloating)
+    check(assert_type(i - left, "pd.Series[complex]"), pd.Series, np.complexfloating)
+    check(assert_type(f - left, "pd.Series[complex]"), pd.Series, np.complexfloating)
+    check(assert_type(c - left, "pd.Series[complex]"), pd.Series, np.complexfloating)
+
+    check(assert_type(left.sub(b), "pd.Series[complex]"), pd.Series, np.complexfloating)
+    check(assert_type(left.sub(i), "pd.Series[complex]"), pd.Series, np.complexfloating)
+    check(assert_type(left.sub(f), "pd.Series[complex]"), pd.Series, np.complexfloating)
+    check(assert_type(left.sub(c), "pd.Series[complex]"), pd.Series, np.complexfloating)
+
+    check(
+        assert_type(left.rsub(b), "pd.Series[complex]"), pd.Series, np.complexfloating
+    )
+    check(
+        assert_type(left.rsub(i), "pd.Series[complex]"), pd.Series, np.complexfloating
+    )
+    check(
+        assert_type(left.rsub(f), "pd.Series[complex]"), pd.Series, np.complexfloating
+    )
+    check(
+        assert_type(left.rsub(c), "pd.Series[complex]"), pd.Series, np.complexfloating
+    )
+
+
+def test_sub_py_sequence() -> None:
+    """Test pd.Series[complex] - Python native sequence"""
+    b, i, f, c = [True, False, True], [2, 3, 5], [1.0, 2.0, 3.0], [1j, 1j, 4j]
+
+    check(assert_type(left - b, "pd.Series[complex]"), pd.Series, np.complexfloating)
+    check(assert_type(left - i, "pd.Series[complex]"), pd.Series, np.complexfloating)
+    check(assert_type(left - f, "pd.Series[complex]"), pd.Series, np.complexfloating)
+    check(assert_type(left - c, "pd.Series[complex]"), pd.Series, np.complexfloating)
+
+    check(assert_type(b - left, "pd.Series[complex]"), pd.Series, np.complexfloating)
+    check(assert_type(i - left, "pd.Series[complex]"), pd.Series, np.complexfloating)
+    check(assert_type(f - left, "pd.Series[complex]"), pd.Series, np.complexfloating)
+    check(assert_type(c - left, "pd.Series[complex]"), pd.Series, np.complexfloating)
+
+    check(assert_type(left.sub(b), "pd.Series[complex]"), pd.Series, np.complexfloating)
+    check(assert_type(left.sub(i), "pd.Series[complex]"), pd.Series, np.complexfloating)
+    check(assert_type(left.sub(f), "pd.Series[complex]"), pd.Series, np.complexfloating)
+    check(assert_type(left.sub(c), "pd.Series[complex]"), pd.Series, np.complexfloating)
+
+    check(
+        assert_type(left.rsub(b), "pd.Series[complex]"), pd.Series, np.complexfloating
+    )
+    check(
+        assert_type(left.rsub(i), "pd.Series[complex]"), pd.Series, np.complexfloating
+    )
+    check(
+        assert_type(left.rsub(f), "pd.Series[complex]"), pd.Series, np.complexfloating
+    )
+    check(
+        assert_type(left.rsub(c), "pd.Series[complex]"), pd.Series, np.complexfloating
+    )
+
+
+def test_sub_numpy_array() -> None:
+    """Test pd.Series[complex] - numpy array"""
+    b = np.array([True, False, True], np.bool_)
+    i = np.array([2, 3, 5], np.int64)
+    f = np.array([1.0, 2.0, 3.0], np.float64)
+    c = np.array([1.1j, 2.2j, 4.1j], np.complex128)
+
+    check(assert_type(left - b, "pd.Series[complex]"), pd.Series, np.complexfloating)
+    check(assert_type(left - i, "pd.Series[complex]"), pd.Series, np.complexfloating)
+    check(assert_type(left - f, "pd.Series[complex]"), pd.Series, np.complexfloating)
+    check(assert_type(left - c, "pd.Series[complex]"), pd.Series, np.complexfloating)
+
+    # `numpy` typing gives the corresponding `ndarray`s in the static type
+    # checking, where our `__rsub__` cannot override. At runtime, they return
+    # `Series`s with the correct element type.
+    check(assert_type(b - left, NoReturn), pd.Series, np.complexfloating)
+    check(assert_type(i - left, "npt.NDArray[np.int64]"), pd.Series, np.complexfloating)
+    check(
+        assert_type(f - left, "npt.NDArray[np.float64]"), pd.Series, np.complexfloating
+    )
+    check(
+        assert_type(c - left, "npt.NDArray[np.complex128]"),
+        pd.Series,
+        np.complexfloating,
+    )
+
+    check(assert_type(left.sub(b), "pd.Series[complex]"), pd.Series, np.complexfloating)
+    check(assert_type(left.sub(i), "pd.Series[complex]"), pd.Series, np.complexfloating)
+    check(assert_type(left.sub(f), "pd.Series[complex]"), pd.Series, np.complexfloating)
+    check(assert_type(left.sub(c), "pd.Series[complex]"), pd.Series, np.complexfloating)
+
+    check(
+        assert_type(left.rsub(b), "pd.Series[complex]"), pd.Series, np.complexfloating
+    )
+    check(
+        assert_type(left.rsub(i), "pd.Series[complex]"), pd.Series, np.complexfloating
+    )
+    check(
+        assert_type(left.rsub(f), "pd.Series[complex]"), pd.Series, np.complexfloating
+    )
+    check(
+        assert_type(left.rsub(c), "pd.Series[complex]"), pd.Series, np.complexfloating
+    )
+
+
+def test_sub_pd_series() -> None:
+    """Test pd.Series[complex] - pandas series"""
+    b = pd.Series([True, False, True])
+    i = pd.Series([2, 3, 5])
+    f = pd.Series([1.0, 2.0, 3.0])
+    c = pd.Series([1.1j, 2.2j, 4.1j])
+
+    check(assert_type(left - b, "pd.Series[complex]"), pd.Series, np.complexfloating)
+    check(assert_type(left - i, "pd.Series[complex]"), pd.Series, np.complexfloating)
+    check(assert_type(left - f, "pd.Series[complex]"), pd.Series, np.complexfloating)
+    check(assert_type(left - c, "pd.Series[complex]"), pd.Series, np.complexfloating)
+
+    check(assert_type(b - left, "pd.Series[complex]"), pd.Series, np.complexfloating)
+    check(assert_type(i - left, "pd.Series[complex]"), pd.Series, np.complexfloating)
+    check(assert_type(f - left, "pd.Series[complex]"), pd.Series, np.complexfloating)
+    check(assert_type(c - left, "pd.Series[complex]"), pd.Series, np.complexfloating)
+
+    check(assert_type(left.sub(b), "pd.Series[complex]"), pd.Series, np.complexfloating)
+    check(assert_type(left.sub(i), "pd.Series[complex]"), pd.Series, np.complexfloating)
+    check(assert_type(left.sub(f), "pd.Series[complex]"), pd.Series, np.complexfloating)
+    check(assert_type(left.sub(c), "pd.Series[complex]"), pd.Series, np.complexfloating)
+
+    check(
+        assert_type(left.rsub(b), "pd.Series[complex]"), pd.Series, np.complexfloating
+    )
+    check(
+        assert_type(left.rsub(i), "pd.Series[complex]"), pd.Series, np.complexfloating
+    )
+    check(
+        assert_type(left.rsub(f), "pd.Series[complex]"), pd.Series, np.complexfloating
+    )
+    check(
+        assert_type(left.rsub(c), "pd.Series[complex]"), pd.Series, np.complexfloating
+    )
diff --git a/tests/series/arithmetic/float/test_sub.py b/tests/series/arithmetic/float/test_sub.py
new file mode 100644
index 000000000..c7095765e
--- /dev/null
+++ b/tests/series/arithmetic/float/test_sub.py
@@ -0,0 +1,131 @@
+from typing import NoReturn
+
+import numpy as np
+from numpy import typing as npt  # noqa: F401
+import pandas as pd
+from typing_extensions import assert_type
+
+from tests import check
+
+left = pd.Series([1.0, 2.0, 3.0])  # left operand
+
+
+def test_sub_py_scalar() -> None:
+    """Test pd.Series[float] - Python native scalars"""
+    b, i, f, c = True, 1, 1.0, 1j
+
+    check(assert_type(left - b, "pd.Series[float]"), pd.Series, np.floating)
+    check(assert_type(left - i, "pd.Series[float]"), pd.Series, np.floating)
+    check(assert_type(left - f, "pd.Series[float]"), pd.Series, np.floating)
+    check(assert_type(left - c, "pd.Series[complex]"), pd.Series, np.complexfloating)
+
+    check(assert_type(b - left, "pd.Series[float]"), pd.Series, np.floating)
+    check(assert_type(i - left, "pd.Series[float]"), pd.Series, np.floating)
+    check(assert_type(f - left, "pd.Series[float]"), pd.Series, np.floating)
+    check(assert_type(c - left, "pd.Series[complex]"), pd.Series, np.complexfloating)
+
+    check(assert_type(left.sub(b), "pd.Series[float]"), pd.Series, np.floating)
+    check(assert_type(left.sub(i), "pd.Series[float]"), pd.Series, np.floating)
+    check(assert_type(left.sub(f), "pd.Series[float]"), pd.Series, np.floating)
+    check(assert_type(left.sub(c), "pd.Series[complex]"), pd.Series, np.complexfloating)
+
+    check(assert_type(left.rsub(b), "pd.Series[float]"), pd.Series, np.floating)
+    check(assert_type(left.rsub(i), "pd.Series[float]"), pd.Series, np.floating)
+    check(assert_type(left.rsub(f), "pd.Series[float]"), pd.Series, np.floating)
+    check(
+        assert_type(left.rsub(c), "pd.Series[complex]"), pd.Series, np.complexfloating
+    )
+
+
+def test_sub_py_sequence() -> None:
+    """Test pd.Series[float] - Python native sequence"""
+    b, i, f, c = [True, False, True], [2, 3, 5], [1.0, 2.0, 3.0], [1j, 1j, 4j]
+
+    check(assert_type(left - b, "pd.Series[float]"), pd.Series, np.floating)
+    check(assert_type(left - i, "pd.Series[float]"), pd.Series, np.floating)
+    check(assert_type(left - f, "pd.Series[float]"), pd.Series, np.floating)
+    check(assert_type(left - c, "pd.Series[complex]"), pd.Series, np.complexfloating)
+
+    check(assert_type(b - left, "pd.Series[float]"), pd.Series, np.floating)
+    check(assert_type(i - left, "pd.Series[float]"), pd.Series, np.floating)
+    check(assert_type(f - left, "pd.Series[float]"), pd.Series, np.floating)
+    check(assert_type(c - left, "pd.Series[complex]"), pd.Series, np.complexfloating)
+
+    check(assert_type(left.sub(b), "pd.Series[float]"), pd.Series, np.floating)
+    check(assert_type(left.sub(i), "pd.Series[float]"), pd.Series, np.floating)
+    check(assert_type(left.sub(f), "pd.Series[float]"), pd.Series, np.floating)
+    check(assert_type(left.sub(c), "pd.Series[complex]"), pd.Series, np.complexfloating)
+
+    check(assert_type(left.rsub(b), "pd.Series[float]"), pd.Series, np.floating)
+    check(assert_type(left.rsub(i), "pd.Series[float]"), pd.Series, np.floating)
+    check(assert_type(left.rsub(f), "pd.Series[float]"), pd.Series, np.floating)
+    check(
+        assert_type(left.rsub(c), "pd.Series[complex]"), pd.Series, np.complexfloating
+    )
+
+
+def test_sub_numpy_array() -> None:
+    """Test pd.Series[float] - numpy array"""
+    b = np.array([True, False, True], np.bool_)
+    i = np.array([2, 3, 5], np.int64)
+    f = np.array([1.0, 2.0, 3.0], np.float64)
+    c = np.array([1.1j, 2.2j, 4.1j], np.complex128)
+
+    check(assert_type(left - b, "pd.Series[float]"), pd.Series, np.floating)
+    check(assert_type(left - i, "pd.Series[float]"), pd.Series, np.floating)
+    check(assert_type(left - f, "pd.Series[float]"), pd.Series, np.floating)
+    check(assert_type(left - c, "pd.Series[complex]"), pd.Series, np.complexfloating)
+
+    # `numpy` typing gives the corresponding `ndarray`s in the static type
+    # checking, where our `__rsub__` cannot override. At runtime, they return
+    # `Series`s with the correct element type.
+    check(assert_type(b - left, NoReturn), pd.Series, np.floating)
+    check(assert_type(i - left, "npt.NDArray[np.int64]"), pd.Series, np.floating)
+    check(assert_type(f - left, "npt.NDArray[np.float64]"), pd.Series, np.floating)
+    check(
+        assert_type(c - left, "npt.NDArray[np.complex128]"),
+        pd.Series,
+        np.complexfloating,
+    )
+
+    check(assert_type(left.sub(b), "pd.Series[float]"), pd.Series, np.floating)
+    check(assert_type(left.sub(i), "pd.Series[float]"), pd.Series, np.floating)
+    check(assert_type(left.sub(f), "pd.Series[float]"), pd.Series, np.floating)
+    check(assert_type(left.sub(c), "pd.Series[complex]"), pd.Series, np.complexfloating)
+
+    check(assert_type(left.rsub(b), "pd.Series[float]"), pd.Series, np.floating)
+    check(assert_type(left.rsub(i), "pd.Series[float]"), pd.Series, np.floating)
+    check(assert_type(left.rsub(f), "pd.Series[float]"), pd.Series, np.floating)
+    check(
+        assert_type(left.rsub(c), "pd.Series[complex]"), pd.Series, np.complexfloating
+    )
+
+
+def test_sub_pd_series() -> None:
+    """Test pd.Series[float] - pandas series"""
+    b = pd.Series([True, False, True])
+    i = pd.Series([2, 3, 5])
+    f = pd.Series([1.0, 2.0, 3.0])
+    c = pd.Series([1.1j, 2.2j, 4.1j])
+
+    check(assert_type(left - b, "pd.Series[float]"), pd.Series, np.floating)
+    check(assert_type(left - i, "pd.Series[float]"), pd.Series, np.floating)
+    check(assert_type(left - f, "pd.Series[float]"), pd.Series, np.floating)
+    check(assert_type(left - c, "pd.Series[complex]"), pd.Series, np.complexfloating)
+
+    check(assert_type(b - left, "pd.Series[float]"), pd.Series, np.floating)
+    check(assert_type(i - left, "pd.Series[float]"), pd.Series, np.floating)
+    check(assert_type(f - left, "pd.Series[float]"), pd.Series, np.floating)
+    check(assert_type(c - left, "pd.Series[complex]"), pd.Series, np.complexfloating)
+
+    check(assert_type(left.sub(b), "pd.Series[float]"), pd.Series, np.floating)
+    check(assert_type(left.sub(i), "pd.Series[float]"), pd.Series, np.floating)
+    check(assert_type(left.sub(f), "pd.Series[float]"), pd.Series, np.floating)
+    check(assert_type(left.sub(c), "pd.Series[complex]"), pd.Series, np.complexfloating)
+
+    check(assert_type(left.rsub(b), "pd.Series[float]"), pd.Series, np.floating)
+    check(assert_type(left.rsub(i), "pd.Series[float]"), pd.Series, np.floating)
+    check(assert_type(left.rsub(f), "pd.Series[float]"), pd.Series, np.floating)
+    check(
+        assert_type(left.rsub(c), "pd.Series[complex]"), pd.Series, np.complexfloating
+    )
diff --git a/tests/series/arithmetic/int/test_sub.py b/tests/series/arithmetic/int/test_sub.py
new file mode 100644
index 000000000..57056e939
--- /dev/null
+++ b/tests/series/arithmetic/int/test_sub.py
@@ -0,0 +1,131 @@
+from typing import NoReturn
+
+import numpy as np
+from numpy import typing as npt  # noqa: F401
+import pandas as pd
+from typing_extensions import assert_type
+
+from tests import check
+
+left = pd.Series([1, 2, 3])  # left operand
+
+
+def test_sub_py_scalar() -> None:
+    """Test pd.Series[int] - Python native scalars"""
+    b, i, f, c = True, 1, 1.0, 1j
+
+    check(assert_type(left - b, "pd.Series[int]"), pd.Series, np.integer)
+    check(assert_type(left - i, "pd.Series[int]"), pd.Series, np.integer)
+    check(assert_type(left - f, "pd.Series[float]"), pd.Series, np.floating)
+    check(assert_type(left - c, "pd.Series[complex]"), pd.Series, np.complexfloating)
+
+    check(assert_type(b - left, "pd.Series[int]"), pd.Series, np.integer)
+    check(assert_type(i - left, "pd.Series[int]"), pd.Series, np.integer)
+    check(assert_type(f - left, "pd.Series[float]"), pd.Series, np.floating)
+    check(assert_type(c - left, "pd.Series[complex]"), pd.Series, np.complexfloating)
+
+    check(assert_type(left.sub(b), "pd.Series[int]"), pd.Series, np.integer)
+    check(assert_type(left.sub(i), "pd.Series[int]"), pd.Series, np.integer)
+    check(assert_type(left.sub(f), "pd.Series[float]"), pd.Series, np.floating)
+    check(assert_type(left.sub(c), "pd.Series[complex]"), pd.Series, np.complexfloating)
+
+    check(assert_type(left.rsub(b), "pd.Series[int]"), pd.Series, np.integer)
+    check(assert_type(left.rsub(i), "pd.Series[int]"), pd.Series, np.integer)
+    check(assert_type(left.rsub(f), "pd.Series[float]"), pd.Series, np.floating)
+    check(
+        assert_type(left.rsub(c), "pd.Series[complex]"), pd.Series, np.complexfloating
+    )
+
+
+def test_sub_py_sequence() -> None:
+    """Test pd.Series[int] - Python native sequence"""
+    b, i, f, c = [True, False, True], [2, 3, 5], [1.0, 2.0, 3.0], [1j, 1j, 4j]
+
+    check(assert_type(left - b, "pd.Series[int]"), pd.Series, np.integer)
+    check(assert_type(left - i, "pd.Series[int]"), pd.Series, np.integer)
+    check(assert_type(left - f, "pd.Series[float]"), pd.Series, np.floating)
+    check(assert_type(left - c, "pd.Series[complex]"), pd.Series, np.complexfloating)
+
+    check(assert_type(b - left, "pd.Series[int]"), pd.Series, np.integer)
+    check(assert_type(i - left, "pd.Series[int]"), pd.Series, np.integer)
+    check(assert_type(f - left, "pd.Series[float]"), pd.Series, np.floating)
+    check(assert_type(c - left, "pd.Series[complex]"), pd.Series, np.complexfloating)
+
+    check(assert_type(left.sub(b), "pd.Series[int]"), pd.Series, np.integer)
+    check(assert_type(left.sub(i), "pd.Series[int]"), pd.Series, np.integer)
+    check(assert_type(left.sub(f), "pd.Series[float]"), pd.Series, np.floating)
+    check(assert_type(left.sub(c), "pd.Series[complex]"), pd.Series, np.complexfloating)
+
+    check(assert_type(left.rsub(b), "pd.Series[int]"), pd.Series, np.integer)
+    check(assert_type(left.rsub(i), "pd.Series[int]"), pd.Series, np.integer)
+    check(assert_type(left.rsub(f), "pd.Series[float]"), pd.Series, np.floating)
+    check(
+        assert_type(left.rsub(c), "pd.Series[complex]"), pd.Series, np.complexfloating
+    )
+
+
+def test_sub_numpy_array() -> None:
+    """Test pd.Series[int] - numpy array"""
+    b = np.array([True, False, True], np.bool_)
+    i = np.array([2, 3, 5], np.int64)
+    f = np.array([1.0, 2.0, 3.0], np.float64)
+    c = np.array([1.1j, 2.2j, 4.1j], np.complex128)
+
+    check(assert_type(left - b, "pd.Series[int]"), pd.Series, np.integer)
+    check(assert_type(left - i, "pd.Series[int]"), pd.Series, np.integer)
+    check(assert_type(left - f, "pd.Series[float]"), pd.Series, np.floating)
+    check(assert_type(left - c, "pd.Series[complex]"), pd.Series, np.complexfloating)
+
+    # `numpy` typing gives the corresponding `ndarray`s in the static type
+    # checking, where our `__rsub__` cannot override. At runtime, they return
+    # `Series`s with the correct element type.
+    check(assert_type(b - left, NoReturn), pd.Series, np.integer)
+    check(assert_type(i - left, "npt.NDArray[np.int64]"), pd.Series, np.integer)
+    check(assert_type(f - left, "npt.NDArray[np.float64]"), pd.Series, np.floating)
+    check(
+        assert_type(c - left, "npt.NDArray[np.complex128]"),
+        pd.Series,
+        np.complexfloating,
+    )
+
+    check(assert_type(left.sub(b), "pd.Series[int]"), pd.Series, np.integer)
+    check(assert_type(left.sub(i), "pd.Series[int]"), pd.Series, np.integer)
+    check(assert_type(left.sub(f), "pd.Series[float]"), pd.Series, np.floating)
+    check(assert_type(left.sub(c), "pd.Series[complex]"), pd.Series, np.complexfloating)
+
+    check(assert_type(left.rsub(b), "pd.Series[int]"), pd.Series, np.integer)
+    check(assert_type(left.rsub(i), "pd.Series[int]"), pd.Series, np.integer)
+    check(assert_type(left.rsub(f), "pd.Series[float]"), pd.Series, np.floating)
+    check(
+        assert_type(left.rsub(c), "pd.Series[complex]"), pd.Series, np.complexfloating
+    )
+
+
+def test_sub_pd_series() -> None:
+    """Test pd.Series[int] - pandas series"""
+    b = pd.Series([True, False, True])
+    i = pd.Series([2, 3, 5])
+    f = pd.Series([1.0, 2.0, 3.0])
+    c = pd.Series([1.1j, 2.2j, 4.1j])
+
+    check(assert_type(left - b, "pd.Series[int]"), pd.Series, np.integer)
+    check(assert_type(left - i, "pd.Series[int]"), pd.Series, np.integer)
+    check(assert_type(left - f, "pd.Series[float]"), pd.Series, np.floating)
+    check(assert_type(left - c, "pd.Series[complex]"), pd.Series, np.complexfloating)
+
+    check(assert_type(b - left, "pd.Series[int]"), pd.Series, np.integer)
+    check(assert_type(i - left, "pd.Series[int]"), pd.Series, np.integer)
+    check(assert_type(f - left, "pd.Series[float]"), pd.Series, np.floating)
+    check(assert_type(c - left, "pd.Series[complex]"), pd.Series, np.complexfloating)
+
+    check(assert_type(left.sub(b), "pd.Series[int]"), pd.Series, np.integer)
+    check(assert_type(left.sub(i), "pd.Series[int]"), pd.Series, np.integer)
+    check(assert_type(left.sub(f), "pd.Series[float]"), pd.Series, np.floating)
+    check(assert_type(left.sub(c), "pd.Series[complex]"), pd.Series, np.complexfloating)
+
+    check(assert_type(left.rsub(b), "pd.Series[int]"), pd.Series, np.integer)
+    check(assert_type(left.rsub(i), "pd.Series[int]"), pd.Series, np.integer)
+    check(assert_type(left.rsub(f), "pd.Series[float]"), pd.Series, np.floating)
+    check(
+        assert_type(left.rsub(c), "pd.Series[complex]"), pd.Series, np.complexfloating
+    )
diff --git a/tests/series/arithmetic/test_sub.py b/tests/series/arithmetic/test_sub.py
new file mode 100644
index 000000000..636aa242c
--- /dev/null
+++ b/tests/series/arithmetic/test_sub.py
@@ -0,0 +1,192 @@
+from datetime import (
+    datetime,
+    timedelta,
+)
+from typing import (
+    TYPE_CHECKING,
+    NoReturn,
+)
+
+import numpy as np
+from numpy import typing as npt  # noqa: F401
+import pandas as pd
+from typing_extensions import assert_type
+
+from tests import check
+
+if TYPE_CHECKING:
+    from pandas.core.series import TimedeltaSeries  # noqa: F401
+
+left_i = pd.DataFrame({"a": [1, 2, 3]})["a"]  # left operand
+
+
+def test_sub_py_scalar() -> None:
+    """Test pd.Series[Any] - Python native scalars"""
+    b, i, f, c = True, 1, 1.0, 1j
+
+    check(assert_type(left_i - b, pd.Series), pd.Series)
+    check(assert_type(left_i - i, pd.Series), pd.Series)
+    check(assert_type(left_i - f, pd.Series), pd.Series)
+    check(assert_type(left_i - c, pd.Series), pd.Series)
+
+    check(assert_type(b - left_i, pd.Series), pd.Series)
+    check(assert_type(i - left_i, pd.Series), pd.Series)
+    check(assert_type(f - left_i, pd.Series), pd.Series)
+    check(assert_type(c - left_i, pd.Series), pd.Series)
+
+    check(assert_type(left_i.sub(b), pd.Series), pd.Series)
+    check(assert_type(left_i.sub(i), pd.Series), pd.Series)
+    check(assert_type(left_i.sub(f), pd.Series), pd.Series)
+    check(assert_type(left_i.sub(c), pd.Series), pd.Series)
+
+    check(assert_type(left_i.rsub(b), pd.Series), pd.Series)
+    check(assert_type(left_i.rsub(i), pd.Series), pd.Series)
+    check(assert_type(left_i.rsub(f), pd.Series), pd.Series)
+    check(assert_type(left_i.rsub(c), pd.Series), pd.Series)
+
+
+def test_sub_py_sequence() -> None:
+    """Test pd.Series[Any] - Python native sequence"""
+    b, i, f, c = [True, False, True], [2, 3, 5], [1.0, 2.0, 3.0], [1j, 1j, 4j]
+
+    check(assert_type(left_i - b, pd.Series), pd.Series)
+    check(assert_type(left_i - i, pd.Series), pd.Series)
+    check(assert_type(left_i - f, pd.Series), pd.Series)
+    check(assert_type(left_i - c, pd.Series), pd.Series)
+
+    check(assert_type(b - left_i, pd.Series), pd.Series)
+    check(assert_type(i - left_i, pd.Series), pd.Series)
+    check(assert_type(f - left_i, pd.Series), pd.Series)
+    check(assert_type(c - left_i, pd.Series), pd.Series)
+
+    check(assert_type(left_i.sub(b), pd.Series), pd.Series)
+    check(assert_type(left_i.sub(i), pd.Series), pd.Series)
+    check(assert_type(left_i.sub(f), pd.Series), pd.Series)
+    check(assert_type(left_i.sub(c), pd.Series), pd.Series)
+
+    check(assert_type(left_i.rsub(b), pd.Series), pd.Series)
+    check(assert_type(left_i.rsub(i), pd.Series), pd.Series)
+    check(assert_type(left_i.rsub(f), pd.Series), pd.Series)
+    check(assert_type(left_i.rsub(c), pd.Series), pd.Series)
+
+
+def test_sub_numpy_array() -> None:
+    """Test pd.Series[Any] - numpy array"""
+    b = np.array([True, False, True], np.bool_)
+    i = np.array([2, 3, 5], np.int64)
+    f = np.array([1.0, 2.0, 3.0], np.float64)
+    c = np.array([1.1j, 2.2j, 4.1j], np.complex128)
+
+    check(assert_type(left_i - b, pd.Series), pd.Series)
+    check(assert_type(left_i - i, pd.Series), pd.Series)
+    check(assert_type(left_i - f, pd.Series), pd.Series)
+    check(assert_type(left_i - c, pd.Series), pd.Series)
+
+    # `numpy` typing gives the corresponding `ndarray`s in the static type
+    # checking, where our `__rsub__` cannot override. At runtime, they return
+    # `Series`s.
+    # `mypy` thinks the return types are `Any`, which is a bug.
+    check(assert_type(b - left_i, NoReturn), pd.Series)  # type: ignore[assert-type]
+    check(
+        assert_type(i - left_i, "npt.NDArray[np.int64]"), pd.Series  # type: ignore[assert-type]
+    )
+    check(
+        assert_type(f - left_i, "npt.NDArray[np.float64]"), pd.Series  # type: ignore[assert-type]
+    )
+    check(
+        assert_type(c - left_i, "npt.NDArray[np.complex128]"), pd.Series  # type: ignore[assert-type]
+    )
+
+    check(assert_type(left_i.sub(b), pd.Series), pd.Series)
+    check(assert_type(left_i.sub(i), pd.Series), pd.Series)
+    check(assert_type(left_i.sub(f), pd.Series), pd.Series)
+    check(assert_type(left_i.sub(c), pd.Series), pd.Series)
+
+    check(assert_type(left_i.rsub(b), pd.Series), pd.Series)
+    check(assert_type(left_i.rsub(i), pd.Series), pd.Series)
+    check(assert_type(left_i.rsub(f), pd.Series), pd.Series)
+    check(assert_type(left_i.rsub(c), pd.Series), pd.Series)
+
+
+def test_sub_pd_series() -> None:
+    """Test pd.Series[Any] - pandas series"""
+    b = pd.Series([True, False, True])
+    i = pd.Series([2, 3, 5])
+    f = pd.Series([1.0, 2.0, 3.0])
+    c = pd.Series([1.1j, 2.2j, 4.1j])
+
+    check(assert_type(left_i - b, pd.Series), pd.Series)
+    check(assert_type(left_i - i, pd.Series), pd.Series)
+    check(assert_type(left_i - f, pd.Series), pd.Series)
+    check(assert_type(left_i - c, pd.Series), pd.Series)
+
+    check(assert_type(b - left_i, pd.Series), pd.Series)
+    check(assert_type(i - left_i, pd.Series), pd.Series)
+    check(assert_type(f - left_i, pd.Series), pd.Series)
+    check(assert_type(c - left_i, pd.Series), pd.Series)
+
+    check(assert_type(left_i.sub(b), pd.Series), pd.Series)
+    check(assert_type(left_i.sub(i), pd.Series), pd.Series)
+    check(assert_type(left_i.sub(f), pd.Series), pd.Series)
+    check(assert_type(left_i.sub(c), pd.Series), pd.Series)
+
+    check(assert_type(left_i.rsub(b), pd.Series), pd.Series)
+    check(assert_type(left_i.rsub(i), pd.Series), pd.Series)
+    check(assert_type(left_i.rsub(f), pd.Series), pd.Series)
+    check(assert_type(left_i.rsub(c), pd.Series), pd.Series)
+
+
+anchor = datetime(2025, 8, 18)
+left_ts = pd.DataFrame({"a": [anchor + timedelta(hours=h + 1) for h in range(3)]})["a"]
+
+
+def test_sub_py_datetime() -> None:
+    """Test pd.Series[Any] - Python native datetime"""
+    s = anchor
+
+    check(assert_type(left_ts - s, "TimedeltaSeries"), pd.Series, pd.Timedelta)
+
+    check(assert_type(s - left_ts, "TimedeltaSeries"), pd.Series, pd.Timedelta)
+
+    check(assert_type(left_ts.sub(s), "TimedeltaSeries"), pd.Series, pd.Timedelta)
+
+    check(assert_type(left_ts.rsub(s), "TimedeltaSeries"), pd.Series, pd.Timedelta)
+
+
+def test_sub_numpy_datetime() -> None:
+    """Test pd.Series[Any] - numpy datetime(s)"""
+    s = np.datetime64(anchor)
+    a = np.array([s + np.timedelta64(m, "m") for m in range(3)], np.datetime64)
+
+    check(assert_type(left_ts - s, "TimedeltaSeries"), pd.Series, pd.Timedelta)
+    check(assert_type(left_ts - a, "TimedeltaSeries"), pd.Series, pd.Timedelta)
+
+    check(assert_type(s - left_ts, "TimedeltaSeries"), pd.Series, pd.Timedelta)
+    # `numpy` typing gives the corresponding `ndarray`s in the static type
+    # checking, where our `__rsub__` cannot override. At runtime, they return
+    # `Series`s.
+    check(assert_type(a - left_ts, "npt.NDArray[np.datetime64]"), pd.Series, pd.Timedelta)  # type: ignore[assert-type]
+
+    check(assert_type(left_ts.sub(s), "TimedeltaSeries"), pd.Series, pd.Timedelta)
+    check(assert_type(left_ts.sub(a), "TimedeltaSeries"), pd.Series, pd.Timedelta)
+
+    check(assert_type(left_ts.rsub(s), "TimedeltaSeries"), pd.Series, pd.Timedelta)
+    check(assert_type(left_ts.rsub(a), "TimedeltaSeries"), pd.Series, pd.Timedelta)
+
+
+def test_sub_pd_datetime() -> None:
+    """Test pd.Series[Any] - Pandas datetime(s)"""
+    s = pd.Timestamp(anchor)
+    a = pd.Series([s + pd.Timedelta(minutes=m) for m in range(3)])
+
+    check(assert_type(left_ts - s, "TimedeltaSeries"), pd.Series, pd.Timedelta)
+    check(assert_type(left_ts - a, "TimedeltaSeries"), pd.Series, pd.Timedelta)
+
+    check(assert_type(s - left_ts, "TimedeltaSeries"), pd.Series, pd.Timedelta)
+    check(assert_type(a - left_ts, "TimedeltaSeries"), pd.Series, pd.Timedelta)
+
+    check(assert_type(left_ts.sub(s), "TimedeltaSeries"), pd.Series, pd.Timedelta)
+    check(assert_type(left_ts.sub(a), "TimedeltaSeries"), pd.Series, pd.Timedelta)
+
+    check(assert_type(left_ts.rsub(s), "TimedeltaSeries"), pd.Series, pd.Timedelta)
+    check(assert_type(left_ts.rsub(a), "TimedeltaSeries"), pd.Series, pd.Timedelta)
diff --git a/tests/series/test_series.py b/tests/series/test_series.py
index 4070b4b3d..de18d40ac 100644
--- a/tests/series/test_series.py
+++ b/tests/series/test_series.py
@@ -822,8 +822,7 @@ def test_types_element_wise_arithmetic() -> None:
     check(assert_type(s + s2, "pd.Series[int]"), pd.Series, np.integer)
     check(assert_type(s.add(s2, fill_value=0), "pd.Series[int]"), pd.Series, np.integer)
 
-    # TODO this one below should type pd.Series[int]
-    check(assert_type(s - s2, pd.Series), pd.Series, np.integer)
+    check(assert_type(s - s2, "pd.Series[int]"), pd.Series, np.integer)
     check(assert_type(s.sub(s2, fill_value=0), "pd.Series[int]"), pd.Series, np.integer)
 
     # TODO these two below should type pd.Series[int]
@@ -1647,7 +1646,7 @@ def test_series_mul() -> None:
     sm = s * 4
     check(assert_type(sm, pd.Series), pd.Series)
     ss = s - 4
-    check(assert_type(ss, pd.Series), pd.Series)
+    check(assert_type(ss, "pd.Series[int]"), pd.Series, np.integer)
     sm2 = s * s
     check(assert_type(sm2, pd.Series), pd.Series)
     sp = s + 4