diff --git a/base/utils/NEWS.md b/base/utils/NEWS.md
index e6effc30ee..c5fd5c1675 100644
--- a/base/utils/NEWS.md
+++ b/base/utils/NEWS.md
@@ -1,5 +1,9 @@
 # PEcAn.utils 1.8.2
 
+## Changed
+
+* `read.output()`: `start.year` and `end.year` now default to `NULL` (read all years) instead of `NA`, and both `NULL` and `NA` are accepted as "read all years" sentinels (#3987).
+
 ## Added
 
 * New function `nc_write_varfiles()` creates text summaries of the variables in all netCDF files in an output directory, writing either one file named `nc_vars.txt` if `output_mode = "combined"` or a separate `[filename].nc.var` alongside each netCDF if `output_mode = "paired"` (#3611).
diff --git a/base/utils/R/read.output.R b/base/utils/R/read.output.R
index d26449297e..db72e7658c 100644
--- a/base/utils/R/read.output.R
+++ b/base/utils/R/read.output.R
@@ -40,7 +40,7 @@
 #'   print a summary of the means of each variable for each year.
 #' @param start.year,end.year first and last year of output to read.
 #'   Specify as a date-time (only the year portion is used) or as a
-#'   four-digit number or string. If `NA`, reads all years found in
+#'   four-digit number or string. If `NULL`, reads all years found in
 #'   `outdir`.
 #' @return If `dataframe = FALSE`, a vector of output variables. If
 #'   `dataframe = TRUE`, a `data.frame` of output variables with
@@ -50,8 +50,8 @@
 #' @export
 #' @author Michael Dietze, David LeBauer, Alexey Shiklomanov
 read.output <- function(runid, outdir,
-                        start.year = NA,
-                        end.year = NA,
+                        start.year = NULL,
+                        end.year = NULL,
                         variables = "GPP",
                         dataframe = FALSE,
                         pft.name = NULL,
@@ -63,6 +63,11 @@ read.output <- function(runid, outdir,
   ## cflux = c('GPP', 'NPP', 'NEE', 'TotalResp', 'AutoResp', 'HeteroResp', 'DOC_flux', 'Fire_flux') # kgC m-2 s-1
   ## wflux = c('Evap', 'TVeg', 'Qs', 'Qsb', 'Rainf') # kgH20 m-2 d-1
 
+  # Accept NULL as a synonym for NA for start.year and end.year,
+  # consistent with how other params (variables, ncfiles, pft.name) signal "use all"
+  if (is.null(start.year)) start.year <- NA
+  if (is.null(end.year)) end.year <- NA
+
   if ((missing(runid) || missing(outdir)) && is.null(ncfiles)) {
     PEcAn.logger::logger.severe(
       "`runid` or `outdir` is missing, and `ncfiles` is NULL.",
diff --git a/base/utils/man/read.output.Rd b/base/utils/man/read.output.Rd
index efc77553a2..00deaf428f 100644
--- a/base/utils/man/read.output.Rd
+++ b/base/utils/man/read.output.Rd
@@ -7,8 +7,8 @@
 read.output(
   runid,
   outdir,
-  start.year = NA,
-  end.year = NA,
+  start.year = NULL,
+  end.year = NULL,
   variables = "GPP",
   dataframe = FALSE,
   pft.name = NULL,
@@ -26,7 +26,7 @@ Can be omitted if \code{ncfiles} is set.}
 
 \item{start.year, end.year}{first and last year of output to read.
 Specify as a date-time (only the year portion is used) or as a
-four-digit number or string. If \code{NA}, reads all years found in
+four-digit number or string. If \code{NULL}, reads all years found in
 \code{outdir}.}
 
 \item{variables}{Character vector of variables to be read from
diff --git a/base/utils/tests/testthat/test-read.output.R b/base/utils/tests/testthat/test-read.output.R
index 00a805311c..4fcc8667dc 100644
--- a/base/utils/tests/testthat/test-read.output.R
+++ b/base/utils/tests/testthat/test-read.output.R
@@ -57,6 +57,21 @@ test_that("accepts start and end years as string, number, datetime", {
         expect_length(res_end[["posix"]], 365 * 2)
 })
 
+test_that("accepts NULL as synonym for NA in start.year and end.year", {
+  nulldir <- withr::local_tempdir()
+  times_to_netcdf(0:364, "days since 2001-01-01", nulldir, "2001.nc")
+  times_to_netcdf(0:364, "days since 2002-01-01", nulldir, "2002.nc")
+
+  out_log <- capture.output(type = "message", {
+    res_na   <- read.output(runid = "", outdir = nulldir, variables = "Y",
+                            start.year = NA, end.year = NA, dataframe = TRUE)
+    res_null <- read.output(runid = "", outdir = nulldir, variables = "Y",
+                            start.year = NULL, end.year = NULL, dataframe = TRUE)
+  })
+  expect_equivalent(res_na, res_null)
+  expect_length(res_null$posix, 365 * 2)
+})
+
 test_that("handles arbitrary time offsets", {
 	times_to_netcdf(365:730, "days since 2003-01-01", testdir, "2004.nc")
 	times_to_netcdf( ((0:364)+916) * 24, "hours since 2002-06-30", testdir, "2005.nc")