c and d (the variables returned from boxplot) are matrices of handles to the lines created by boxplot. In legend, you can specify exactly which lines to make a legend for, which in this case is certain elements from c and d, along with a and b.
Also, you can modify properties of the boxplot lines (such as LineWidth) by setting properties of the lines in c and d.
bothgrouped0hr = readtable('both_grouped_0hr.xlsx');
ROSdata = readtable('ROS_data.xlsx');
yyaxis left
ticks = 0:2:26;
xticks(ticks);
a = plot(ROSdata.Time,ROSdata.ROSpercell,'LineWidth',1);
a.Color = '#0072BD';
hold on
b = plot(ROSdata.Time,ROSdata.Expulsion,'-','LineWidth',1);
b.Color = '#D95319';
ax = gca;
ax.YAxis(1).Color = 'k';
yl = yline(1.41e-14, 'LineWidth',1); %threshold
yl.Color = [0.9290 0.6940 0.1250];
yyaxis right
c = boxplot(bothgrouped0hr.Turbinaria,bothgrouped0hr.Day,'colors','k','widths',1,'Positions',[1, 13, 18, 26]);
hold on
% here's how you can set the LineWidth:
set(c,'LineWidth',2)
d = boxplot(bothgrouped0hr.Duncan,bothgrouped0hr.Day,'colors','m','widths',1, 'Symbol', '','Positions',[0, 13, 18, 26]);
% here's how you can set the LineWidth:
set(d,'LineWidth',2)
xticks(ticks)
xticklabels(sprintfc('%d',ticks))
ax = gca;
ax.YAxis(2).Color = 'k';
% c and d, for reference.
% Each box (there are 4 boxes) contains 7 lines,
% so 7 rows, 4 columns. With some trial-and-error,
% you can figure out which element is which line.
c
d
legend( ...
[a b c(3,1) d(3,1)], ...
["ROS", "norm expulsion", "Turbinaria", "Duncan"], ...
'location', 'southeast')
