diff --git a/DESCRIPTION b/DESCRIPTION
index 3b673f2..7d14541 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -10,10 +10,13 @@ Authors@R: c(
         person("Alan", "O'Callaghan", role=c("ctb", "cre"), email="alan.ocallaghan@outlook.com"),
         person("Yun", "Peng", role=c("ctb"), email="yunyunp96@gmail.com"),
         person("Leo", "Lahti", role=c("ctb"), email="leo.lahti@utu.fi",
-	             comment = c(ORCID = "0000-0001-5537-637X"))
+	             comment = c(ORCID = "0000-0001-5537-637X")),
+	    person("Tuomas", "Borman", role = c("ctb"),
+	        email = "tuomas.v.borman@utu.fi",
+            comment = c(ORCID = "0000-0002-8563-8884"))
 	)
-Version: 1.31.2
-Date: 2024-01-28
+Version: 1.31.3
+Date: 2024-06-04
 License: GPL-3
 Title: Single-Cell Analysis Toolkit for Gene Expression Data in R
 Description: A collection of tools for doing various analyses of
diff --git a/R/plotExpression.R b/R/plotExpression.R
index a2f3a8c..b7e7902 100644
--- a/R/plotExpression.R
+++ b/R/plotExpression.R
@@ -87,6 +87,10 @@
 #' plotExpression(example_sce, rownames(example_sce)[1:6],
 #'     colour_by = "Mutation_Status", shape_by = "Treatment",
 #'     size_by = "Gene_0010")
+#' 
+#' ## use boxplot as well as violin plot
+#' plotExpression(example_sce, rownames(example_sce)[1:6],
+#'     show_boxplot = TRUE, show_violin = FALSE)
 #'
 #' ## plot expression against expression values for Gene_0004
 #' plotExpression(example_sce, rownames(example_sce)[1:4],
diff --git a/R/plot_central.R b/R/plot_central.R
index 23b9001..9ebee1f 100644
--- a/R/plot_central.R
+++ b/R/plot_central.R
@@ -30,6 +30,7 @@
 #' Defaults to \code{FALSE}.}
 #' \item{\code{show_se}:}{Logical, should standard errors for the fitted line be shown on a scatter plot when \code{show_smooth=TRUE}?
 #' Defaults to \code{TRUE}.}
+#' \item{\code{show_boxplot}:}{Logical, should a box plot be shown? Defaults to \code{FALSE}.}
 #' }
 #'
 #' @section Miscellaneous fields: Addititional fields can be added to the
@@ -63,7 +64,7 @@ NULL
 #' @importFrom ggbeeswarm geom_quasirandom
 #' @importFrom ggplot2 ggplot geom_violin xlab ylab stat_summary geom_jitter
 #'   position_jitter coord_flip geom_point stat_smooth geom_tile theme_bw theme
-#'   geom_bin2d geom_hex stat_summary_2d stat_summary_hex
+#'   geom_bin2d geom_hex stat_summary_2d stat_summary_hex geom_boxplot
 .central_plotter <- function(object, xlab = NULL, ylab = NULL,
                              colour_by = NULL, shape_by = NULL, size_by = NULL, fill_by = NULL,
                              show_median = FALSE, show_violin = TRUE, show_smooth = FALSE, show_se = TRUE,
@@ -71,7 +72,7 @@ NULL
                              theme_size = 10, point_alpha = 0.6, point_size = NULL, point_shape = 19, add_legend = TRUE,
                              point_FUN = NULL, jitter_type = "swarm",
                              rasterise = FALSE, scattermore = FALSE, bins = NULL,
-                             summary_fun = "sum", hex = FALSE)
+                             summary_fun = "sum", hex = FALSE, show_boxplot = FALSE)
 # Internal ggplot-creating function to plot anything that involves points.
 # Creates either a scatter plot, (horizontal) violin plots, or a rectangle plot.
 {
@@ -102,6 +103,28 @@ NULL
                     c(viol_args, list(colour = "gray60", alpha = 0.2, scale = "width", width = 0.8))
                 )
         }
+        # Adding box plot
+        if (show_boxplot) {
+            if (is.null(fill_by)) {
+                box_args <- list(fill="grey90")
+            } else {
+                box_args <- list(mapping=aes(fill=.data[[fill_by]]))
+            }
+            # If violin plot is plotted, make the width of box plot smaller to
+            # improve readability.
+            if (show_violin) {
+                box_args[["width"]] <- 0.25
+            }
+            box_args <- c(box_args, list(colour = "black", alpha = 0.2))
+            # If user wants that jitter plot is not added, add outliers.
+            # Otherwise remove outliers since then they would be plotted twice;
+            # once as outlier and once as part of jitter plot.
+            if (!is.na(point_shape)) {
+                box_args[["outlier.shape"]] <- NA
+            }
+            plot_out <- plot_out +
+                do.call(geom_boxplot, box_args)
+        }
 
         # Adding median, if requested.
         if (show_median) {
diff --git a/man/plotExpression.Rd b/man/plotExpression.Rd
index d10d69d..c33b429 100644
--- a/man/plotExpression.Rd
+++ b/man/plotExpression.Rd
@@ -153,6 +153,10 @@ plotExpression(example_sce, rownames(example_sce)[1:6],
     colour_by = "Mutation_Status", shape_by = "Treatment",
     size_by = "Gene_0010")
 
+## use boxplot instead of violin plot
+plotExpression(example_sce, rownames(example_sce)[1:6],
+    show_boxplot = TRUE, show_violin = FALSE)
+
 ## plot expression against expression values for Gene_0004
 plotExpression(example_sce, rownames(example_sce)[1:4],
     "Gene_0004", show_smooth = TRUE)
diff --git a/man/scater-plot-args.Rd b/man/scater-plot-args.Rd
index a007d16..461e154 100644
--- a/man/scater-plot-args.Rd
+++ b/man/scater-plot-args.Rd
@@ -37,6 +37,7 @@ Defaults to \code{TRUE}.}
 Defaults to \code{FALSE}.}
 \item{\code{show_se}:}{Logical, should standard errors for the fitted line be shown on a scatter plot when \code{show_smooth=TRUE}?
 Defaults to \code{TRUE}.}
+\item{\code{show_boxplot}:}{Logical, should a box plot be shown? Defaults to \code{FALSE}.}
 }
 }