diff --git a/scipy/stats/_resampling.py b/scipy/stats/_resampling.py
index 7d640fe21efe..2e8ec860cdcd 100644
--- a/scipy/stats/_resampling.py
+++ b/scipy/stats/_resampling.py
@@ -674,6 +674,13 @@ def monte_carlo_test(sample, rvs, statistic, *, vectorized=None,
     null_distribution : ndarray
         The values of the test statistic generated under the null hypothesis.
 
+    References
+    ----------
+
+    .. [1] B. Phipson and G. K. Smyth. "Permutation P-values Should Never Be
+       Zero: Calculating Exact P-values When Permutations Are Randomly Drawn."
+       Statistical Applications in Genetics and Molecular Biology 9.1 (2010).
+
     Examples
     --------
 
@@ -709,7 +716,7 @@ def monte_carlo_test(sample, rvs, statistic, *, vectorized=None,
     >>> print(res.statistic)
     0.12457412450240658
     >>> print(res.pvalue)
-    0.701070107010701
+    0.7012
 
     The probability of obtaining a test statistic less than or equal to the
     observed value under the null hypothesis is ~70%. This is greater than
@@ -762,12 +769,12 @@ def monte_carlo_test(sample, rvs, statistic, *, vectorized=None,
 
     def less(null_distribution, observed):
         cmps = null_distribution <= observed
-        pvalues = cmps.sum(axis=0) / n_resamples
+        pvalues = (cmps.sum(axis=0) + 1) / (n_resamples + 1)  # see [1]
         return pvalues
 
     def greater(null_distribution, observed):
         cmps = null_distribution >= observed
-        pvalues = cmps.sum(axis=0) / n_resamples
+        pvalues = (cmps.sum(axis=0) + 1) / (n_resamples + 1)  # see [1]
         return pvalues
 
     def two_sided(null_distribution, observed):
diff --git a/scipy/stats/tests/test_resampling.py b/scipy/stats/tests/test_resampling.py
index 7829e7c85a18..c786dc3553c8 100644
--- a/scipy/stats/tests/test_resampling.py
+++ b/scipy/stats/tests/test_resampling.py
@@ -763,6 +763,15 @@ def statistic1d(x):
         assert_allclose(res.statistic, expected_stat)
         assert_allclose(res.pvalue, expected_p, atol=2*self.atol)
 
+    def test_p_never_zero(self):
+        # Use biased estimate of p-value to ensure that p-value is never zero
+        # per monte_carlo_test reference [1]
+        rng = np.random.default_rng(2190176673029737545)
+        x = np.zeros(100)
+        res = monte_carlo_test(x, rng.random, np.mean,
+                               vectorized=True, alternative='less')
+        assert res.pvalue == 0.0001
+
 
 class TestPermutationTest: