diff --git a/qiskit/algorithms/gradients/base_estimator_gradient.py b/qiskit/algorithms/gradients/base_estimator_gradient.py
index 7a2819c8b566..fe71441b45d0 100644
--- a/qiskit/algorithms/gradients/base_estimator_gradient.py
+++ b/qiskit/algorithms/gradients/base_estimator_gradient.py
@@ -23,6 +23,7 @@
 from qiskit.circuit import Parameter, QuantumCircuit
 from qiskit.opflow import PauliSumOp
 from qiskit.primitives import BaseEstimator
+from qiskit.providers import Options
 from qiskit.algorithms import AlgorithmJob
 from qiskit.quantum_info.operators.base_operator import BaseOperator
 
@@ -35,7 +36,7 @@ class BaseEstimatorGradient(ABC):
     def __init__(
         self,
         estimator: BaseEstimator,
-        **run_options,
+        run_options: dict | None = None,
     ):
         """
         Args:
@@ -45,7 +46,9 @@ def __init__(
                 setting. Higher priority setting overrides lower priority setting.
         """
         self._estimator: BaseEstimator = estimator
-        self._default_run_options = run_options
+        self._default_run_options = Options()
+        if run_options is not None:
+            self._default_run_options.update_options(**run_options)
 
     def run(
         self,
@@ -86,10 +89,9 @@ def run(
         # The priority of run option is as follows:
         # run_options in ``run`` method > gradient's default run_options > primitive's default setting.
         run_opts = copy(self._default_run_options)
-        run_opts.update(run_options)
-
+        run_opts.update_options(**run_options)
         job = AlgorithmJob(
-            self._run, circuits, observables, parameter_values, parameters, **run_opts
+            self._run, circuits, observables, parameter_values, parameters, **run_opts.__dict__
         )
         job.submit()
         return job
@@ -163,3 +165,16 @@ def _validate_arguments(
                     f"not match the number of qubits of the {i}-th observable "
                     f"({observable.num_qubits})."
                 )
+
+    def _get_local_run_options(self, run_options: dict) -> Options:
+        """Update the run options in the results.
+
+        Args:
+            run_options: The run options to update.
+
+        Returns:
+            The updated run options.
+        """
+        run_opts = copy(self._estimator.run_options)
+        run_opts.update_options(**run_options)
+        return run_opts
diff --git a/qiskit/algorithms/gradients/base_sampler_gradient.py b/qiskit/algorithms/gradients/base_sampler_gradient.py
index f6cc4611096e..91e3ca23f8fa 100644
--- a/qiskit/algorithms/gradients/base_sampler_gradient.py
+++ b/qiskit/algorithms/gradients/base_sampler_gradient.py
@@ -22,6 +22,7 @@
 
 from qiskit.circuit import QuantumCircuit, Parameter
 from qiskit.primitives import BaseSampler
+from qiskit.providers import Options
 from qiskit.algorithms import AlgorithmJob
 from .sampler_gradient_result import SamplerGradientResult
 
@@ -29,7 +30,7 @@
 class BaseSamplerGradient(ABC):
     """Base class for a ``SamplerGradient`` to compute the gradients of the sampling probability."""
 
-    def __init__(self, sampler: BaseSampler, **run_options):
+    def __init__(self, sampler: BaseSampler, run_options: dict | None = None):
         """
         Args:
             sampler: The sampler used to compute the gradients.
@@ -38,7 +39,9 @@ def __init__(self, sampler: BaseSampler, **run_options):
                 setting. Higher priority setting overrides lower priority setting.
         """
         self._sampler: BaseSampler = sampler
-        self._default_run_options = run_options
+        self._default_run_options = Options()
+        if run_options is not None:
+            self._default_run_options.update_options(**run_options)
 
     def run(
         self,
@@ -77,8 +80,8 @@ def run(
         # The priority of run option is as follows:
         # run_options in `run` method > gradient's default run_options > primitive's default run_options.
         run_opts = copy(self._default_run_options)
-        run_opts.update(run_options)
-        job = AlgorithmJob(self._run, circuits, parameter_values, parameters, **run_opts)
+        run_opts.update_options(**run_options)
+        job = AlgorithmJob(self._run, circuits, parameter_values, parameters, **run_opts.__dict__)
         job.submit()
         return job
 
@@ -134,3 +137,16 @@ def _validate_arguments(
                     f"The number of values ({len(parameter_value)}) does not match "
                     f"the number of parameters ({circuit.num_parameters}) for the {i}-th circuit."
                 )
+
+    def _get_local_run_options(self, run_options: dict) -> dict:
+        """Update the run options in the results.
+
+        Args:
+            run_options: The run options to update.
+
+        Returns:
+            The updated run options.
+        """
+        run_opts = copy(self._sampler.run_options)
+        run_opts.update_options(**run_options)
+        return run_opts
diff --git a/qiskit/algorithms/gradients/estimator_gradient_result.py b/qiskit/algorithms/gradients/estimator_gradient_result.py
index 10c6d74555f6..910e64e69122 100644
--- a/qiskit/algorithms/gradients/estimator_gradient_result.py
+++ b/qiskit/algorithms/gradients/estimator_gradient_result.py
@@ -20,6 +20,8 @@
 
 import numpy as np
 
+from qiskit.providers import Options
+
 
 @dataclass(frozen=True)
 class EstimatorGradientResult:
@@ -29,6 +31,5 @@ class EstimatorGradientResult:
     """The gradients of the expectation values."""
     metadata: list[dict[str, Any]]
     """Additional information about the job."""
-    run_options: dict[str, Any]
-    """run_options for the estimator. Currently, estimator's default run_options is not
-        included."""
+    run_options: Options
+    """run_options for the job."""
diff --git a/qiskit/algorithms/gradients/finite_diff_estimator_gradient.py b/qiskit/algorithms/gradients/finite_diff_estimator_gradient.py
index ecb72288e4bc..a2a52446e4ae 100644
--- a/qiskit/algorithms/gradients/finite_diff_estimator_gradient.py
+++ b/qiskit/algorithms/gradients/finite_diff_estimator_gradient.py
@@ -75,7 +75,6 @@ def _run(
             plus = parameter_values_ + self._epsilon * offset
             minus = parameter_values_ - self._epsilon * offset
             n = 2 * len(indices)
-
             job = self._estimator.run(
                 [circuit] * n, [observable] * n, plus.tolist() + minus.tolist(), **run_options
             )
@@ -92,8 +91,5 @@ def _run(
             n = len(result.values) // 2  # is always a multiple of 2
             gradient_ = (result.values[:n] - result.values[n:]) / (2 * self._epsilon)
             gradients.append(gradient_)
-
-        # TODO: include primitive's run_options as well
-        return EstimatorGradientResult(
-            gradients=gradients, metadata=metadata_, run_options=run_options
-        )
+        run_opt = self._get_local_run_options(run_options)
+        return EstimatorGradientResult(gradients=gradients, metadata=metadata_, run_options=run_opt)
diff --git a/qiskit/algorithms/gradients/finite_diff_sampler_gradient.py b/qiskit/algorithms/gradients/finite_diff_sampler_gradient.py
index 9ad28cfc8a09..bdf79cbf382c 100644
--- a/qiskit/algorithms/gradients/finite_diff_sampler_gradient.py
+++ b/qiskit/algorithms/gradients/finite_diff_sampler_gradient.py
@@ -92,7 +92,5 @@ def _run(
                 gradient_.append(dict(enumerate(grad_dist)))
             gradients.append(gradient_)
 
-        # TODO: include primitive's run_options as well
-        return SamplerGradientResult(
-            gradients=gradients, metadata=metadata_, run_options=run_options
-        )
+        run_opt = self._get_local_run_options(run_options)
+        return SamplerGradientResult(gradients=gradients, metadata=metadata_, run_options=run_opt)
diff --git a/qiskit/algorithms/gradients/lin_comb_estimator_gradient.py b/qiskit/algorithms/gradients/lin_comb_estimator_gradient.py
index 72d4b9c59cd5..357db3c602bc 100644
--- a/qiskit/algorithms/gradients/lin_comb_estimator_gradient.py
+++ b/qiskit/algorithms/gradients/lin_comb_estimator_gradient.py
@@ -129,7 +129,5 @@ def _run(
                 gradient_[idx] += coeff * grad_
             gradients.append(gradient_)
 
-        # TODO: include primitive's run_options as well
-        return EstimatorGradientResult(
-            gradients=gradients, metadata=metadata_, run_options=run_options
-        )
+        run_opt = self._get_local_run_options(run_options)
+        return EstimatorGradientResult(gradients=gradients, metadata=metadata_, run_options=run_opt)
diff --git a/qiskit/algorithms/gradients/lin_comb_sampler_gradient.py b/qiskit/algorithms/gradients/lin_comb_sampler_gradient.py
index 56a7953d05c8..440e083dda50 100644
--- a/qiskit/algorithms/gradients/lin_comb_sampler_gradient.py
+++ b/qiskit/algorithms/gradients/lin_comb_sampler_gradient.py
@@ -120,7 +120,5 @@ def _run(
                 gradient_.append(dict(enumerate(grad_dist)))
             gradients.append(gradient_)
 
-        # TODO: include primitive's run_options as well
-        return SamplerGradientResult(
-            gradients=gradients, metadata=metadata_, run_options=run_options
-        )
+        run_opt = self._get_local_run_options(run_options)
+        return SamplerGradientResult(gradients=gradients, metadata=metadata_, run_options=run_opt)
diff --git a/qiskit/algorithms/gradients/param_shift_estimator_gradient.py b/qiskit/algorithms/gradients/param_shift_estimator_gradient.py
index 15eebf7f9987..e1a6d300385b 100644
--- a/qiskit/algorithms/gradients/param_shift_estimator_gradient.py
+++ b/qiskit/algorithms/gradients/param_shift_estimator_gradient.py
@@ -110,7 +110,5 @@ def _run(
                 values[idx] += coeff * grad_
             gradients.append(values)
 
-        # TODO: include primitive's run_options as well
-        return EstimatorGradientResult(
-            gradients=gradients, metadata=metadata_, run_options=run_options
-        )
+        run_opt = self._get_local_run_options(run_options)
+        return EstimatorGradientResult(gradients=gradients, metadata=metadata_, run_options=run_opt)
diff --git a/qiskit/algorithms/gradients/param_shift_sampler_gradient.py b/qiskit/algorithms/gradients/param_shift_sampler_gradient.py
index 6265fc6f61b2..1f0e254b2b0f 100644
--- a/qiskit/algorithms/gradients/param_shift_sampler_gradient.py
+++ b/qiskit/algorithms/gradients/param_shift_sampler_gradient.py
@@ -115,7 +115,5 @@ def _run(
                 gradient_.append(dict(enumerate(grad_dist)))
             gradients.append(gradient_)
 
-        # TODO: include primitive's run_options as well
-        return SamplerGradientResult(
-            gradients=gradients, metadata=metadata_, run_options=run_options
-        )
+        run_opt = self._get_local_run_options(run_options)
+        return SamplerGradientResult(gradients=gradients, metadata=metadata_, run_options=run_opt)
diff --git a/qiskit/algorithms/gradients/sampler_gradient_result.py b/qiskit/algorithms/gradients/sampler_gradient_result.py
index 8cbb0c9ea205..89b94bf10bf7 100644
--- a/qiskit/algorithms/gradients/sampler_gradient_result.py
+++ b/qiskit/algorithms/gradients/sampler_gradient_result.py
@@ -18,6 +18,8 @@
 from typing import Any
 from dataclasses import dataclass
 
+from qiskit.providers import Options
+
 
 @dataclass(frozen=True)
 class SamplerGradientResult:
@@ -27,5 +29,5 @@ class SamplerGradientResult:
     """The gradients of the sample probabilities."""
     metadata: list[dict[str, Any]]
     """Additional information about the job."""
-    run_options: dict[str, Any]
-    """run_options for the sampler. Currently, sampler's default run_options is not included"""
+    run_options: Options
+    """run_options for the job."""
diff --git a/qiskit/algorithms/gradients/spsa_estimator_gradient.py b/qiskit/algorithms/gradients/spsa_estimator_gradient.py
index 37828c346697..84681dfb3634 100644
--- a/qiskit/algorithms/gradients/spsa_estimator_gradient.py
+++ b/qiskit/algorithms/gradients/spsa_estimator_gradient.py
@@ -119,7 +119,5 @@ def _run(
             indices = [circuits[i].parameters.data.index(p) for p in metadata_[i]["parameters"]]
             gradients.append(gradient[indices])
 
-        # TODO: include primitive's run_options as well
-        return EstimatorGradientResult(
-            gradients=gradients, metadata=metadata_, run_options=run_options
-        )
+        run_opt = self._get_local_run_options(run_options)
+        return EstimatorGradientResult(gradients=gradients, metadata=metadata_, run_options=run_opt)
diff --git a/qiskit/algorithms/gradients/spsa_sampler_gradient.py b/qiskit/algorithms/gradients/spsa_sampler_gradient.py
index 578872db2554..b5426b0b3221 100644
--- a/qiskit/algorithms/gradients/spsa_sampler_gradient.py
+++ b/qiskit/algorithms/gradients/spsa_sampler_gradient.py
@@ -120,7 +120,5 @@ def _run(
                 gradient.append(dict(enumerate(gradient_j)))
             gradients.append(gradient)
 
-        # TODO: include primitive's run_options as well
-        return SamplerGradientResult(
-            gradients=gradients, metadata=metadata_, run_options=run_options
-        )
+        run_opt = self._get_local_run_options(run_options)
+        return SamplerGradientResult(gradients=gradients, metadata=metadata_, run_options=run_opt)
diff --git a/releasenotes/notes/add-gradients-with-primitives-561cf9cf75a7ccb8.yaml b/releasenotes/notes/add-gradients-with-primitives-561cf9cf75a7ccb8.yaml
index 93bd3b5d155a..da0bc9dd6497 100644
--- a/releasenotes/notes/add-gradients-with-primitives-561cf9cf75a7ccb8.yaml
+++ b/releasenotes/notes/add-gradients-with-primitives-561cf9cf75a7ccb8.yaml
@@ -2,9 +2,9 @@
 features:
   - |
     New gradient Algorithms using the primitives have been added. They internally
-    use the primitives to calculate the gradients. There are 3 types of
-    gradient classes (Finite Difference, Parameter Shift, and
-    Linear Combination of Unitary) for a sampler and estimator.
+    use the primitives to calculate the gradients. There are 4 types of
+    gradient classes (Finite Difference, Parameter Shift,
+    Linear Combination of Unitary, and SPSA) for a sampler and estimator.
 
     Example::
     .. code-block:: python
diff --git a/test/python/algorithms/test_estimator_gradient.py b/test/python/algorithms/test_estimator_gradient.py
index a139350e27b5..d859adcf1d30 100644
--- a/test/python/algorithms/test_estimator_gradient.py
+++ b/test/python/algorithms/test_estimator_gradient.py
@@ -374,6 +374,45 @@ def test_gradient_random_parameters(self, grad):
             rtol=1e-4,
         )
 
+    @combine(
+        grad=[
+            FiniteDiffEstimatorGradient,
+            ParamShiftEstimatorGradient,
+            LinCombEstimatorGradient,
+            SPSAEstimatorGradient,
+        ],
+    )
+    def test_run_options(self, grad):
+        """Test estimator gradient's run options"""
+        a = Parameter("a")
+        qc = QuantumCircuit(1)
+        qc.rx(a, 0)
+        op = SparsePauliOp.from_list([("Z", 1)])
+        estimator = Estimator(run_options={"shots": 100})
+        with self.subTest("estimator"):
+            if grad is FiniteDiffEstimatorGradient or grad is SPSAEstimatorGradient:
+                gradient = grad(estimator, epsilon=1e-6)
+            else:
+                gradient = grad(estimator)
+            result = gradient.run([qc], [op], [[1]]).result()
+            self.assertEqual(result.run_options.get("shots"), 100)
+
+        with self.subTest("gradient init"):
+            if grad is FiniteDiffEstimatorGradient or grad is SPSAEstimatorGradient:
+                gradient = grad(estimator, epsilon=1e-6, run_options={"shots": 200})
+            else:
+                gradient = grad(estimator, run_options={"shots": 200})
+            result = gradient.run([qc], [op], [[1]]).result()
+            self.assertEqual(result.run_options.get("shots"), 200)
+
+        with self.subTest("gradient run"):
+            if grad is FiniteDiffEstimatorGradient or grad is SPSAEstimatorGradient:
+                gradient = grad(estimator, epsilon=1e-6, run_options={"shots": 200})
+            else:
+                gradient = grad(estimator, run_options={"shots": 200})
+            result = gradient.run([qc], [op], [[1]], shots=300).result()
+            self.assertEqual(result.run_options.get("shots"), 300)
+
 
 if __name__ == "__main__":
     unittest.main()
diff --git a/test/python/algorithms/test_sampler_gradient.py b/test/python/algorithms/test_sampler_gradient.py
index 51e0246aee6c..e1a4bfef1c30 100644
--- a/test/python/algorithms/test_sampler_gradient.py
+++ b/test/python/algorithms/test_sampler_gradient.py
@@ -507,6 +507,45 @@ def test_gradient_random_parameters(self, grad):
             array2 = _quasi2array(res2, num_qubits)
             np.testing.assert_allclose(array1, array2, rtol=1e-4)
 
+    @combine(
+        grad=[
+            FiniteDiffSamplerGradient,
+            ParamShiftSamplerGradient,
+            LinCombSamplerGradient,
+            SPSASamplerGradient,
+        ],
+    )
+    def test_run_options(self, grad):
+        """Test sampler gradient's run options"""
+        a = Parameter("a")
+        qc = QuantumCircuit(1)
+        qc.rx(a, 0)
+        qc.measure_all()
+        sampler = Sampler(run_options={"shots": 100})
+        with self.subTest("sampler"):
+            if grad is FiniteDiffSamplerGradient or grad is SPSASamplerGradient:
+                gradient = grad(sampler, epsilon=1e-6)
+            else:
+                gradient = grad(sampler)
+            result = gradient.run([qc], [[1]]).result()
+            self.assertEqual(result.run_options.get("shots"), 100)
+
+        with self.subTest("gradient init"):
+            if grad is FiniteDiffSamplerGradient or grad is SPSASamplerGradient:
+                gradient = grad(sampler, epsilon=1e-6, run_options={"shots": 200})
+            else:
+                gradient = grad(sampler, run_options={"shots": 200})
+            result = gradient.run([qc], [[1]]).result()
+            self.assertEqual(result.run_options.get("shots"), 200)
+
+        with self.subTest("gradient run"):
+            if grad is FiniteDiffSamplerGradient or grad is SPSASamplerGradient:
+                gradient = grad(sampler, epsilon=1e-6, run_options={"shots": 200})
+            else:
+                gradient = grad(sampler, run_options={"shots": 200})
+            result = gradient.run([qc], [[1]], shots=300).result()
+            self.assertEqual(result.run_options.get("shots"), 300)
+
 
 def _quasi2array(quasis: List[QuasiDistribution], num_qubits: int) -> np.ndarray:
     ret = np.zeros((len(quasis), 2**num_qubits))