Is there a way to plot exponential decay curves with different decay constants on one graph?

I don't know if this is actually helpful because it's so specific, but this is how I would do it (with ggplot2). First, you need data for the function you want to plot. Take the x for all the values you want to display and apply your function with your coefficients to the data. You need to have data points, not just a function, to plot data.

df_simulated <- data.frame("x" = rep(1:100, 2),
                           "class"= rep(c("DNA", "RNA"), each = 100))
df_simulated$y <- c(1683.7 * exp(-0.103 * 1:100),  # DNA
                    578.7455 * exp(-0.156 * 1:100))  # RNA

However, since I never used the packages you used, I don't know how to extract the values from the models, so I took the values in your example plot. It's important that the "simulated" values for both groups are within one dataframe, and that you have a column which attributes the points to the respective group (RNA or DNA). At least it's easier if you do it like this. Then you need a data frame with your actual observations for the dots. I invented data again:

df_observed <- data.frame("x" = c(12, 13, 25, 26, 50, 51),
                          "y" = c(500, 50, 250, 25, 0, 5),
                          "class" = rep(c("DNA", "RNA"), 3))

Then you can create the plot. With color=class you specify that the data points will be grouped by "class" and will be colored accordingly. ("apple" and "banana" are just dummy words to demonstrate linebreaks)

ggplot() +
  geom_line(data = df_simulated, aes(x = x, y = y, color = class), size = 1, linetype = "dashed") +
  geom_point(data = df_observed, aes(x = x, y = y, color = class), size = 4, pch = 1) +
  annotate("text", x = 50, y = 1250, label = "DNA\napple", color = "tomato", hjust = 0) +
  annotate("text", x = 50, y = 750, label  ="RNA\nbanana", color = "steelblue", hjust = 0) +
  ggtitle(expression(~italic("Styela clava")~"(isolated)")) +
  ylab("COI copies per 1ml") +
  xlab("Time since removal of organisms (hours)") +
  theme_classic() +
  theme(legend.position = "none") +
  scale_color_manual(values = c("DNA" = "tomato", "RNA" = "steelblue"))

This is the output:

First note that R squared is normally used for linear models and not for nonlinear models so the use of this statistic is suspect here; however, below we show it anyways since it seems that is what was asked for. A different goodness of fit measurement that is often used is residual standard error. If fm is the fitted model from nls then sigma(fm) is the residual standard error. Smaller values are more favorable. summary(fm) also reports this value.

For each of df4 and r4 we use lm to get starting values (taking log of both sides we get a model that is linear in log(a) and b), run nls fits and get the coefficients.

Now plot the points and add the fitted lines and legend. (Note that in setting up the graph we use rbind which assumes that df4 and r4 have the same column names, which they do.)

Note that the data provided in the question is much different than that shown in the question's image.

The code below does not need starting values since it uses lm to get them, runs nls and automatically extracts whatever information is needed for the graph.

1) Classic graphics In this alternative no packages are used.

r2 <- function(fm, digits = 3) {
  y <- fitted(fm) + resid(fm) 
  r2 <- 1 - deviance(fm) / sum((y - mean(y))^2)
  if (is.numeric(digits)) r2 <- round(r2, digits)
  r2
}

fo <- copies ~ a * exp(b * hours)  # formula used in nls

# get nls fitted model and coefficients for df4
co_d0 <- coef(lm(log(copies) ~ hours, df4, subset = copies > 0))
fmd <- nls(fo, df4, start = list(a = exp(co_d0[[1]]), b = co_d0[[2]]))
co_d <- round(coef(fmd), 4)

# get nls fitted model and coefficients for r4
co_r0 <- coef(lm(log(copies) ~ hours, r4, subset = copies > 0))
fmr <- nls(fo, r4, start = list(a = exp(co_r0[[1]]), b = co_r0[[2]]))
co_r <- round(coef(fmr), 4)

both <- rbind(cbind(df4, col = "red"), cbind(r4, col = "blue"))
plot(both[1:2], col = both$col,
  xlab = "Time since removal of organisms", ylab = "COI copies per 1ml",
  main = "C)" ~ italic("Styela clava") ~ "(isolated)", adj = 0)
lines(fitted(fmd) ~ hours, df4, col = "red", lty = 2)
lines(fitted(fmr) ~ hours, r4, col = "blue", lty = 2)

legend <- c(bquote(DNA), 
            bquote(y == .(co_d[[1]]) * e ^ {.(co_d[[2]])*x}),
            bquote(R^2 == .(r2(fmd))),
            bquote(),
            bquote(RNA), 
            bquote(y == .(co_r[[1]]) * e ^ {.(co_r[[2]])*x}),
            bquote(R^2 == .(r2(fmr))))
legend("right", legend = as.expression(legend), bty = "n",
  text.col = c("red", "red", "red", NA, "blue", "blue", "blue"))

2) ggplot2 This uses ggplot2 and gridtext. r2, fmd, fmr, co_d and co_r are all taken from (1). We use richtest_grob from gridtext to create a custom grob for the legend and pass it using annotate_custom.

library(gridtext)
library(ggplot2)

txt <- sprintf(
  "<span style='color:red'>DNA
    <br>y = %.3f*e<sup>%.3fx</sup> 
    <br>R<sup>2</sup> = %.3f</span>
    <br><br><span style='color:blue'>RNA
    <br>y = %.3f*e<sup>%.3fx</sup>
    <br>R<sup>2</sup> = %.3f</span>",
  co_d[[1]], co_d[[2]], r2(fmd), co_r[[1]], co_r[[2]], r2(fmr))

both2 <- rbind(cbind(df4, col = "red", fitted = fitted(fmd)),
               cbind(r4, col = "blue", fitted = fitted(fmr)))
ggplot(both2, aes(hours, copies, col = I(col))) +
  geom_point() +
  geom_line(aes(y = fitted), linetype = 2) +
  annotation_custom(richtext_grob(txt, hjust = 0)) +
  theme(legend.position = "none") +
  labs(x = "Time since removal of organisms", y = "COI copies per 1ml") +
  ggtitle(("C)" ~ italic("Styela clava") ~ "(isolated)"))

3) lattice

This uses legend from (1) and both2 from (2). First create a plot for the data points. It will also contain the legend, axes and labels. Then add a layer for the fitted lines. main.settings specifies that the main title should be left justified and bold and is adapted from this page.

library(latticeExtra)

main.settings <- list(par.main.text = list(font = 2, just = "left", 
  x = grid::unit(25, "mm")))

xyplot(copies ~ hours, both2, col = both2$col,
    xlab = "Time since removal of organisms", ylab = "COI copies per 1ml",
    main = "C)" ~ italic("Styela clava") ~ "(isolated)", adj = 0,
    key = list(text = list(as.expression(legend), 
      col = c("red", "red", "red", NA, "blue", "blue", "blue")),
        x = 0.65, y = 0.65, columns = 1),
    par.settings = main.settings) +
  as.layer(xyplot(fitted ~ hours, both2, groups = col, type = "l", lty = 2))