GAM: Plotting gam models of different treatments?

I want to figure out the best way to plot gams, but I am getting confused about the best way to do this for my data. Please note that I am not a statistician, so 'stats/coding for dummies' answers are much appreciated.

Briefly, my experiment looked at how flower production changed over time in response to 4 different flower removal treatments that differed by the timing of removal (Early, Peak, Late, Control/None). My data looks like this, where Trt is treatment, RIL is genetic line and Total is total flower production (response).

> head(data2)
DateNum  Date Tray TrtNum   Trt Treat Plant RIL Primary Total NonPrim Dist Dist_num_tray Dist_num_plant
1        1 43426    1      1 Early    E1     1 128       0     0       0    2            50           66.8
3        3 43428    1      1 Early    E1     1 128       0     0       0    2            50           66.8
5        5 43430    1      1 Early    E1     1 128       0     0       0    2            50           66.8
7        7 43432    1      1 Early    E1     1 128       0     0       0    2            50           66.8
9        9 43434    1      1 Early           1 128       0     0       0    2            50           66.8
11      11 43436    1      1 Early           1 128       0     0       0    2            50           66.8


During the times of flower removal (3 day period), flower count for that entire treatment is 0 for 3 days. I knew this would obviously cause there to be a treatment effect, but I am interested in the response after this removal. So I divided my data into different subsets to compare flowering after damage of one of the treatments to the control (no damage).

Here is an example of a simple gam comparing Control vs Early after damage occurred:

gam1 <- gam(Total ~ Trt * RIL + s(DateNum, k = 9, bs = "fs") + s(Plant, bs = "re"), data=ce1230)


I needed to make this to see the results of anova(gam1) to determine if there is a treatment (Trt) effect. However the plot() function (and some other things I've played around with like plotGAM from vortex package) shows only one curve, I think because it is combining the control and the early data into one curve. What I really want to see is two separate gam curves, one for early and one for control.

Should I just use gam1 for testing if there is a treatment effect, but make separate gams for each treatment, then extract the values (predict() function?) and graph those separately? There must be an easier way, or no?

Thank you! Sorry about the length, I just wanted to make sure I was giving enough info

1. See @GavinSimpson 's answer to a related question. Note how he is nesting using the by term within mgcv::gam.
2. Specify the interaction term using ti smooth. From mgcv documentation: