Which method is correct? (generalized additive model, mgcv) If possible could you please let me know if I am on the right track? I don't know anyone who works with GAMs who I can ask and I would be so appreciative of any help I receive. I use mgcv in R.
My data is daily flower counts ('Total'), with different plant genetic lines ('RIL') and 4 different treatments ('Trt'). I basically want to compare the 4 treatments and see if there is a difference in flower production over time between them. Are the shapes of the flowering curves different? 
An example of my gam: 
gam1 <- gam(Total ~ Trt * RIL + s(DateNum, k = 9, bs = "fs") + s(Plant, bs = "re"), data=ce1230)

I understand that using anova.gam(gam1) I can see if there is a significant difference between treatments and rils. However when I plot gam1, it shows one curve that represents everything. I was hoping to see a plot with separate curves for each of my 4 treatments. Does this mean that I need to make gamtrt1, gamtrt2, gamtrt3, gamtrt4 and somehow compare the separate gams to each other? or use the original anova(gam1) and just plot each treatment separately? Or am I looking at this the entirely wrong way?
 A: I would first try an ANOVA and look at that F-test. If it rejects that, there is evidence that there is variation between treatments. After that, since your outcome if count data, I would do a poisson regression. Unlike GAMs where it can be hard to test the important of effects, Poisson regression gives you your usually p-values which makes it simple to see which things are important to causing the differences in your data. 
https://stats.idre.ucla.edu/r/dae/poisson-regression/
If there is significnet difference between your mean and variance, (poisson says that they should be the same), modify the poisson regression with a Poisson-Gamma model. 
A: I'd guess you want the by parameter for s.  The resulting set of plots then has one for each level of the by factor.  In this example from ?gam.models, there are three levels of fac, so there are three plots for x2 and one for x0.
dat <- gamSim(4)
b <- gam(y ~ fac + s(x2, by=fac) + s(x0), data=dat)
plot(b, pages=1)


A: You have some problems with your data; you can't have a random effect of Plant as a factor and get the plot you showed in the link in comment to @Aaron's Answer. Are you sure that Plant is a factor with more than 1 level? If not, you need to code your data correctly to get a random intercept per plant. Also, can you include both RIL and Plant level effects? Once you've accounted for the separate Plant effects (intercepts), won't that also logically account for the genetic line effects also?
Second, if you are using bs = 'fs', you need to pass in a continuous variable and a factor; so far you only pass in DateNum. At the moment you have
s(DateNum, k = 9, bs = "fs")

and I think you wanted
s(DateNum, Trt, k = 9, bs = "fs")

This model is similar to the one proposed by @Aaron, but is somewhat different in detail. The full model might be
gam(Total ~ Trt * RIL + s(DateNum, Trt, k = 9, bs = "fs") + s(Plant, bs = "re"),
    data=ce1230)

but I suspect even that is wrong? (For fs smooths you don't need the parametric Trt.) So I think,
gam(Total ~ RIL + Trt:RIL + s(DateNum, Trt, k = 9, bs = "fs") + 
    s(Plant, bs = "re"), data=ce1230)

Where the main effect of Trt is actually contained in the fs smooth, so we don't specify it parametrically. 
The main difference between this model and @Aaron's in that here, there is a single smoothness parameter for the smooths of DateNum by Trt, whereas in @Aaron's answer each of the smooths gets it's own smoothness parameter. This boils down to a choice between whether you expect similar wiggliness (the shapes of the smooths can be different) for each smooth or whether you expect some of the four smooths to be a lot wigglier than others? 
What you seem to want is interactions between RIL and Trt plus separate smooths for each Trt. But do you want separate smooths for each RIL and Trt combination? That would require a separate variable formed by the combinations of RIL and Trt available in your data:
ce1230 <- transform(ce1230, 
                    RILTrt = interaction(RIL, Trt, drop = TRUE))

And then you could fit
gam(Total ~ s(DateNum, RILTrt, k = 9, bs = "fs") + 
    s(Plant, bs = "re"), data=ce1230)

or 
gam(Total ~ s(DateNum, k = 9, by = RILTrt) + 
    s(Plant, bs = "re"), data=ce1230)

depending on whether you wanted similar wigliness (use bs = 'fs') or different wigglinesses (use by) for each estimated smooth of DateNum.
The model also needs to respect the non-negative nature of the response; you have counts and you can't have negative counts of anything. Using family = poisson or family = nb would be reasonable starting points.
